A warning in solving one SDP

1

Hello, everyone.
When the number of transmit antennas “N =4”, the program runs normally, but when N >=5, e.g., “N = 5” or “N = 6” , sometimes matlab will display the warning:

[
Warning: This linear matrix inequality appears to be unsymmetric. This is
very likely an error that will produce unexpected results. Please check
the LMI; and, if necessary, re-enter the model.

-> In cvxprob.newcnstr at 230
In cvx.ge at 21
In paper_bern_cvx at 53

Warning: This linear matrix inequality appears to be unsymmetric. This is
very likely an error that will produce unexpected results. Please check
the LMI; and, if necessary, re-enter the model.

-> In cvxprob.newcnstr at 230
In cvx.ge at 21
In paper_bern_cvx at 63
]

I don’t know why?
Thanks!

the code is here:

 clear;
 num_CSI = 5; 
 N = 6; % number of transmit antennas                        
 M = 3;                        
 K = 3;                         
 gamma_b = 10;                 
 gamma_b = 10 ^ (gamma_b / 10); 
 gamma_e = -2;                   
 gamma_e = 10 ^ (gamma_e / 10); 
 rho_b = 0.1;                   
 rho_e = 0.01;                 
 alpha_b = -log(rho_b);
 alpha_e = -log(rho_e);
 sigma2 = 0.1;                  
 epsilon2 = 0.002;             
 j = 1;
 while (j <= num_CSI)
   H(:,:,j) = (randn(N,M) + sqrt(-1) * randn(N,M)) / sqrt(2); %
   G(:,:,j) = (randn(N,K) + sqrt(-1) * randn(N,K)) / sqrt(2); % 
   j = j + 1;
 end
 cvx_solver sedumi;
 for j = 1:num_CSI 
     cvx_begin sdp quiet
       variable W(N,N,M) hermitian;
       variable P_z;
       variables mu_b(M) nu_b(M);
       variables mu_e(M,K) nu_e(M,K);
       expression sum_W;
       expression U(N,N);
       expression D_b(N,N);
       expression d_b(N);
       expression c_b;
       expression D_e(N,N);
       expression d_e(N);
       expression c_e;
       sum_W = 0;
       for m = 1:M
           sum_W = sum_W + W(:,:,m);
       end
       V = null(H(:,:,j)');
       U = (P_z / (N - M)) * V * V';    
       minimize(trace(sum_W) + P_z);  % object function 
       subject to     
       for m = 1:M 
           % the constraints at the legitimate users
           D_b = epsilon2 * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W - U);
           d_b = sqrt(epsilon2) * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j);
           c_b = sigma2 - real(H(:,m,j)' * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j));        
           real(trace(D_b)) - sqrt(2 * alpha_b) * mu_b(m) - alpha_b * nu_b(m) - c_b >= 0;
           norm(vec([D_b , sqrt(2) * d_b])) <= mu_b(m);
           nu_b(m) * eye(N) + D_b >= 0;
           nu_b(m) >= 0;
           W(:,:,m) >= 0;
           % the constraints at the eavesdroppers
           for k = 1:K           
               D_e = epsilon2 * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U);
               d_e = sqrt(epsilon2) * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j);
               c_e = sigma2 - real(G(:,k,j)' * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j));
               real(trace(D_e)) + sqrt(2 * alpha_e) * mu_e(m,k) + alpha_e * nu_e(m,k) - c_e <= 0;
               norm(vec([D_e , sqrt(2) * d_e])) <= mu_e(m,k);
               nu_e(m,k) * eye(N) - D_e >= 0;
               nu_e(m,k) >= 0;
           end                
       end 
       P_z >= 0;
     cvx_end
     if strfind(cvx_status,'Solved')
        % get the power value
        b = cvx_optval      
     end
 end

Blockquote

2

Sometimes block matrices that look symmetric are actually slightly asymmetric due to roundoff errors. CVX has a tolerance that it uses for such cases, but sometimes the tolerance isn’t enough. So you need to explicitly symmetrize the expression. So for instance, in your case, I think the problem may be with your c_b term. So, after you first define it, do this:

c_b = 0.5 * ( c_b + c_b' );

Depending upon which version of CVX you’re running, you can also use sym(c_b) to do the same thing.

You’ll want to check your other expressions as well.

3

Thanks for your suggestions, I have made changes for D_b = 0.5 * ( D_b + D_b’ ); and D_e = 0.5 * ( D_e + D_e’ );

4

Hi,mcg. Sorry to bother you again.
I have written programs using CVX ans YALMIP, respectively, to solve the SDP. Runing the programs with the same parametres, I got different optimal objective results. Especially, sometimes it is not feasible in YALMIP, but it is feasible in CVX, and the optimal value if strange, i.e., too large.
Here I commented the codes that generate the matrices, because I loaded the variable use the command ‘load HG.mat’, so they have the same H and G.
I try the case where H(:,:,134) and G(:,:,134) are used, the objective value gained using CVX is very large compared with other H and G, but it is not feasible if I use YALMIIP.
Why? Thanks.
the CVX version

%clear;
 num_CSI = 200; 
 N = 4;                         
 M = 3;                        
 K = 5;                         
 gamma_b = 8;                 
 gamma_b = 10 ^ (gamma_b / 10); 
 gamma_e = -3;                   
 gamma_e = 10 ^ (gamma_e / 10); 
 rho_b = 0.1;                   
 rho_e = 0.01;                 
 alpha_b = -log(rho_b);
 alpha_e = -log(rho_e);
 sigma2 = 0.1;                  
 epsilon2 = 0.002;             
%  j = 1;
%  while (j <= num_CSI)
%    H(:,:,j) = (randn(N,M) + sqrt(-1) * randn(N,M)) / sqrt(2); %
%    G(:,:,j) = (randn(N,K) + sqrt(-1) * randn(N,K)) / sqrt(2); % 
%    j = j + 1;
%  end
 cvx_solver sedumi;%sdpt3;%
 for j = 134:134 %1:num_CSI %
     cvx_begin sdp quiet
       variable W(N,N,M) hermitian;
       variable P_z;
       variables mu_b(M) nu_b(M);
       variables mu_e(M,K) nu_e(M,K);
       expression sum_W;
       expression U(N,N);
       expression D_b(N,N);
       expression d_b(N);
       expression c_b;
       expression D_e(N,N);
       expression d_e(N);
       expression c_e;
       sum_W = 0;
       for m = 1:M
           sum_W = sum_W + W(:,:,m);
       end
       V = null(H(:,:,j)');
       U = (P_z / (N - M)) * V * V';    
       minimize(trace(sum_W) + P_z);  % object function 
       subject to     
       for m = 1:M 
           % the constraints at the legitimate users
           D_b = epsilon2 * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W - U);
           D_b = 0.5 * ( D_b + D_b' );%sym(D_b);
           d_b = sqrt(epsilon2) * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j);
           c_b = sigma2 - real(H(:,m,j)' * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j));        
           real(trace(D_b)) - sqrt(2 * alpha_b) * mu_b(m) - alpha_b * nu_b(m) - c_b >= 0;
           norm(vec([D_b , sqrt(2) * d_b])) <= mu_b(m);
           nu_b(m) * eye(N) + D_b >= 0;
           nu_b(m) >= 0;
           W(:,:,m) >= 0;
           % the constraints at the eavesdroppers
           for k = 1:K           
               D_e = epsilon2 * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U);
               D_e = 0.5 * ( D_e + D_e' );%sym(D_e);
               d_e = sqrt(epsilon2) * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j);
               c_e = sigma2 - real(G(:,k,j)' * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j));
               real(trace(D_e)) + sqrt(2 * alpha_e) * mu_e(m,k) + alpha_e * nu_e(m,k) - c_e <= 0;
               norm(vec([D_e , sqrt(2) * d_e])) <= mu_e(m,k);
               nu_e(m,k) * eye(N) - D_e >= 0;
               nu_e(m,k) >= 0;
           end                
       end 
       P_z >= 0;
     cvx_end
     if strfind(cvx_status,'Solved')
        % get the power value
        cvx_optval 
     else
         unfeasible = 1
     end
 end

the YALMIP version

%clear;
 num_CSI = 2000;                   % number of channel realization
 N = 4;                         
 M = 3;                        
 K = 5;                         
 
 gamma_b = 8;                 
 gamma_b = 10 ^ (gamma_b / 10); 
 gamma_e = -3;                   
 gamma_e = 10 ^ (gamma_e / 10); 
 rho_b = 0.1;                  
 alpha_b = -log(rho_b);
 rho_e = 0.01;                   
 alpha_e = -log(rho_e);
 sigma2 = 0.1;                 
 epsilon2 = 0.002;              
 %generate the estimated channel vectors
 %random and independent,standard complex Gaussian distribution CN(0,I)
%  j = 1;
%  while (j <= num_CSI)
%    H(:,:,j) = (randn(N,M) + sqrt(-1) * randn(N,M)) / sqrt(2); %
%    %G(:,:,j) = (randn(N,K) + sqrt(-1) * randn(N,K)) / sqrt(2); % 
%    j = j + 1;
%  end
 
 for j = 134:134 %1:num_CSI 
     mu_b = sdpvar(M,1);
     nu_b = sdpvar(M,1);
     mu_e = sdpvar(M,K);
     nu_e = sdpvar(M,K);
     P_z = sdpvar(1,1);
     W = sdpvar(N,N,M,'hermitian','complex');     
     sum_W = zeros(N,N);
     for m = 1:M
         sum_W = sum_W + W(:,:,m);
     end
     V = null(H(:,:,j)');
     U = (P_z / (N - M)) * V * V';
     C = [ ];
     for m = 1:M
         D_b = epsilon2 * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W - U);
         d_b = sqrt(epsilon2) * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j);
         c_b = sigma2 - real(H(:,m,j)' * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j));
         C = C + [real(trace(D_b)) - sqrt(2 * alpha_b) * mu_b(m) - alpha_b * nu_b(m) - c_b >= 0];
         C = C + [norm(vec([D_b , sqrt(2) * d_b])) <= mu_b(m)];
         C = C + [nu_b(m) * eye(N) + D_b >= 0];
         C = C + [nu_b(m) >= 0];
         C = C + [W(:,:,m) >= 0];
         for k = 1:K
             D_e = epsilon2 * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U);
             d_e = sqrt(epsilon2) * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j);
             c_e = sigma2 - real(G(:,k,j)' * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j));
             C = C + [real(trace(D_e)) + sqrt(2 * alpha_e) * mu_e(m,k) + alpha_e * nu_e(m,k) - c_e <= 0];
             C = C + [norm(vec([D_e , sqrt(2) * d_e])) <= mu_e(m,k)];
             C = C + [nu_e(m,k) * eye(N) - D_e >= 0];
             C = C + [nu_e(m,k) >= 0];
         end
     end
     C = C + [P_z >= 0];
     
     Obj = trace(sum_W + P_z);
     ops = sdpsettings('solver','sedumi','verbose',0);%,'debug',1sdpt3
     diagnostics = solvesdp(C,Obj,ops) ;   
     if (diagnostics.problem == 0) % if feasible
         % get the power value
         pow = trace(double(Obj))
         
     else
         unfeasible = 1
     end
 end