Unexpected Results in Semidefinite Constraint Optimization

1

I tried running the following code:

A = [
1 0 0 1
0 0 0 0
0 0 0 0
1 0 0 1]

B = [
1 0 0 0 1 0 0 0 1
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
1 0 0 0 1 0 0 0 1
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
1 0 0 0 1 0 0 0 1]

cvx_begin SDP
variable p
expression C(4*4*9*9, 4*4*9*9)

C = p * Tensor(A, A, B, B) % Tensor() from QETLAB

subject to
  p <= 1
  C == semidefinite(4*4*9*9)
maximize p
cvx_end

When executed, I got the following output:

Calling SDPT3 4.0: 840457 variables, 1 equality constraint
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints =  1
 dim. of sdp    var  = 1296,   num. of sdp  blk  =  1
 dim. of linear var  =  1
*******************************************************************
   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|6.4e+02|3.6e+01|3.4e+06| 1.296000e+03  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|0.989|1.0e-11|5.3e-01|4.9e+04| 7.224939e+00 -7.388049e-01| 0:0:02| chol  1  1 
 2|0.995|1.000|6.9e-14|4.5e-02|3.5e+03|-9.177506e-01 -9.785768e-01| 0:0:02| chol  1  1 
 3|1.000|1.000|9.7e-14|1.3e-02|5.2e+02|-9.765279e-01 -9.937538e-01| 0:0:02| chol  1  1 
 4|1.000|1.000|2.2e-13|4.0e-03|7.8e+01|-9.963402e-01 -9.981769e-01| 0:0:02| chol  1  1 
 5|1.000|1.000|1.4e-13|4.0e-04|6.3e+00|-9.999532e-01 -9.998184e-01| 0:0:02| chol  1  1 
 6|1.000|1.000|6.4e-14|4.0e-05|5.1e-01|-9.999959e-01 -9.999819e-01| 0:0:02| chol  1  1 
 7|1.000|1.000|9.0e-14|4.0e-06|4.2e-02|-9.999997e-01 -9.999982e-01| 0:0:02| chol  1  1 
 8|1.000|1.000|1.3e-13|4.0e-07|3.4e-03|-1.000000e+00 -9.999998e-01| 0:0:02| chol  1  1 
 9|1.000|1.000|8.0e-14|4.0e-08|2.7e-04|-1.000000e+00 -1.000000e+00| 0:0:02| chol  1  1 
10|1.000|0.980|3.0e-14|7.9e-10|5.4e-06|-1.000000e+00 -1.000000e+00| 0:0:02| chol  1  1 
11|1.000|0.980|9.3e-15|1.7e-11|1.2e-07|-1.000000e+00 -1.000000e+00| 0:0:02| chol  1  1 
12|1.000|0.982|6.1e-15|3.4e-13|2.4e-09|-1.000000e+00 -1.000000e+00| 0:0:02|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 12
 primal objective value = -1.00000000e+00
 dual   objective value = -1.00000000e+00
 gap := trace(XZ)       = 2.36e-09
 relative gap           = 7.87e-10
 actual relative gap    = 3.68e-13
 rel. primal infeas (scaled problem)   = 6.05e-15
 rel. dual     "        "       "      = 3.36e-13
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 9.7e+01, 1.0e+00, 1.0e+00
 norm(A), norm(b), norm(C) = 7.2e+01, 2.0e+00, 7.2e+01
 Total CPU time (secs)  = 1.84  
 CPU time per iteration = 0.15  
 termination code       =  0
 DIMACS: 6.1e-15  0.0e+00  8.1e-12  0.0e+00  3.7e-13  7.9e-10
-------------------------------------------------------------------

Status: Solved
Optimal value (cvx_optval): +1.70974e-14

However, this result seems incorrect because C should be positive semidefinite for any positive real number p, and thus the maximum value of p should be 1.

What could be the cause of this issue, and how can it be fixed?

2

For p = 1, C has 1295 “zero” eigenvalues and one eigenvalue of 36. That is, if they were computed in exact arithmetic. According to eig` on my computer, 361 of these eigenvalue are actually (ever so slightly) negative (at a magnitude much smaller than what I thought would present a difficulty for a solver)

When I ran your program with Mosek and with SeDumi as solver, the optimal value of p = 1 was obtained.

Running with SDPT3, but changing the semidefinite constraint to
C + m*eye(1296)== semidefinite(4*4*9*9)
for various constant values of m (all the way up to 1) I tried, the optimal value of p was calculated as (about) 1/m. So it seems like SDPT3 is determining the minimum eigenvalue of Tensor(A, A, B, B) as (about) -1. This doesn’t seem to be just the numerical issue I thought it might be with ever so slightly negative eigenvalues causing difficulties in the solver. SDPT3 uses an infeasible path-following algorithm. I think I have read that SDPT3 might be worse at determining feasibility of some problems (or perhaps dealing with borderline feasible problems) than other solvers used by CVX. I will leave it to others to figure out whether this is a bug (in SDPT3 or the CVX/SDPT3 interface), and to make any further assessment.

I ran your original problem, converted to YALMIP, with SDPT3 as solver. It produced the optimal value of p = 1. CVX provided the (supposed) dual to SDPT3. YALMIP dualized. However, I don’t know what differences, if any, there were between the problem provided by CVX vs. YALMIP.

I suggest you use Mosek if available to you, otherwise use SeDumI.

3

Thank you for your support.

Using SeDumi, I was able to obtain the optimal value of p = 1. It seems that SDPT3 does indeed encounter numerical difficulties with certain problems. I will consider using Mosek if it is available.