error using Yalmip

[Dan Iancu]

Hi,
I am encountering a very weird error while using Yalmip. I have a very simple convex QP that I am trying to solve. It is trivially feasible, and I am quite convinced that the optimal value is finite. I have tried four different solvers under Yalmip, with different outcomes: CPLEX reports an unbounded problem, MOSEK reports finding an optimal value (but actually returns something infeasible), and SDPT3 and SeDuMi both run into numerical issues. For comparison, I have also implemented the exact same model in AMPL. Both CPLEX and MOSEK report the same optimal solution, which is different from the one found by MOSEK under Matlab / Yalmip.
I attached an archive containing code for both methods. The results are in 'conclusions.txt'. I am wondering if anyone can replicate the results under Matlab / Yalmip. Any thoughts on the matter would be highly appreciated!
Thank you,
  Dan

[Johan Löfberg]

Confirmed here. Thank you for reporting this. It is due to an unfortunate combination of numerical issues when writing the quadratic functions as an SOCP.

x'*Q*x +c'*x + d

only has n-1 positive eigenvalue. One is 0, but is actually numerically negative

eig(Q)

ans =

    -9.112591989459025e-18
     5.624471244543198e-03
     8.071685114660276e-03
     1.122453694574583e-02
     1.883316296070654e-02
     2.023763538350374e-02
     2.577445794487629e-02
     3.466901169060292e-02
     3.881238805692980e-02
     4.672138321763209e-02
     6.243829119109319e-02

Despite this, chol returns a factorization Q=R'*R (I detect low-rank quadratics through failure of chol, which thus doesn't work here). This factorization is of course ill-conditioned, and then it goes downhill. At one point, I try to detect a quadratic form ||Rx-y||^2 in order to generate a sliiightly smaller SOCP, requiring an R\y computation, and this computation yields values in the order of 10^15, creating a lousy model.

I'll have to find a way to robustify this. You can deactivate the if-case on line 78 in convertquadratics (i.e., always use general case on line 82) , or model manually as

SC = [g>=0, 0 <= x <= 1,norm([2*sqrt(gamma)*aux;1-alpha*g]) <= 1+alpha*g,r*x'*s <= g];

or

SC = [g>=0, 0 <= x <= 1,cone([2*sqrt(gamma)*aux;1-alpha*g],1+alpha*g),r*x'*s <= g];

(Note, SET is obsolete. You just do standard MATLAB concatenation)

[Dan Iancu]

Thank you for the fast reply, Johan. Please keep me posted if / when this issue is resolved (in the meanwhile, I'll resort to one of your suggestions).
  -Dan