Solver not applicable (Sedumi error)

[ramakant singh]

Hii everyone,

            I am getting the solver not applicable(sedumi) error when I am solving a QCQP(quadratically constrained quadratic program). Please help me to understand this: 

Matrix=csvread('breast_cancer.dat');

[M, N]=size(Matrix);

 Y = Matrix(:,N) ;

 Matrix = Matrix(:,2:N-1); 

 [M, N]=size(Matrix);

 

% display(Matrix) ;

 

 for i = drange(1:M)

    

     if (Y(i) == 2)

         Y(i) = -1 ;

     

     else

         Y(i) = 1 ;

         

     end

     

 end

 

 training = M * 0.8 ;

 %it may be a floor value, dont forget it to make an integer

 

 training = floor(training) ;

 

 test = M * 0.2 ; 

 

 test = ceil(test) ;

 

 y = Y(1:training,:) ;

 

 % The same thing applies here as well.

 

% Since we have to learn the kernel matrix in transduction setting, we

% jusn need to be aware of trhis fact that starrting 80 per data are of

% training data and 20 per data is test data.

% now we have read the data .Now we will calculate kernel matrix

KM1=zeros(M,M);

sigma = 0.01 ;

for i=drange(1:M)

  for j=drange(1:M)

   KM1(i,j)=exp(-norm((Matrix(i,:)-Matrix(j,:)))/(2*sigma*sigma));

  end   

end

sigma = sigma * 10 ;

KM2 = zeros(M,M) ;

for i=drange(1:M)

  for j=drange(1:M)

   KM2(i,j)=exp(-norm((Matrix(i,:)-Matrix(j,:)))/(2*sigma*sigma));

  end   

end

sigma = sigma * 10 ;

KM3 = zeros(M,M) ;

for i=drange(1:M)

  for j=drange(1:M)

   KM3(i,j)=exp(-norm((Matrix(i,:)-Matrix(j,:)))/(2*sigma*sigma));

  end   

end

sigma = sigma * 10 ;

KM4 = zeros(M,M) ;

for i=drange(1:M)

  for j=drange(1:M)

   KM4(i,j)=exp(-norm((Matrix(i,:)-Matrix(j,:)))/(2*sigma*sigma));

  end   

end

sigma = sigma * 10 ; 

KM5 = zeros(M, M) ;

for i=drange(1:M)

  for j=drange(1:M)

   KM5(i,j)=exp(-norm((Matrix(i,:)-Matrix(j,:)))/(2*sigma*sigma));

  end   

end

Tau = sdpvar(1,1) ;

t = sdpvar(1,1) ;

mu = sdpvar(5,1) ;

KM = mu(1) * KM1 + mu(2) * KM2 + mu(3) * KM3 + mu(4) * KM4 + mu(5) * KM5 ;

alpha = sdpvar(training, 1) ;   

I = eye(M) ;

c = trace(KM + Tau * I) ;

% How to write a vector of 1's in Matlab

e = ones(training,1) ;

objective = 2*alpha'*e - c * t;  

