SDP not working, infeasible problem (mosek)

[S Hellan]

Hi!

I'm trying to implement an SDP problem using mosek, but am getting the error message 'Infeasible problem (MOSEK)'. This despite the fact that I am feeding it with a feasible solution, to the precision I have attempted to specify.
I am now very confused as to why it is not working. Any help would be appreciated.

    clear;
   
    Amin = [0, -3.5];
    Amax = [6, -3.5];
    Bmin = [-4, 3.5];
    Bmax = [4, 3.5];
    N = length(Amin);

    Asums = sdpvar(N, 1);
    Bsums = sdpvar(N, 1);
    A = sdpvar(N, N, 'full');
    B = sdpvar(N, N, 'full');
    C = sdpvar(N, N, 'full', 'complex');
    D = sdpvar(N, N, 'full', 'complex');
    Constraints = [D >= 0];
    Constraints = [Constraints, ...
        imag(D(1,1)) == 0, imag(D(2,2)) == 0];
    Constraints = [Constraints, ...
        real(C) == A, imag(C) == B];
    e1 = 0.961538461538461;
    e2 = 4.80769230769231i;
    E = [e1 - e2, -e1 + e2;
            -e1 + e2, e1 - e2];
   
    Constraints = [Constraints, Asums == A * ones(N,1), ...
        Bsums == B * ones(N,1), ...
        Amin' <= Asums <= Amax', ...
        Bmin' <= Bsums <= Bmax'];
   
    Constraints = [Constraints, C(1,1) == 0, C(2,2) == 0, ...
        C(1,2) == (D(1,1) - D(1,2)) * conj(-E(1,2)), ...
        C(2,1) == (D(2,2) - D(2,1)) * conj(-E(2,1))];
   
    % Solution
    assign(Asums,[4.388889; -3.5]);
    assign(Bsums,[0.944444; 3.5]);
    assign(C,[0, 4.388889+0.944444j; -3.5+3.5j, 0]);
    assign(A, real(value(C)));
    assign(B, imag(value(C)))
    assign(D, [0.9523^2, 0.5425+0.84j; 0.5425-0.84j, 1.05^2]);
   
    objective = Asums;
    options = sdpsettings('solver', 'mosek', ...
        'mosek.MSK_DPAR_ANA_SOL_INFEAS_TOL', 0.1);
    optimize(Constraints, objective, options)
   
    % Extract the solution
    Sol.A = value(A);
    Sol.B = value(B);
    Sol.D = value(D);
   
    check(Constraints)

[Johan Löfberg]

Your claimed feasible solution is far from feasible

    Constraints = [Constraints, Asums == A * ones(N,1), ...

        Bsums == B * ones(N,1), ...

        Amin' <= Asums <= Amax', ...

        Bmin' <= Bsums <= Bmax'];

    

    Constraints = [Constraints, C(1,1) == 0, C(2,2) == 0, ...

        C(1,2) == (D(1,1) - D(1,2)) * conj(-E(1,2)), ...

        C(2,1) == (D(2,2) - D(2,1)) * conj(-E(2,1))];

    

    % Solution

    assign(Asums,[4.388889; -3.5]);

    assign(Bsums,[0.944444; 3.5]);

    assign(C,[0, 4.388889+0.944444j; -3.5+3.5j, 0]);

    assign(A, real(value(C)));

    assign(B, imag(value(C)))

    assign(D, [0.9523^2, 0.5425+0.84j; 0.5425-0.84j, 1.05^2]);

 check(Constraints)

 

+++++++++++++++++++++++++++++++++++++++++++++++++++++

|    ID|               Constraint|   Primal residual|

+++++++++++++++++++++++++++++++++++++++++++++++++++++

|    #1|   Elementwise inequality|           0.90688|

|    #2|      Equality constraint|                 0|

|    #3|      Equality constraint|                 0|

|    #4|      Equality constraint|                 0|

|    #5|      Equality constraint|                 0|

|    #6|      Equality constraint|                 0|

|    #7|      Equality constraint|                 0|

|    #8|   Elementwise inequality|                 0|

|    #9|   Elementwise inequality|                 0|

|   #10|      Equality constraint|                 0|

|   #11|      Equality constraint|                 0|

|   #12|      Equality constraint|       -0.00033864|

|   #13|      Equality constraint|       -1.8971e-15|

+++++++++++++++++++++++++++++++++++++++++++++++++++++

[Johan Löfberg]

Note though, you say you have an SDP, but that is not what I see. D>=0 is not an LMI if you think so. D is not Hermitian, hence it is an elementwise constraint

BTW, a bit redundant to introduce the auxiliary variables Asums and Bsums

[S Hellan]

Specifying it as Hermitian fixed the problem. Thank you!

Would the solution not be feasible within the precision specified?

[Erling D. Andersen]

MOSEK prints a solution summary in the log. If you include that in a reply, I can tell you how likely it is that your problem is infeasible.

See

https://docs.mosek.com/8.1/toolbox/tutorial-solution-analysis.html

[Johan Löfberg]

Not sure what you mean with "Would the solution not be feasible within the precision specified? " Solutions are almost never strictly feasible, as solvers work with infeasible methods.

https://yalmip.github.io/faq/solutionviolated/

[S Hellan]

I think I might have interpreted your response here wrongly
https://groups.google.com/forum/#!topic/yalmip/osvism-XHR4
I thought you meant to use the listed parameter to set it. But looking at the Mosek documentation that's probably not what I was doing?
Would the correct way to change the precision in Yalmip be to add constants either side on the bound, as in
my_error_const = 0.001
0.5-my_error_cons <= X <= 0.5+my_error_cons
instead of X == 0.5?

[S Hellan]

That link really helped, thank you! I'm still quite new to all of this.

The problem was how I was specifying the matrix. After changing that the solver was happy, and I got this report

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -8.7777777134e+02   nrm: 3e+02    Viol.  con: 1e-08    var: 0e+00    barvar: 0e+00 
  Dual.    obj: -8.7777777116e+02   nrm: 5e+00    Viol.  con: 0e+00    var: 2e-08    barvar: 8e-09 

[Johan Löfberg]

YALMIP does not have any notion of precision or anything. If you want to relax your constraints, you simply relax them as you've done (keeping in mind that the magic number you pick, 0.001, will be completely redundant if you pick it smaller than the general tolerance that the solver work with, typically around 1e-8)

[Erling D. Andersen]

That report shows you got a pretty accurate solution. For most practical purposes it should be good enough.