Do I need to apply Monte-carlo simulations for solving optimization problems using CVX?
For example, In the below code I have applied monte-carlo only for channel generations, but not applying monte-carlo to whole code which uses CVX inside. Is it fine? Also, this code is taking very high run-time (even with latest Mosek version), how to reduce this run-time? Thank you
%Paper: Downlink MIMO-RSMA With Successive Null-Space Precoding
%Authors: Aravindh Krishnamoorthy, Robert Schober
%IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 21, NO. 11, NOVEMBER 2022
%Fig. 9. Average WSR for perfect CSI knowledge as a function of PT
%clc;
clear variables; close all;
tic
K = 3; % Number of users
N = 10; % Number of antennas at the BS
M1 = 2; % Number of antennas at user-1
M2 = 4; % Number of antennas at user-2
M3 = 4; % Number of antennas at user-3
M = 2;
d1 = 250; % distance(in meters) of use-1 from the BS
d2 = 150; % distance(in meters) of use-2 from the BS
d3 = 50; % distance(in meters) of use-3 from the BS
eta_k = 1/K; % User weights
eta_c_k = 1/K; % User weights
L_1 = d1^2; % Scalar path loss for user-1
L_2 = d2^2; % Scalar path loss for user-2
L_3 = d3^2; % Scalar path loss for user-1
sigm_sq = 10^((-35-30)/10); %Noise variance in linear scale (of -35dBm)
epsiln = 10^-5; %numerical tolerance
%epsiln = 10^-2; %numerical tolerance
N_itr = 10^4; % Monte-carlo iterations
N1 = N;
N2 = N-M1;
N3 = N-(M1+M2);
P_t_range = 10:5:50; %Transmit power in dBm
R_wsr_LB = zeros(1, length(P_t_range));
for pow_i = 1:length(P_t_range)
PT_dBm = P_t_range(pow_i); % Transmit power in dBm
pt = (10^-3)*10.^(PT_dBm/10); %Transmit power in linear scale
H_1 = zeros(M1,N);
H_2 = zeros(M2,N);
H_3 = zeros(M3,N);
for itr = 1:N_itr
H_1 = H_1 + (1/sqrt(2))(randn(M1,N)+1jrandn(M1,N)) ; % MIMO channel matrix between the BS and user-1
H_2 = H_2 + (1/sqrt(2))(randn(M2,N)+1jrandn(M2,N)) ; % MIMO channel matrix between the BS and user-2
H_3 = H_3 + (1/sqrt(2))(randn(M3,N)+1jrandn(M3,N)) ; % MIMO channel matrix between the BS and user-3
end %end for itr
H_1 = H_1./N_itr;
H_2 = H_2./N_itr;
H_3 = H_3./N_itr;
Sie_1 = eye(N);
% Finding a matrix Sie whose columns are the unit-length basis vectors of
% the null space of the following augmented matrix:
aug_mat_2 = ([H_1.‘]).’;
% Use the null function to compute the null space basis
null_space_basis_2 = null(aug_mat_2, ‘r’);
% Normalize the columns to unit length
Sie_2 = null_space_basis_2 ./ vecnorm(null_space_basis_2); % ----CHECK -----
aug_mat_3 = ([H_1.‘, H_2.’]).';
% Use the null function to compute the null space basis
null_space_basis_3 = null(aug_mat_3, ‘r’);
% Normalize the columns to unit length
Sie_3 = null_space_basis_3 ./ vecnorm(null_space_basis_3);% ----CHECK -----
% Algorithm 1: Algorithm for solving (19-) via SCA
%X_1_h(:,:,1) = 0; X_2_h(:,:,1) = 0; X_3_h(:,:,1) = 0;
%X_1_h(:,:,1) = zeros(10,10); X_2_h(:,:,1) = zeros(8,8); X_3_h(:,:,1) = zeros(4,4);
X_1_h(:,:,1) = (pt/(K*N1))eye(N1); X_2_h(:,:,1) = (pt/(KN2))eye(N2); X_3_h(:,:,1) = (pt/(KN3))*eye(N3);
itr_ind_l = 1; R_wsr_tild_str(itr_ind_l) = -Inf;
while true
% Increment iteration index
itr_ind_l = itr_ind_l + 1;
Q_2_h(:,:,itr_ind_l-1) = Sie_1*X_1_h(:,:,itr_ind_l-1)*Sie_1';
Q_3_h(:,:,itr_ind_l-1) = Sie_1*X_1_h(:,:,itr_ind_l-1)*Sie_1' + Sie_2*X_2_h(:,:,itr_ind_l-1)*Sie_2';
Q_4_h(:,:,itr_ind_l-1) = Sie_1*X_1_h(:,:,itr_ind_l-1)*Sie_1' + Sie_2*X_2_h(:,:,itr_ind_l-1)*Sie_2' + Sie_3*X_3_h(:,:,itr_ind_l-1)*Sie_3';
inver_term_Rc2 = (eye(M2)+1/(L_2*sigm_sq)*H_2*Q_3_h(:,:,itr_ind_l-1)*H_2')^(-1);
inver_term_R2 = (eye(M2)+1/(L_2*sigm_sq)*H_2*Q_2_h(:,:,itr_ind_l-1)*H_2')^(-1);
inver_term_Rc3_a = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_4_h(:,:,itr_ind_l-1)*H_3')^(-1);
inver_term_Rc3_b = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_4_h(:,:,itr_ind_l-1)*H_3')^(-1);
inver_term_R3_a = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3_h(:,:,itr_ind_l-1)*H_3')^(-1);
inver_term_R3_b = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3_h(:,:,itr_ind_l-1)*H_3')^(-1);
cvx_expert true
%cvx_save_prefs
cvx_clear
cvx_begin quiet
% Define variables
%cvx_precision high
% Variables
variable Q_c(10, 10) complex
variable X_1(10, 10) complex
variable X_2(8, 8) complex
variable X_3(4, 4) complex
%expression constraints(1,5);
%variables Q_c X_1 X_2 X_3
Q_1 = Sie_1*X_1*Sie_1'; Q_2 = Sie_2*X_2*Sie_2'; Q_3 = Sie_3*X_3*Sie_3';
R_c_1_tild = log_det(eye(M1)+1/(L_1*sigm_sq)*H_1*Q_c*H_1'+1/(L_1*sigm_sq)*H_1*(Q_1)*H_1')/log(2) - real(log2(det( eye(M1)+1/(L_1*sigm_sq)*H_1*Q_2_h(:,:,itr_ind_l-1)*H_1'))) ;
R_1_tild = log_det(eye(M1)+1/(L_1*sigm_sq)*H_1*Q_1*H_1'+1/(L_1*sigm_sq)*H_1*(0)*H_1')/log(2) - log2(det( eye(M1)+1/(L_1*sigm_sq)*H_1*(0)*H_1')) ;
R_c_2_tild = log_det(eye(M2)+1/(L_2*sigm_sq)*H_2*Q_c*H_2'+1/(L_2*sigm_sq)*H_2*(Q_1+Q_2)*H_2')/log(2) - real(log2(det( eye(M2)+1/(L_2*sigm_sq)*H_2*Q_3_h(:,:,itr_ind_l-1)*H_2'))) - 1/((L_2*sigm_sq)*log(2))*real(trace(Sie_1'*H_2'*inver_term_Rc2*H_2*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) ;
R_2_tild = log_det(eye(M2)+1/(L_2*sigm_sq)*H_2*Q_2*H_2'+1/(L_2*sigm_sq)*H_2*(Q_1)*H_2')/log(2) - real(log2(det( eye(M2)+1/(L_2*sigm_sq)*H_2*Q_2_h(:,:,itr_ind_l-1)*H_2'))) - 1/((L_2*sigm_sq)*log(2))*real(trace(Sie_1'*H_2'*inver_term_R2*H_2*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) ;
R_c_3_tild = log_det(eye(M3)+1/(L_3*sigm_sq)*H_3*Q_c*H_3'+1/(L_3*sigm_sq)*H_3*(Q_1+Q_2+Q_3)*H_3')/log(2) - real(log2(det( eye(M3)+1/(L_3*sigm_sq)*H_3*Q_4_h(:,:,itr_ind_l-1)*H_3'))) - 1/((L_3*sigm_sq)*log(2))*( real(trace(Sie_1'*H_3'*inver_term_Rc3_a*H_3*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) + real(trace(Sie_2'*H_3'*inver_term_Rc3_b*H_3*Sie_2*(X_2-X_2_h(:,:,itr_ind_l-1)) )) );
R_3_tild = log_det(eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3*H_3'+1/(L_3*sigm_sq)*H_3*(Q_1+Q_2)*H_3')/log(2) - real(log2(det( eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3_h(:,:,itr_ind_l-1)*H_3'))) - 1/((L_3*sigm_sq)*log(2))*( real(trace(Sie_1'*H_3'*inver_term_R3_a*H_3*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) + real(trace(Sie_2'*H_3'*inver_term_R3_b*H_3*Sie_2*(X_2-X_2_h(:,:,itr_ind_l-1)) )) );
% Your inner optimization problem objective
inner_objective = ( (eta_k*eta_c_k*min([R_c_1_tild,R_c_2_tild,R_c_3_tild])*3) + eta_k*R_1_tild+eta_k*R_2_tild+eta_k*R_3_tild) ;
%objective function
maximize inner_objective
subject to
real(trace(Q_c)+trace(X_1)+trace(X_2)+trace(X_3)) <= pt;
real(Q_c) >= 0;
real(X_1) >= 0;
real(X_2) >= 0;
real(X_3) >= 0;
cvx_end
% Extract the optimal values
R_wsr_tild_str(itr_ind_l)=inner_objective;
Q_c_str = Q_c; % optimal solution
X_1_str = X_1; % optimal solution
X_2_str = X_2; % optimal solution
X_3_str = X_3; % optimal solution
X_1_h(:,:,itr_ind_l) = X_1_str;
X_2_h(:,:,itr_ind_l) = X_2_str;
X_3_h(:,:,itr_ind_l) = X_3_str;
% Check convergence
if itr_ind_l > 2 && abs(R_wsr_tild_str(itr_ind_l) - R_wsr_tild_str(itr_ind_l-1)) < epsiln
break; % Convergence achieved
end
end
R_wsr_tild_str = R_wsr_tild_str(itr_ind_l);
Q_c_str = Q_c_str;
X_1_str = X_1_str;
X_2_str = X_2_str;
X_3_str = X_3_str;
[Uc,~] = eigs(Q_c_str, 2, ‘lm’);
[U1,~] = eigs(X_1_str, 2, ‘lm’);
[U2,~] = eigs(X_2_str, 4, ‘lm’);
[U3,~] = eigs(X_3_str, 4, ‘lm’);
% [Uc,~,~] = svd(Q_c_str);
% [U1,~,~] = svd(X_1_str);
% [U2,~,~] = svd(X_2_str);
% [U3,~,~] = svd(X_3_str);
% Algorithm 1: Algorithm for solving lower-bound via SCA
%X_1_h(:,:,1) = 0; X_2_h(:,:,1) = 0; X_3_h(:,:,1) = 0;
%X_1_h(:,:,1) = zeros(10,10); X_2_h(:,:,1) = zeros(8,8); X_3_h(:,:,1) = zeros(4,4);
X_1_h(:,:,1) = (pt/(K*N1))eye(N1); X_2_h(:,:,1) = (pt/(KN2))eye(N2); X_3_h(:,:,1) = (pt/(KN3))*eye(N3);
itr_ind_l = 1; R_wsr_lb_str(itr_ind_l) = -Inf;
while true
% Increment iteration index
itr_ind_l = itr_ind_l + 1;
Q_2_h(:,:,itr_ind_l-1) = Sie_1*X_1_h(:,:,itr_ind_l-1)*Sie_1';
Q_3_h(:,:,itr_ind_l-1) = Sie_1*X_1_h(:,:,itr_ind_l-1)*Sie_1' + Sie_2*X_2_h(:,:,itr_ind_l-1)*Sie_2';
Q_4_h(:,:,itr_ind_l-1) = Sie_1*X_1_h(:,:,itr_ind_l-1)*Sie_1' + Sie_2*X_2_h(:,:,itr_ind_l-1)*Sie_2' + Sie_3*X_3_h(:,:,itr_ind_l-1)*Sie_3';
inver_term_Rc2 = (eye(M2)+1/(L_2*sigm_sq)*H_2*Q_3_h(:,:,itr_ind_l-1)*H_2')^(-1);
inver_term_R2 = (eye(M2)+1/(L_2*sigm_sq)*H_2*Q_2_h(:,:,itr_ind_l-1)*H_2')^(-1);
inver_term_Rc3_a = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_4_h(:,:,itr_ind_l-1)*H_3')^(-1);
inver_term_Rc3_b = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_4_h(:,:,itr_ind_l-1)*H_3')^(-1);
inver_term_R3_a = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3_h(:,:,itr_ind_l-1)*H_3')^(-1);
inver_term_R3_b = (eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3_h(:,:,itr_ind_l-1)*H_3')^(-1);
cvx_expert true
%cvx_save_prefs
cvx_clear
cvx_begin quiet
% Define variables
%cvx_precision high
% Variables
variable Xc_tild(2, 2) complex
variable X1_tild(2, 2) complex
variable X2_tild(4, 4) complex
variable X3_tild(4, 4) complex
%variables Xc_tild X1_tild X2_tild X3_tild
Q_c = (Uc*Xc_tild*Uc');
X_1 = U1*X1_tild*U1';
X_2 = U2*X2_tild*U2';
X_3 = U3*X3_tild*U3';
Q_1 = Sie_1*X_1*Sie_1'; Q_2 = Sie_2*X_2*Sie_2'; Q_3 = Sie_3*X_3*Sie_3';
R_c_1_tild = log_det(eye(M1)+1/(L_1*sigm_sq)*H_1*Q_c*H_1'+1/(L_1*sigm_sq)*H_1*(Q_1)*H_1')/log(2) - real(log2(det( eye(M1)+1/(L_1*sigm_sq)*H_1*Q_2_h(:,:,itr_ind_l-1)*H_1'))) ;
R_1_tild = log_det(eye(M1)+1/(L_1*sigm_sq)*H_1*Q_1*H_1'+1/(L_1*sigm_sq)*H_1*(0)*H_1')/log(2) - log2(det( eye(M1)+1/(L_1*sigm_sq)*H_1*(0)*H_1')) ;
R_c_2_tild = log_det(eye(M2)+1/(L_2*sigm_sq)*H_2*Q_c*H_2'+1/(L_2*sigm_sq)*H_2*(Q_1+Q_2)*H_2')/log(2) - real(log2(det( eye(M2)+1/(L_2*sigm_sq)*H_2*Q_3_h(:,:,itr_ind_l-1)*H_2'))) - 1/((L_2*sigm_sq)*log(2))*real(trace(Sie_1'*H_2'*inver_term_Rc2*H_2*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) ;
R_2_tild = log_det(eye(M2)+1/(L_2*sigm_sq)*H_2*Q_2*H_2'+1/(L_2*sigm_sq)*H_2*(Q_1)*H_2')/log(2) - real(log2(det( eye(M2)+1/(L_2*sigm_sq)*H_2*Q_2_h(:,:,itr_ind_l-1)*H_2'))) - 1/((L_2*sigm_sq)*log(2))*real(trace(Sie_1'*H_2'*inver_term_R2*H_2*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) ;
R_c_3_tild = log_det(eye(M3)+1/(L_3*sigm_sq)*H_3*Q_c*H_3'+1/(L_3*sigm_sq)*H_3*(Q_1+Q_2+Q_3)*H_3')/log(2) - real(log2(det( eye(M3)+1/(L_3*sigm_sq)*H_3*Q_4_h(:,:,itr_ind_l-1)*H_3'))) - 1/((L_3*sigm_sq)*log(2))*( real(trace(Sie_1'*H_3'*inver_term_Rc3_a*H_3*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) + real(trace(Sie_2'*H_3'*inver_term_Rc3_b*H_3*Sie_2*(X_2-X_2_h(:,:,itr_ind_l-1)) )) );
R_3_tild = log_det(eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3*H_3'+1/(L_3*sigm_sq)*H_3*(Q_1+Q_2)*H_3')/log(2) - real(log2(det( eye(M3)+1/(L_3*sigm_sq)*H_3*Q_3_h(:,:,itr_ind_l-1)*H_3'))) - 1/((L_3*sigm_sq)*log(2))*( real(trace(Sie_1'*H_3'*inver_term_R3_a*H_3*Sie_1*(X_1-X_1_h(:,:,itr_ind_l-1)) )) + real(trace(Sie_2'*H_3'*inver_term_R3_b*H_3*Sie_2*(X_2-X_2_h(:,:,itr_ind_l-1)) )) );
% Your inner optimization problem objective
inner_objective = ( (eta_k*eta_c_k*min([R_c_1_tild,R_c_2_tild,R_c_3_tild])*3) + eta_k*R_1_tild+eta_k*R_2_tild+eta_k*R_3_tild) ;
%objective function
maximize inner_objective
subject to
real(trace(Xc_tild)+trace(X1_tild)+trace(X2_tild)+trace(X3_tild)) <= pt;
real(Xc_tild) >= 0;
real(X1_tild) >= 0;
real(X2_tild) >= 0;
real(X3_tild) >= 0;
cvx_end
% Extract the optimal values
R_wsr_lb_str(itr_ind_l)=inner_objective;
Xc_tild_str = Xc_tild; % optimal solution
X1_tild_str = X1_tild; % optimal solution
X2_tild_str = X2_tild; % optimal solution
X3_tild_str = X3_tild; % optimal solution
X_1_h(:,:,itr_ind_l) = U1*X1_tild_str*U1';
X_2_h(:,:,itr_ind_l) = U2*X2_tild_str*U2';
X_3_h(:,:,itr_ind_l) = U3*X3_tild_str*U3';
% Check convergence
if itr_ind_l > 2 && abs(R_wsr_lb_str(itr_ind_l) - R_wsr_lb_str(itr_ind_l-1)) < epsiln
break; % Convergence achieved
end
end
R_wsr_lb_str = R_wsr_lb_str(itr_ind_l);
Xc_tild_str = Xc_tild_str;
X1_tild_str = X1_tild_str;
X2_tild_str = X2_tild_str;
X3_tild_str = X3_tild_str;
R_wsr_LB(pow_i) = R_wsr_lb_str; % WSR-LB
end %end for pow_i
figure; hold on; % For Fig-5b
plot(P_t_range, R_wsr_LB, ‘m’,‘Marker’,‘diamond’);
xlabel(‘Pt (in dBm)’); ylabel(‘Average Sum Rate (BPCU)’);
legend(‘Proposed LB’);
title(‘Average WSR for M1=2, M2=M3=4, N=10, d1=250,d2=250,d3=50m’);
grid on
Elapsed_time = seconds(toc);
Elapsed_time.Format = ‘hh:mm:ss’
%}
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------