Outer Approximation Algorithm for Sum Rate Maximization

The outer approximation algorithm is proposed in a IEEE Journal on Selected Areas in Communications journal in 2011. An overview survey is available in Wireless Network Optimization by Perron-Frobenius Theory.

The Problem Formulation

The weighted sum rate maximization problem in a multiuser Gaussian interference channel subject to affine power constraint can be stated as:

$$ \begin{array}[c]{rl} \mbox{maximize} & \displaystyle\sum_{l=1}^L w_l\log\big(1+\mathsf{SINR}_l(\mathbf{p})\big)\\ \mbox{subject to} & \mathbf{a}_l^\top\mathbf{p}\le\bar{p}_l,\,\;\;l=1,\dots,L,\\ \mbox{variables:} & \mathbf{p}, \end{array} $$

where $\mathbf{w}=(w_1,\dots,w_L)^\top \ge 0$ is a given probability vector, and $w_l$ is a weight assigned to the $l$th link to reflect priority (a larger weight reflects a higher priority). The power budget constraint set is modeled by the nonnegative vectors $\mathbf{a}_l,\;l=1,\dots,L$ and the upper bound $\bar{\mathbf{p}}$.

Let us denote $\boldsymbol\gamma$ as the SINR vector of the users , i.e., $\boldsymbol{\gamma}=(\gamma_1,\dots,\gamma_L)^\top>0$. The weighted sum rate maximization problem is equivalent to the following problem:

$$ \begin{array}[c]{rl} \mbox{maximize} & \displaystyle\sum_{l=1}^L w_l\log\big(1+\gamma_l\big)\\ \mbox{subject to} & \rho\big( \mathop{{\rm diag}}(\boldsymbol{\gamma})(\mathbf{F}+(1/\bar{p}_l)\mathbf{v}\mathbf{a}_l^\top) \big) \le 1,\,\;\;l=1,\dots,L,\\ \mbox{variables:} & \boldsymbol{\gamma}, \end{array} $$

where $\rho(\cdot)$ denotes the Perron-Frobenius eigenvalue function and whose optimal $\boldsymbol{\gamma}$ yields the original optimal $\mathbf{p}$ through a Perron-Frobenius eigenvector relationship. Now, let $\tilde{\boldsymbol{\gamma}}=\log{\boldsymbol{\gamma}}$. Then, the weighted sum rate maximization problem can be further rewritten as:

$$ \begin{array}[c]{rl} \mbox{maximize} & \displaystyle\sum_{l=1}^L w_l\log\big(1+e^{\tilde{\gamma}_l}\big)\\ \mbox{subject to} & \rho\big( \mathop{{\rm diag}}(e^{\tilde{\boldsymbol{\gamma}}})(\mathbf{F}+(1/\bar{p}_l)\mathbf{v}\mathbf{a}_l^\top) \big) \le 1,\,\;\;l=1,\dots,L,\\ \mbox{variables:} & \tilde{\boldsymbol{\gamma}}, \end{array} $$

which, notably, maximizes a convex objective function over a closed unbounded convex set.

The MATLAB Code

Our approach is as follows: The feasible region containing the optimal extreme point is first embedded inside a compact polyhedral convex set (the tightest possible that is ensured by fundamental results in nonnegative matrix theory and the Perron-Frobenius theorem). Infeasible regions are then successively removed from this initial polyhedral set. This method generates a nested sequence of polyhedrons approximating the closed unbounded convex set that yields the global optimal solution $\tilde{\boldsymbol\gamma}^\star$ asymptotically from the exterior.

clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem parameter initialization
% L: number of users.
% G: channel gain.
% F: nonnegative matrix. F_{lj} = G_{lj}/G_{ll} if l \ne j, and F_{lj} = 0 if l = j
% var: noise power.
% v: nonnegative vector. v_l = var_l/G_{ll}
% pmax: upper bound of the weighted power constraints.
% a: weight of the power constraints.
% w: a probability vector that is a weight assigned to links to reflect priority.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

L = 2;

G = zeros(L,L);
F = zeros(L,L);
v = zeros(L,1);

var = ones(L,1);

