# Solution for SDP is not positive semi-definite

Hi everyone!

I am new to CVX for Matlab and hope you can help me with an issue.

I want to solve a semi-definite program min(X * M) , where the constraints have the form Tr(A_i X) = a_i for a certain set of matrices A_i and real values a_i, for i goes from 1 to m. The n times n matrix X is required to be positive semi-definite (and Hermitian).

I decided to formulate the proble in the following way;

I used the {A_i} to form an orthnormal system {B_i} and added other matrices, such that in total we hold an orthonormal basis of the set of all complex valued n times n matrices. So, I can compose X as follows:
X = \sum_{i=1}^{m} B_i b_i + \sum_{i=m+1}^{n^2} D_i d_i .

As the first part of the sum describes those parts of X that are fixed by the constraints, and the B_i’s form an orthonormal system, it remains to optimize over the second sum, requiring that X is PSD.

For the sake of an easier notation, let us call the part that remains unchanged X_0 := \sum_{i=1}^{m} B_i b_i.

The advantage of this formulation is that I only have to optimize over vectors instead of matrices.

Unfortunately, when I do this I obtain a matrix that is not positive semi-definite, although this is the only constraint we have. In fact, the resulting matrix has only 4 positive eigenvalues and many (compared small) negative eigenvalues. Although the negative eigenvalues are quite small, I think it should be possible for the solver to find a matrix that is really positive semi-definite, and not only approximately, shouldn’t it?

That was the “context” of my question; now we come to the “minimal code” required to explain what I am doing:

%Preparation
for k = 1:length(D)
y(k) = trace(D{k} * M);
end

%SDP solving
cvx_begin sdp

    variables x(length(D))
dual variable Q

C = 0;
for k  = 1:length(D)
C = C + x(k) .* D{k};
end

C = C + X_0;