constraints = [t >= 1/(trace(KM1)) * alpha' * diag(y) * KM1(1:training,1:training) * diag(y) * alpha ] ;

constraints =  [constraints,t >= 1/(trace(KM2)) * alpha' * diag(y) * KM2(1:training,1:training) * diag(y) * alpha ] ;

constraints =  [constraints,t >= 1/(trace(KM3)) * alpha' * diag(y) * KM3(1:training,1:training) * diag(y) * alpha ] ;

constraints =  [constraints,t >= 1/(trace(KM4)) * alpha' * diag(y) * KM4(1:training,1:training) * diag(y) * alpha ] ;

constraints =  [constraints,t >= 1/(trace(KM5)) * alpha' * diag(y) * KM5(1:training,1:training) * diag(y) * alpha ] ;

constraints =  [constraints,t >= 1/(training) * (alpha' * alpha)] ;

constraints =  [constraints,alpha' * y == 0] ;

constraints =  [constraints,alpha >= 0] ;

options = sdpsettings('solver','sedumi','verbose',1,'cachesolvers',1,'sedumi.eps',1e-3);

sol = optimize(constraints,objective,options);

if sol.problem==0

   % solution = value(mu) ;

   % solution1 = value(Tau) ;

    

else

    display('Hmmm, something went wrong') ;

    sol.info

    yalmiperror(sol.problem)

end

  

[Johan Löfberg]

Your objective is non-convex (bilinear in c and t), hence sedumi, or any other convex solver, is not applicable.

BTW, don't use the cachesolver options. With 99% probability, you have n o reason to do that here

and  2alpha'*e can just as well be written as 2*sum(alpha)

[Johan Löfberg]

and your model looks ill-posed. I don't see what prevents us from picking t = inf, and Tau suffciently large to render c positive, which means optimal objective is -inf

[ramakant singh]

Thanks Johan...

[ramakant singh]

Johan, can you do a little favor for me ??   I am attaching a paper "Learning the kernel matrix with semidefinite Programming". And the problem that I have tried to solve is on page no.50 , equation number 47 and 48. There it is mentioned that above problem is a QCQP which is an instance of SOCP and it can be solved using Sedumi. 

[Johan Löfberg]

c is supposed to be a fixed tuning parameter with the constraint

c == trace(KM + Tau * I)

[Johan Löfberg]

And in the paper, the objective is maximized, not minimized

[ramakant singh]

Thanks a lot  johan, So how to solve maximization problem as there is no options in YALMIP for it??

[Johan Löfberg]

Negate the objective

[ramakant singh]

anyway thanks Johan, you have helped me a lot. 

[ramakant singh]

Johan , When I want to find the value of dual variables by writing dual(constraints(1)),.... I am getting the answer not a number. Why it is show ?? Why i am not being able to retrieve the values of dual variables.

[Johan Löfberg]

Was the problem really solved? I.e., what does the output from optimize (or solvesdp) say (the info field)

[ramakant singh]

Yes I have made few corrections like specifying the value of c and removing the constraint c = trace(KM + Tau * I) and negating the objective. The optimal value is around 16.667. Now I want to get the value of dual variables ?? Which I am getting NaN ??

[Johan Löfberg]

When you solve this QCP using SeDuMi (or any other SOCP/SDP spolver), all quadratic constraints are reformulated to SOCPs (since that's how SeDuMi solve QCPs). An SOCP has a vector dual, and the scalar dual to the original QCP is not available. If you want to use SeDuMi, you have to manually reformulate you model to an SOCP, solve the SOCP, extract duals, and figure out how those duals correspond to you original model

Some solvers (gurobi, mosek, cplex I think) try to extract QCP duals a-posteriori from the SOCP duals. I think YALMIP only supports this through gurobi

sdpvar x

C = x^2 <= 3;

optimize(C,x,sdpsettings('solver','gurobi','gurobi.QCPdual',1));

dual(C)

ans =

    0.2887

[ramakant singh]

Johan, Now I have changed my solver to mosek instead of sedumi. Still , when I used dual(constraints(1))... I got NaN. So, How to find the value of dual variables ??

[Johan Löfberg]

You still haven't told us what problem you are solving, and what you see when you solve the problem.

And as I said above, if you have a QCP, it will be reformulated to an SOCP, and your duals will be lost with all solvers except gurobi

[ramakant singh]

Johan , I have told you about the problem that I am solving. I have tried to use gurobi now, and I got a trial license for that> now when I am solving the problem using it. I am getting the error:

Warning: adding variables with small (< 1e-13) coefficients, ignored

Hmmm, something went wrong

ans =

Unknown problem in solver (try using 'debug'-flag in sdpsettings) (Error using gurobi

Gurobi error 10010: Model too large for current Gurobi license

)

[Johan Löfberg]

If this is the data you use
http://pages.cs.wisc.edu/~brecht/cs726docs/breast_cancer.dat

it is way too large to be implemented in this naive form. Things like  alpha' * diag(y) * KM1(1:training,1:training) * diag(y) * alpha will introduce an awful amount of symbolic monomials for YALMIP to handle

Instead of writing x'*Q*x, introduce a new variable z, factorize data Q as M'M, add the constraint z == M*x, and use z'*z instead of x'Qx. Alternatively, manually rewrite the quadratic constraint t >= x'*Q*x to an socp constraint norm([t-1;2*M*x])<= 1+t

[ramakant singh]

Johan, I managed to get academic license of gurobi. Now, for the first 6 constraints I get dual(constraints(i)) , i =1,2,...6 . but for constraint 7 I am getting NaN and for constraint 8, a vector of NaNs. why it is so ??

[Johan Löfberg]

What data, and how did you implement the constraints

[Johan Löfberg]

what does your complete model look like (as your model above was wrong according to earlier discussions)

[Johan Löfberg]

The recovery of duals when solving a problem where YALMIP has converted QCPs to SOCPs, and then using dual recovery in gurobi, is currently not extracting the duals for the original linear constraints. I'll add that later.

For now, you would have to manually look at the data returned by gurobi. The only tip I can give you is to look at the callgurobi file and figure out what is happening with the returned duals

[ramakant singh]

Johan, I have implemented equation 47 and 48 on page num 50 in the document that I have provided you earlier. I have provided the value of c explicitly and removed the constraints c==trace(KM + Tau *I).

[ramakant singh]

I am attaching the file as well as snap shot of output.

[Johan Löfberg]

easiest work-around is that you implement the SOCPs manually (using e.g. the cone operator), then all duals will be extracted as expected if you use QCPdual recovery

[Johan Löfberg]

On line 458 in solvesdp, change to

if length(output.qcDual) == length(ForiginalQuadratics)              
 Ftemp = F + ForiginalQuadratics;
 K = interfacedata.K;
 Ktemp = K;
 Ktemp.l = Ktemp.l + length(ForiginalQuadratics);
 tempdual = output.Dual;
 tempdual = [tempdual(1:K.f + K.l);-output.qcDual;tempdual(1+K.f+K.l:end)];
  setduals(Ftemp,tempdual,Ktemp);
end

However, as I said before, if you really have the big data I asked you about, you should implement it using a cone operator, to avoid the massive performance hit from creating the quadratic functions

[ramakant singh]

It is working now. Thanks Johan.....happy V day. And I am not very good in optimization theory that is why I could not cooperate with you(i.e the instructions that you have told me to do.)

[ramakant singh]

Hi Johan, I am now using three different kernel matrices. 

KM1(i,j) = (1 + Matrix(i,:)*Matrix(j,:)')^d ;    

KM2(i,j) = exp(-(Matrix(i,:) - Matrix(j,:))*(Matrix(i,:) - Matrix(j,:))'/2*sigma) ;

KM3(i,j) = KM3(i,j)/sqrt(KM3(i,i)*KM3(j,j)) ;

Now , in the paper it is written that we should normalize the matrices. so for normalizing the matirx, they have suggested :

KM1(i,j) = KM1(i,j)/sqrt(KM1(i,i)*KM1(j,j)) ; 

Now, after normalizing the matrices...... I am getting the error : Solver not applicable(gurobi). If I don't normalize the  matrices, the code runs fine, but after normalization it gives the  error.

[ramakant singh]

Sorry, the third matrix KM3 is :

 KM3(i,j) =  Matrix(i,:) * Matrix(j,:)' ;

[Johan Löfberg]

I would assume you run into numerical issues (or you simply did something wrong) and KM is no longer symmetric and/or psd after this manipulation. Hence, you have to symmetrize it, or even reformulate your problem from t >= x'*K*x to t >= z'*z where z==Mx and M is a low-rank factorization of K, K = M'*M

[ramakant singh]

Yes, the psd property is lost. So how to handle it ?? Because they have mentioned it in their experimental section as well.

[Johan Löfberg]

As I said, if you've lost psd (which only should be marginally as the matix theoretically should be psd), use SVD or similar to compute a natural low-rank decomposition K = M'*M etc. Or simply add a diagonal perturbation to make it psd