无线通信代码搬运/复现系列(1)

“Revisiting the MIMO Capacity with Per-antenna Power Constraint: Fixed-point Iteration and Alternating Optimization,” IEEE Trans. Wireless Commun., vol. 18, no. 1, pp. 388-401, Jan. 2019 by T. M. Pham, R. Farrell, and L.-N. Tran

摘要— 在本文中，我们重新审视了在每天线功率约束（PAPC）下计算MIMO容量的基本问题。与可能接受类似注水解决方案的和功率约束对应物不同，PAPC的MIMO容量主要是在通用凸优化框架下研究的。这些方法的两个主要缺点是：1）它们的复杂性与问题大小迅速成比例，这对于大规模天线系统没有吸引力;和/或2）它们的收敛特性对问题数据很敏感。作为起点，我们首先考虑单用户MIMO场景，并提出了两种可证明收敛迭代算法来求其容量，第一种方法基于定点迭代，另一种方法基于交替优化和极小最大对偶性。具体而言，与现有方法相比，所提出的两种方法在每次迭代中都能利用填充水算法，收敛速度更快。然后，我们将所提出的解决方案扩展到多用户MIMO系统，这些系统具有基于脏纸编码的传输策略。在这方面，使用PAPC的高斯广播信道的容量区域也是使用闭合形式表达式计算的。数值结果证明了所提出的解决方案优于现有方法。
索引术语 — MIMO、定点迭代、交替优化、最小最大值对偶性、填充水、脏纸编码。

固定点算法

function [Sopt, nIterations, err_seq, obj_seq] = Algorithm1_FixedPoint(H, PAPC, errtol, maxIterations)
% 初始化
[m,n] = size(H);