pmax = 2.*ones(L,1)+2.*rand(L,1);
a = rand(L,L);
w = rand(L,1);
w = w./sum(w);

for l = 1:1:L
    for j = 1:1:L
        if l~=j
            G(l,j)=0.1+0.1*rand(1,1);
        else
            G(l,j)=0.6+0.3*rand(1,1);
        end
    end
end

for l=1:1:L
    for j=1:1:L
        if l ~= j
            F(l,j)=G(l,j)/G(l,l);
        else
            F(l,j)=0;
        end
    end
    v(l) = var(l)/G(l,l);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% max_iteration: maximal number of iterations.
% epsilon: stopping criterion if ||p(k+1)-p(k)|| <= epsilon.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

max_iteration = 100;
epsilon = 0.01;

% Initialize the polyhedral convex set.
% A*gamma_tilde <= b is the polyhedral convex set.
A = []; b = [];
for l = 1:1:L
    B(:,:,l) = F+(1/pmax(l)).*(v*a(:,l)');
    
    [vec lamb] = eig(B(:,:,l));
    [val index]=max(max(lamb));
    x = vec(:,index);
    [vec lamb] = eig(B(:,:,l).');
    [val index]=max(max(lamb));
    y = vec(:,index);
    xy = (x.*y./sum(x.*y));
    A = [A;xy'];
    b = [b;-log(max(eig(B(:,:,l))))];
end
A = [A; -eye(L)]; b=[b; 100*ones(L,1)];

gamma_tilde = rand(L,1);
tolerance = 1;
k = 0;
power_evolution = [];

while tolerance >= epsilon || k <= max_iteration
    
    gamma_tilde_k = gamma_tilde;
    
    % Step 1 in Algorithm 1: compute the vertices.
    [V,nr] = con2vert(A,b);
    
    % Step 2 in Algorithm 1: select max{sumrate} from the vertices.
    sumrate = -Inf;
    s = size(V,1);
    for t = 1:1:s
        temp = sum(w.*log(ones(L,1)+exp(V(t,:)')));
        if temp > sumrate
            sumrate = temp;
            gamma_tilde = V(t,:)';
        end
    end
    
    % Step 3 in Algorithm 1: compute power.
    power = inv( eye(L) - diag(exp(gamma_tilde)) * F ) * diag(exp(gamma_tilde)) * v;
    power_evolution = [power_evolution ; power'];
    
    % Step 4 in Algorithm 1: let J = max{ diag(exp(gamma_tilde))*B(:,:,l) }.
    J=1;
    rho = 0;
    for l = 1:1:L
        Btemp = diag(exp(gamma_tilde))*B(:,:,l);
        if max(eig(Btemp)) > rho
            rho = max(eig(Btemp));
            J = l;
        end
    end
    
    % Step 5 in Algorithm 1: compute the Perron right and left eigenvectors
    % of diag(exp(gamma_tilde))*B(:,:,J).
    Btemp = diag(exp(gamma_tilde))*B(:,:,J);
    [vec lamb] = eig(Btemp);
    [val index]=max(max(lamb));
    x = vec(:,index);
    [vec lamb] = eig(Btemp.');
    [val index]=max(max(lamb));
    y = vec(:,index);
    xy = x.*y./sum(x.*y);
    % Step 6 in Algorithm 1: add a hyperplane to the polyhedral convex set.
    A = [A;xy'];
    b = [b;xy'*gamma_tilde-log(max(eig(Btemp)))];
    
    % Step 7 in Algorithm 1: update the iteration number.
    tolerance = norm(gamma_tilde_k-gamma_tilde);
    k = k+1;
end

% plot the power evolution.
set(gca, 'Fontname', 'Times newman', 'Fontsize', 15);
plot(1:1:k,power_evolution(:,1),'-o',1:1:k,power_evolution(:,2),'-^','linewidth',1.5);
legend('User 1','User 2');
xlim([0 k]);
ylim([min(power)-2 max(power)+4]);

Attention: The function con2vert() used in our code to compute the vertices is from matlab central.

A Tiny Numerical Example

The problem parameters are set randomly in the code above. The following figure, obtained by running the code directly, plots the power evolution for a 2-user case.