I’m solving a robust SDP problem in CVX (MATLAB) for beamforming with three types of streams (stream n, stream c, and stream p). Each stream has its own Schur-complement PSD constraint written as a (M+1)×(M+1) block matrix blk == semidefinite(M+1).
The cvx code is given as follows.
cvx_begin %quiet
variable Q_n(M,M) hermitian semidefinite;
variable Q_c(M,M) hermitian semidefinite;
variable Q_p(M,M,N_user) hermitian semidefinite;
variable common_rate(N_user,1) nonnegative;
variable rho_p(N_user,1) nonnegative;
variable Gamma_n(N_user,1) nonnegative;
variable Gamma_c(N_user,1) nonnegative;
variable Gamma_p(N_user,1) nonnegative;
% lamda introduced by S-Procedure
variable lambda_n(N_user, M) nonnegative;
variable lambda_c(N_user, M) nonnegative;
variable lambda_p(N_user, M) nonnegative;
% power
expression tmp_power;
tmp_power = trace(Q_n) + trace(Q_c);
for kk = 1:N_user
tmp_power = tmp_power + trace(Q_p(:,:,kk));
end
expression Q_p_sum;
Q_p_sum = zeros(M,M);
for j = 1:N_user
Q_p_sum = Q_p_sum + Q_p(:,:,j);
end
minimize ( (1-regulation_parameter) * sum_square( common_rate + rho_p - user_demand ) ...
+ regulation_parameter * tmp_power )
subject to
% power constraint
tmp_power <= transmit_power_W;
% SDP
Q_n == hermitian_semidefinite(M);
Q_c == hermitian_semidefinite(M);
for kk = 1:N_user
Q_p(:,:,kk) == hermitian_semidefinite(M);
end
% SINR/rate constraint
for k = 1:N_user
sum(common_rate) <= log(1 + Gamma_c(k)) / log(2);
rho_p(k) <= log(1 + Gamma_p(k)) / log(2);
SINR_n_th <= Gamma_n(k);
end
%% PSD constraint
% stream n
for k = 1:N_user
A_k = scaled_path_loss_linear(k) * M * exp(-1j * 2 * random_phase(k)) ...
* (Q_n/Gamma_n_old(k) - (Q_n_old/(Gamma_n_old(k)^2))*Gamma_n(k) ...
+ Q_n_old/Gamma_n_old(k) - (Q_c + Q_p_sum));
sum_E = zeros(M,M); sum_r2 = 0;
for l = 1:M
E_l = zeros(M,M); E_l(l,l) = 1;
sum_E = sum_E + lambda_n(k,l)*E_l;
sum_r2 = sum_r2 + lambda_n(k,l)*max_uncertainty_radius(k, l)^2;
end
off = A_k * bar_ResponseVector(:, k);
bottom = bar_ResponseVector(:, k)' * A_k * bar_ResponseVector(:, k) ...
- scaled_noise_power - sum_r2;
A_plus = A_k + sum_E;
top = [A_plus, off];
bottom_row = [ctranspose(off), bottom];
blk = [top; bottom_row];
blk == hermitian_semidefinite(M+1);
end
% stream c
for k = 1:N_user
A_k = scaled_path_loss_linear(k) * M * exp(-1j * 2 * random_phase(k)) ...
* (Q_c/Gamma_c_old(k) - (Q_c_old/(Gamma_c_old(k)^2))*Gamma_c(k) ...
+ Q_c_old/Gamma_c_old(k) - Q_p_sum);
sum_E = zeros(M,M); sum_r2 = 0;
for l = 1:M
E_l = zeros(M,M); E_l(l,l) = 1;
sum_E = sum_E + lambda_c(k,l)*E_l;
sum_r2 = sum_r2 + lambda_c(k,l)*max_uncertainty_radius(k, l)^2;
end
off = A_k * bar_ResponseVector(:, k);
bottom = bar_ResponseVector(:, k)' * A_k * bar_ResponseVector(:, k) ...
- scaled_noise_power - sum_r2;
A_plus = A_k + sum_E;
top = [A_plus, off];
bottom_row = [ctranspose(off), bottom];
blk = [top; bottom_row];
blk == hermitian_semidefinite(M+1);
end
% stream p
for k = 1:N_user
A_k = scaled_path_loss_linear(k) * M * exp(-1j * 2 * random_phase(k)) ...
* (Q_p(:,:,k)/Gamma_p_old(k) - (Q_p_old(:,:,k)/(Gamma_p_old(k)^2))*Gamma_p(k) ...
+ Q_p_old(:,:,k)/Gamma_p_old(k) - (Q_p_sum - Q_p(:,:,k)));
sum_E = zeros(M,M); sum_r2 = 0;
for l = 1:M
E_l = zeros(M,M); E_l(l,l) = 1;
sum_E = sum_E + lambda_p(k,l)*E_l;
sum_r2 = sum_r2 + lambda_p(k,l)*max_uncertainty_radius(k, l)^2;
end
off = A_k * bar_ResponseVector(:, k);
bottom = bar_ResponseVector(:, k)' * A_k * bar_ResponseVector(:, k) ...
- scaled_noise_power - sum_r2;
A_plus = A_k + sum_E;
top = [A_plus, off];
bottom_row = [ctranspose(off), bottom];
blk = [top; bottom_row];
blk == hermitian_semidefinite(M+1);
end
cvx_end
When all three PSD constraint blocks are included, I get the following error:
Index exceeds array bounds.
Error in cvxprob/eliminate (line 179)
P = P( :, colX ) + P( :, cols ) * temp;
Error in cvxprob/solve (line 18)
[ At, cones, sgn, Q, P, dualized ] = eliminate( prob, true, shim.dualize );
Error in cvx_end (line 88)
solve( prob );
Error in beamforming_robust (line 207)
cvx_end
Error in MAIN_for_robust (line 206)
beamforming_robust(ChannelGain, user_demand,
transmit_power_W, ...
When I comment out any two of the three PSD for k loops (only one active), it solves successfully and command window shows:
Calling Mosek 9.1.9: 10316 variables, 4132 equality constraints
For improved efficiency, Mosek is solving the dual problem.
------------------------------------------------------------
MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (1) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (2) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (3) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (4) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (5) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (6) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (7) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (8) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (9) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (10) of matrix 'A'.
Warning number 710 is disabled.
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 4132
Cones : 15
Scalar variables : 3685
Matrix variables : 25
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 16
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.03
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 4132
Cones : 15
Scalar variables : 3685
Matrix variables : 25
Integer variables : 0
Optimizer - threads : 16
Optimizer - solved problem : the primal
Optimizer - Constraints : 4116
Optimizer - Cones : 16
Optimizer - Scalar variables : 3662 conic : 3530
Optimizer - Semi-definite variables: 25 scalarized : 13669
Factor - setup time : 0.20 dense det. time : 0.00
Factor - ML order time : 0.05 GP order time : 0.00
Factor - nonzeros before factor : 5.14e+06 after factor : 5.15e+06
Factor - dense dim. : 7 flops : 9.03e+09
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 4.0e+00 4.0e+02 1.9e+01 0.00e+00 1.823841826e+01 0.000000000e+00 1.0e+00 0.28
1 4.5e-01 4.5e+01 2.2e+00 -2.66e-01 2.644278447e+00 -1.230315611e+00 1.1e-01 0.73
2 1.4e-01 1.4e+01 4.1e-01 6.98e-01 2.112961329e-01 -1.105779588e+00 3.5e-02 1.08
3 6.7e-02 6.6e+00 1.3e-01 9.01e-01 -3.009242869e-01 -9.651191290e-01 1.7e-02 1.42
4 3.8e-02 3.8e+00 5.5e-02 9.53e-01 -6.184502691e-01 -1.004291433e+00 9.4e-03 1.77
5 1.5e-02 1.5e+00 1.3e-02 9.73e-01 -8.704470768e-01 -1.027755626e+00 3.8e-03 2.11
6 4.9e-03 4.8e-01 2.3e-03 9.89e-01 -9.842240369e-01 -1.036296406e+00 1.2e-03 2.44
7 6.0e-04 6.0e-02 1.0e-04 9.97e-01 -1.032563370e+00 -1.039000152e+00 1.5e-04 2.78
8 1.9e-04 1.9e-02 1.8e-05 1.00e+00 -1.036952556e+00 -1.039032164e+00 4.8e-05 3.13
9 1.5e-04 1.5e-02 9.8e-06 1.16e+00 -1.028311272e+00 -1.029700425e+00 3.7e-05 3.47
10 6.2e-05 6.1e-03 2.2e-06 3.73e+00 -1.000245880e+00 -1.000409919e+00 1.5e-05 3.81
11 1.0e-06 1.0e-04 3.9e-09 1.05e+00 -9.999391682e-01 -9.999422146e-01 2.6e-07 4.16
12 7.7e-08 2.2e-06 1.3e-11 9.94e-01 -9.999010855e-01 -9.999011464e-01 5.5e-09 4.48
13 2.2e-08 6.3e-07 4.5e-12 3.11e-01 -9.999016880e-01 -9.999017146e-01 1.6e-09 4.86
Optimizer terminated. Time: 4.91
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: -9.9990168802e-01 nrm: 8e+01 Viol. con: 4e-08 var: 2e-09 barvar: 0e+00 cones: 0e+00
Dual. obj: -9.9990171464e-01 nrm: 4e+02 Viol. con: 0e+00 var: 1e-06 barvar: 3e-09 cones: 0e+00
Optimizer summary
Optimizer - time: 4.91
Interior-point - iterations : 13 time: 4.86
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.71464e-06
I have read similar topics in the forum but still don’t know why eliminate only fails when multiple identical PSD blocks are present? And how to modify the code to ensure the optimization problem solves successfully under PSD constraints for all three types of streams?
Thank you very much for any insights!