lambda_tilde = (ones(n,1)); 
lambda = 1./lambda_tilde;
[~, R] = qr(H); % 第2步：进行QR分解
err_seq = zeros(maxIterations,1);
obj_seq = zeros(maxIterations,1);
nIterations = maxIterations;
for iIter=1:maxIterations    
    % 第4步
    [~, Sigma_bar, V] = svd(R*diag(lambda.^(-1/2)));
    Sigma_bar(Sigma_bar<0) = 0;
    Sv = Sigma_bar(Sigma_bar>0).^(-2);
    Svs = 1-Sv;
    Svs(Svs <0) = 0;
    % 第6步
    Phi = V*diag(1-[Svs;zeros(n-length(Svs),1)])*V';
    % 第7步
    S = diag((lambda).^(-1)) - diag((lambda).^(-1/2))*Phi*diag((lambda).^(-1/2));
    obj_seq(iIter) = real(log(det(eye(m)+H*S*H')));
    % 第8步
    err_seq(iIter) = abs(sum((((lambda))).*(diag(S)-PAPC)));
    % 第9步
    lambda_tilde = real(PAPC + diag(Phi).*lambda_tilde);
    % 第10步
    lambda = ((lambda_tilde).^(-1));
    
    % 如果超过误差限则提前终止
    if (err_seq(iIter) < errtol)
        err_seq(iIter+1:maxIterations)=[];
        obj_seq(iIter+1:maxIterations)=[];
        nIterations = iIter;
        break;
    end
end

Sopt = S;
end

交替优化

function [Sopt, nIterations, err_seq, obj_seq] = Algorithm2_AlternatingOptimization(H, PAPC, errtol, maxIterations)
P_sum = sum(PAPC);
%initialization
[~, nTx] = size(H);
[G,R] = qr(H);
q = ones(nTx,1);
err_seq = zeros(maxIterations,1);
obj_seq = zeros(maxIterations,1);
nIterations = maxIterations;
for iIter =1:maxIterations    
    % step 3
    [U_tilde, Sigma_tilde, V_tilde] = svd(R*diag(q.^(-1/2)),'econ');
    Sigma_tilde = diag(Sigma_tilde); 
    [Sigma_tilde, ind] = sort(Sigma_tilde,'descend'); % sort eigen channels for water filling
    U_tilde = U_tilde(:, ind); % rearrange the columns of U after sorting
    V_tilde = V_tilde(:,ind); % rearrange the columns of V after sorting

    zeroeigchan = (Sigma_tilde < 1e-7); % ignore very small eigen channels
    U_tilde(:,zeroeigchan) = []; % remove columns of U accordingly
    V_tilde(:,zeroeigchan) = []; % remove columns of V accordingly
    U = (G*U_tilde); % U from SVD of H*diag(q.^(-1))
    Sigma = Sigma_tilde.^2;
    power = wf(Sigma,P_sum); % perform water filling
    chan_pos = (power>0); % get stricly positive eigen channels
    power_pos = power(chan_pos);

    U(:,~chan_pos) = []; % remove columns of U that have no power
    Sigma = Sigma(chan_pos); % remove eigen channels that have no power
    S_bar = U*diag(power_pos)*U'; % the optimal S_bar for given Q
    % end of step 3

    % calculate the objective
    obj_seq(iIter) = real(log((det(diag(q) + H'*S_bar*H)))) - sum(log(q));
    if (iIter>1)
        err_seq(iIter-1)=abs(obj_seq(iIter)-obj_seq(iIter-1));
        if(err_seq(iIter-1)<errtol)
            err_seq(iIter:end) = [];
            obj_seq(iIter+1:end) = [];
            nIterations = iIter;
            break
        end
    end
    % step 6
    % Newton method for solving (26)
    V_dot = V_tilde(:,chan_pos);
    % step 5
    phi_inv = 1./q - real(diag(diag(q.^-0.5)*V_dot*diag(1./(1+1./(Sigma.*power_pos)))*(V_dot')*diag(q.^-0.5)));
    gamma = 0.01;
    fgamma=1;
    while(abs(fgamma) > errtol)
        fgamma = sum(1./(gamma+phi_inv./PAPC)) - P_sum;
        fgamma_diff = -sum(1./((phi_inv./PAPC + gamma).^2));
        gamma = gamma - fgamma/fgamma_diff;
    end
    q = 1./(phi_inv+ gamma*PAPC);
    % end of newton's method

    
end

[U, ~, V] = svd(H*diag(q.^(-1/2)), 'econ');
Sopt = diag(q.^(-1/2))*V*U'*S_bar*U*V'*diag(q.^(-1/2));
end

function power= wf(eigchan,P)
% water filling algorithm for solving the problem max sum(log(1+eigchan_i.*p_i)) s.t. sum(p) == P
% NOTE: eigen must be shorted in descending order
waterlevel = 1;
nEigchans = length(eigchan);
igamma = waterlevel./(eigchan);
temp = 0;
% water filling algorithm
for k = nEigchans:-1:1
    temp = (P+sum(1./eigchan(1:k))*waterlevel)/k;
    if ((temp-igamma(k))>0)
        break;
    end
end
power = max(temp-igamma,0);
end

主程序

clear ;
clc
rng(1)
%initialization
SNRdB = 0;
P = 10.^(SNRdB/10);
eps = 1e-6;

maxIterations = 50;


nTx = 5;
nRx = 3;
H = (randn(nRx, nTx) + 1i*randn(nRx, nTx))/sqrt(2);
% 
% H =[0.2581+1i*0.6535i  0.2623+1i*0.9434i;
% 0.4385+1i*0.3081 0.4090-1i*0.2288];
% [nRx,nTx]=size(H);

PAPC = (P/nTx)*ones(nTx,1); % equal power constraint
% cvx_solver mosek
cvx_expert true
cvx_begin quiet
variable X(nTx,nTx) complex semidefinite
maximize(log_det(eye(nRx)+H*X*H'))
diag(X) <= PAPC
X == hermitian_semidefinite(nTx)
cvx_end
cvx_optval

%Alg1, fixed point
[Sopt_fp, nIterations_fp, err_seq_fp,obj_seq_fp] = Algorithm1_FixedPoint(H, PAPC, eps, maxIterations);

CMIMO_Alg1 = real(log(det(eye(nRx) + H*Sopt_fp*H')))

%Alg2, alternating optimization
[Sopt_ao, nIterations_ao, err_seq_ao,obj_seq_ao] = Algorithm2_AlternatingOptimization(H, PAPC, eps, maxIterations);
CMIMO_Alg2 = real(log(det(eye(nRx) + H*Sopt_ao*H')))



%plot duality
subplot(2,1,1)
semilogy(1:nIterations_fp,err_seq_fp,'--b','LineWidth',1.5);
hold on
semilogy(1:nIterations_ao-1,err_seq_ao,'-k','LineWidth',1.5);
legend('Algorithm 1', 'Algorithm 2','Location','Best');
xlabel('Iteration Index','FontSize',12,'FontWeight','bold');
ylabel('Duality gap','FontSize',12,'FontWeight','bold');
title('Residual error')
subplot(2,1,2)
plot(obj_seq_fp,'--b')
hold on
plot(obj_seq_ao,'-k')
plot(cvx_optval*ones(length(obj_seq_ao),1),'-r')
title('Convergence of the objective')
legend('Algorithm 1', 'Algorithm 2','Optimal Objective (CVX)','Location','Best');
saveas(gcf,'../results/convergence.png')