## Outer Approximation Algorithm for Sum Rate MaximizationThe 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 FormulationThe 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 CodeOur 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 ## A Tiny Numerical ExampleThe 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. |