My solution is somehow unstabile

#1

Hello all,
I’m solving a feasibility problem using SDP. My problem is convex and cvx is solving it, but there is some instability in the solution.
Mainly cvx gives a different optimal solution when I scale my constraints. I mean by scale the following:


% This is a constraint:
(10^11 * term1{i,j} ) >= (10^11 * sigma_n_square_b_mwatt);

I’ve tried using cvx_precision best or cvx_precision high and I still face the same issue. The scaling factor here is 10^11, and when I change it it gives different solution. I can’t understand why this is happening. My problem has many constraints that I’m also using this scaling to solve it because without scaling it gives infeasible.
Note that sigma_n_square_b_mwatt = 2.2607e-11
Please, can anybody tell me why this is happening, and how to solve it.
Best Regards,

(Mark L. Stone) #2

If the problem is poorly scaled, it might be beyond the capability of the solver to produce a reliable answer. If all calculations were exact (in exact arithmetic), multiplying all terms in a constraint by some positive number would not affect the feasibility or not of a given problem. But he reality is that the solvers called by CVX use double precision, and are affected by such multiplication.

You should try to get all coefficients (when multiplied together in a given term) to be either exactly zero or within a small number of orders of magnitude of one. Similarly, for (optimal) variable values. You can start by choosing units wisely to accomplish this.

If you provide a reproducible example, complete with all input data, perhaps more specific advice can be provided.

#3

Thank you for your reply.
I have figured out how to get the optimal solution by using different approach, but still I would like to know why this scaling thing produces different values, and what is the correct approach to follow in the future.
It will be very lengthy to write all the terms here, so I have saved the needed inputs as “Inputs.mat” and they are found at:

And here is my CVX testbed (you can try it for different scaling values, e.g., scalingTerm1 = 10^-11;scalingTerm2 = 10^-4; then at scalingTerm1 = 10^-11;scalingTerm2 = 1;). Note that all my tried values didn’t give the optimal solution.

%% Matlab code
clear all

load(‘Inputs.mat’);

% scaling terms
scalingTerm1 = 10^11; % 1
scalingTerm2 = 10^-4;% 1

t_min = 10^(-6);
t_max = 100;%1000;

convg_interval = 0.0001;%10^(-5);

counter = 0;

while( (t_max - t_min) > convg_interval)

if(counter > 1000)
    disp('No solution found after the counter limit set.')
    break;
end

t = (t_min + t_max) ./ 2;

%>>clear X_q_optimized

temp = cell(Q,N);
temp2 = cell(Q,N);

term1 = cell(Q,N);

% initialize
for i = 1:Q
    for j = 1 : N
        
        temp{i,j} = zeros(RperG,RperG);
        temp2{i,j} = zeros(RperG,RperG);
        
    end
end

clear X_q_mycvx

cvx_begin sdp quiet %quiet
cvx_precision best
variable X_q_mycvx(RperG, RperG, Q, Group_N) Hermitian
%variable t

% Feasibility test problem
minimize(0)
subject to
for i = 1 : Q
    
    for j = 1 : N
        
        serving_group = min( ceil(j/RperG) , Group_N);
        
        for myIndex = 1 : Q
            for g = 1 : Group_N
                if(i ~= myIndex || g ~= serving_group)
                    
                    if(myIndex ~= i)
                        temp{i,j} = temp{i,j} + A_iqn{(i-1)*N + j,myIndex,g} * X_q_mycvx(:,:,myIndex,g);
                    else % if g == serving_group -> it won't reach here becasue of the if condition: if(i ~= myIndex || g ~= serving_group)
                        temp2{i,j} = temp2{i,j} + A_iqn{(i-1)*N + j,myIndex,g} * X_q_mycvx(:,:,myIndex,g);
                    end
                end
                
            end
            
        end
        
        term1{i,j} = real( trace( t * A_qqn_Gr{(i-1)*N + j} * X_q_mycvx(:,:,i,serving_group) - temp{i,j} - temp2{i,j} ) );
        
        (scalingTerm1 * term1{i,j} ) >= (scalingTerm1 * s);
        
    end
    
    for g = 1 : Group_N
        (scalingTerm2 * real( trace( B_q_Gr{(i-1)*Group_N + g} * X_q_mycvx(:,:,i,g) ) )) <=  (scalingTerm2 * P);
        
        
    end
end
X_q_mycvx >= 0; % +ve semidefinite
cvx_end


if (cvx_status == "Solved" || cvx_status == "Inaccurate/Solved") % isnan(cvx_optval) % || isinf(cvx_optval)
    t_max = t; %i.e. we do not need the upper half
else
    %break;
    t_min = t;
end

counter = counter + 1;

end

fprintf(‘Value of t = %d\n’, t)

fprintf(‘Trace: [’)
for i = 1 : Q

fprintf("%e,", real( trace( B_q_Gr{(i-1)*Group_N + g} * X_q_mycvx(:,:,i,g) ) ))

end
fprintf(’]’)

(Mark L. Stone) #4

Rea this http://www.krwalters.com/numerical-issues/ . Also, read sections 7.3 and 7.4 of the MOSEK Modeling Cookbook https://docs.mosek.com/modeling-cookbook/practical.html#scaling