minimize(  real( y' * x) )
C >= 0 : Q;


cvx_end

%Calculating X
s=0;
for k = 1:length(D)
s = s + x(k) .* D{k};
end

X = X_0 + s;

The eigenvalues that I obtain are mostly negative, as can be seen below.
In some cases the solver gives me a perfectly positive semi-definite matrix, but in many cases it doesn’t.
Is there any way how I can avoid that the solver violates the PSD requirement, or is this just the best that I can obtain within the numerial method?

-0.000000516155921
-0.000000516017584
-0.000000515702941
-0.000000515425191
-0.000000450301896
-0.000000450288575
-0.000000450148228
-0.000000450127882
-0.000000398222341
-0.000000398070259
-0.000000398055811
-0.000000397997038
-0.000000353923196
-0.000000353906003
-0.000000353833455
-0.000000353685267
-0.000000317047886
-0.000000316949589
-0.000000316870178
-0.000000316638173
-0.000000284905313
-0.000000284840041
-0.000000284630373
-0.000000284564555
-0.000000250103285
-0.000000250003089
-0.000000249892010
-0.000000249543980
-0.000000210383219
-0.000000210181552
-0.000000209457046
-0.000000208865466
-0.000000173907113
-0.000000173803027
-0.000000173642420
-0.000000173432423
-0.000000126522490
-0.000000126056356
-0.000000125631958
-0.000000125268480
0.009709750768604
0.064618961999724
0.286773403261360
0.638910197995383

This is very likely an effect of that computations are done in finite precision floating point numbers and hence the solution is only an approximation to an optimal and feasible solution.

I would post

• The log output so we have more info. to work with.
• Or try another optimizer e.g. Mosek, Sdpt3, Sedumi. Maybe you have better luck with one of them.

Here’s the solver output:

Calling SDPT3 4.0: 2184 variables, 280 equality constraints
------------------------------------------------------------

num. of constraints = 280
dim. of sdp    var  = 88,   num. of sdp  blk  =  1
dim. of free   var  = 248 *** convert ublk to lblk
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version  predcorr  gam  expon  scale_data
HKM      1      0.000   1        0
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
0|0.000|0.000|3.0e+02|1.7e+03|1.2e+08| 3.815853e-08  0.000000e+00| 0:0:00| chol  1  1
1|0.963|0.986|1.1e+01|2.5e+01|3.5e+05|-1.833399e+01 -1.942466e+01| 0:0:00| chol  1  1
2|0.962|0.990|4.3e-01|3.9e-01|1.2e+03|-5.603797e-01 -1.788959e+01| 0:0:00| chol  1  1
3|0.739|0.713|1.1e-01|1.2e-01|8.8e+01| 1.255950e-01 -2.189920e+00| 0:0:00| chol  1  1
4|0.423|0.894|6.4e-02|1.4e-02|2.0e+01| 4.395906e-01 -1.870282e+00| 0:0:00| chol  1  1
5|0.741|0.726|1.7e-02|1.2e-02|5.7e+00| 6.118455e-01 -2.269299e-01| 0:0:00| chol  1  1
6|0.913|0.883|1.5e-03|3.9e-03|5.5e-01| 6.596966e-01  5.441641e-01| 0:0:00| chol  1  1
7|0.926|0.856|1.1e-04|8.5e-04|4.6e-02| 6.657450e-01  6.535479e-01| 0:0:00| chol  2  2
8|0.881|0.960|1.3e-05|1.6e-04|2.8e-03| 6.665677e-01  6.679756e-01| 0:0:00| chol  2  2
9|0.650|0.666|5.4e-06|1.3e-05|9.8e-04| 6.675951e-01  6.689782e-01| 0:0:00| chol  2  2
10|0.958|0.893|2.6e-06|4.1e-06|7.0e-05| 6.678742e-01  6.691224e-01| 0:0:00| chol  3  3
11|0.866|0.841|2.2e-06|7.9e-07|1.4e-05| 6.679923e-01  6.691720e-01| 0:0:00| chol  3  3
12|0.964|0.892|2.1e-06|1.2e-07|1.7e-06| 6.680443e-01  6.691969e-01| 0:0:00| chol  3  3
13|0.937|0.896|2.1e-06|1.8e-08|2.5e-07| 6.680669e-01  6.692085e-01| 0:0:00| chol  3  3
14|0.886|0.873|2.0e-06|3.6e-09|5.2e-08| 6.680790e-01  6.692142e-01| 0:0:00| chol  4  4
15|0.789|0.846|2.0e-06|1.0e-09|1.7e-08| 6.680851e-01  6.692174e-01| 0:0:00| chol  5  5
16|0.346|0.780|2.0e-06|1.0e-09|1.8e-08| 6.680891e-01  6.692201e-01| 0:0:00| chol  5  5
17|0.023|0.412|2.0e-06|1.6e-09|1.0e-07| 6.680924e-01  6.692258e-01| 0:0:00| chol 11  8
18|0.003|0.017|2.0e-06|2.8e-09|2.9e-07| 6.680988e-01  6.692223e-01| 0:0:00|
stop: steps too short consecutively
-------------------------------------------------------------------
number of iterations   = 18
primal objective value =  6.68098796e-01
dual   objective value =  6.69222296e-01
gap := trace(XZ)       = 2.88e-07
relative gap           = 1.23e-07
actual relative gap    = -4.81e-04
rel. primal infeas (scaled problem)   = 2.01e-06
rel. dual     "        "       "      = 2.80e-09
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual     "        "       "      = 0.00e+00
norm(X), norm(y), norm(Z) = 1.1e+00, 4.3e+02, 3.1e+02
norm(A), norm(b), norm(C) = 2.0e+01, 1.3e+00, 6.4e+00
Total CPU time (secs)  = 0.40
CPU time per iteration = 0.02
termination code       = -5
DIMACS: 2.4e-06  0.0e+00  9.2e-09  0.0e+00  -4.8e-04  1.2e-07
-------------------------------------------------------------------

------------------------------------------------------------
Status: Inaccurate/Solved
Optimal value (cvx_optval): +0.665274


If I take the result of this run as new starting value and solve the SDP again, I obtain Status: Solved
So without “Inaccurate”, but the eigenvalues remain negative. I hoped that this might be a way to obtain better results, as the new starting value is “closer” to be positive semi-definite than the old one was.

Calling SDPT3 4.0: 2184 variables, 280 equality constraints
------------------------------------------------------------

num. of constraints = 280
dim. of sdp    var  = 88,   num. of sdp  blk  =  1
dim. of free   var  = 248 *** convert ublk to lblk
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version  predcorr  gam  expon  scale_data
HKM      1      0.000   1        0
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
0|0.000|0.000|2.8e+02|5.0e+02|1.6e+08| 3.446456e-07  0.000000e+00| 0:0:00| chol  1  1
1|0.885|0.980|3.2e+01|1.0e+01|6.9e+05|-5.197614e+02 -2.270980e+01| 0:0:00| chol  1  1
2|0.957|0.977|1.4e+00|3.5e-01|4.8e+03|-1.195382e+01 -2.703754e+01| 0:0:00| chol  1  1
3|0.872|1.000|1.8e-01|1.2e-02|7.1e+01|-9.639782e+00 -1.827315e+01| 0:0:00| chol  1  1
4|0.708|0.785|5.2e-02|3.6e-03|1.3e+01|-1.373441e+01 -1.196621e+01| 0:0:00| chol  1  1
5|0.474|0.641|2.7e-02|3.2e-03|5.1e+00|-1.412982e+01 -1.695008e+00| 0:0:00| chol  1  1
6|0.076|0.404|2.5e-02|4.7e-03|5.9e+00|-1.216408e+01 -1.424416e+01| 0:0:00| chol  1  1
7|0.126|0.210|2.2e-02|7.9e-03|5.2e+00|-1.133930e+01 -3.617445e+00| 0:0:00| chol  1  1
8|0.207|0.496|1.7e-02|8.4e-03|6.9e+00|-7.790153e+00 -5.740730e+00| 0:0:00| chol  1  1
9|0.262|0.207|1.3e-02|1.0e-02|6.7e+00|-5.063383e+00 -3.373085e+00| 0:0:00| chol  1  1
10|0.505|0.182|6.4e-03|1.3e-02|5.6e+00|-1.719121e+00 -2.809829e+00| 0:0:00| chol  1  1
11|0.706|0.405|1.9e-03|9.1e-03|3.0e+00|-1.287328e-01 -1.587342e+00| 0:0:00| chol  1  1
12|0.608|0.484|7.3e-04|5.0e-03|1.4e+00| 8.223002e-02 -7.352969e-01| 0:0:00| chol  1  1
13|0.793|0.425|1.5e-04|3.0e-03|6.5e-01| 1.524868e-01 -3.517447e-01| 0:0:00| chol  1  1
14|0.823|0.831|2.7e-05|5.5e-04|1.1e-01| 1.636328e-01  7.947411e-02| 0:0:00| chol  1  1
15|0.569|0.205|1.2e-05|4.4e-04|8.0e-02| 1.659498e-01  9.755556e-02| 0:0:00| chol  1  1
16|0.524|0.213|5.5e-06|3.5e-04|6.1e-02| 1.668866e-01  1.124965e-01| 0:0:00| chol  1  2
17|0.792|0.292|1.1e-06|2.5e-04|4.1e-02| 1.674422e-01  1.286034e-01| 0:0:00| chol  2  2
18|0.324|0.302|7.7e-07|6.7e-05|2.8e-02| 1.675660e-01  1.404448e-01| 0:0:00| chol  2  2
19|0.789|0.946|1.6e-07|1.9e-05|1.7e-03| 1.676567e-01  1.662788e-01| 0:0:00| chol  3  3
20|0.395|0.695|9.8e-08|1.6e-06|5.9e-04| 1.676946e-01  1.673204e-01| 0:0:00| chol  4  4
21|0.502|0.387|5.0e-08|7.6e-07|3.7e-04| 1.677276e-01  1.674968e-01| 0:0:00| chol  5  5
22|0.853|0.952|1.7e-08|2.5e-07|5.8e-05| 1.677450e-01  1.677525e-01| 0:0:00| chol  5  6
23|0.998|0.975|7.6e-09|4.1e-08|1.8e-06| 1.677443e-01  1.677668e-01| 0:0:00| chol  7  7
24|0.991|0.984|7.4e-09|1.7e-09|3.2e-08| 1.677443e-01  1.677674e-01| 0:0:00| chol  7  7
25|0.978|0.980|7.4e-09|4.8e-11|9.4e-10| 1.677444e-01  1.677674e-01| 0:0:01|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations   = 25
primal objective value =  1.67744354e-01
dual   objective value =  1.67767450e-01
gap := trace(XZ)       = 9.43e-10
relative gap           = 7.06e-10
actual relative gap    = -1.73e-05
rel. primal infeas (scaled problem)   = 7.42e-09
rel. dual     "        "       "      = 4.75e-11
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual     "        "       "      = 0.00e+00
norm(X), norm(y), norm(Z) = 9.8e-01, 2.1e+03, 1.5e+03
norm(A), norm(b), norm(C) = 2.0e+01, 1.5e+00, 3.0e+01
Total CPU time (secs)  = 0.51
CPU time per iteration = 0.02
termination code       =  0
DIMACS: 9.0e-09  0.0e+00  1.6e-10  0.0e+00  -1.7e-05  7.1e-10
-------------------------------------------------------------------

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -0.00010365


I solved with Sdpt3, but the same happens if I solve with Sedumi or Mosek. In both cases, the first try is “Inaccurate/Solved” and the second try is “Solved”.

I have no prior knowledge in how SDP solvers are implemented, but 40 out of 44 eigenvalues being negative with a magnitude of 1e-7 seems to be a lot in my eyes. I could bear with 1e-12 or better, but here it seems to me as if I do something wrong or at least can improve something to obtain better results.

I’m rather confused what you’re doing. What is D? Is C >= 0 actually a semidefinite constraint applied to a symmetric matrix? Are there any CVX warnings?

If you have a symmetric matrix M which you want to constrain to be psd without any slightly negative eigenvalues, impose the constraint
M - small_positive_number*eye(size(M)) >= 0
where small_positive_number is perhaps about 1e-7.

Alternatively, if you are using Mosek, even though @Erling might not approve of the following, you could try
cvx_solver_settings('MSK_DPAR_INTPNT_CO_TOL_PFEAS',1e-12)
which if Mosek succeeds (perhaps there is less chance it will do so?), will result in the maximum magnitude of negative eigenvalue being <= 1e-12.

Probably I wasn’t clear enough in my initial post, so I try to clarify what I aim to do;
B and D together form an orthonormal basis of the space of C^{n \times n} matrices. The first m =#constraints basis elements are collected in the cell array B. I chose the basis in a way such that the first m basis vectors satisfy Tr[B{i}X] = \tilde{a}i, where \tilde{a}i is given. The remaining basis elements are collected in the cell array D. The matrix X looks as follows:
X = \sum
{i=1}^{m} B_i b_i + \sum
{i=m+1}^{n^2} D_i d_i

My approach was the following one;
I can find the coefficients b_i (by solving a system of linear equations), such that all constraints (exept the one requiring that C is PSD) are satisfied. Thus, my problem simplifies as I merely have to vary the d_i’s in the second sum and there is only one constraint, C >= 0 left (because all the other constraints are satisfied by construction).

By the way, I recognized, one may simplify my problem by changing the line minimize( real( y’ * x) ) to minimize(1), as my primal goal is to find a matrix that is PSD and satisfies my constraints.

I haven’t observed any warnings.

You can try small_pos-number = 1e-6
As for cvx_solver_settings`, my suggestion will only work if Mosek is used as the solver.