Python PICOS + CVXOPT solver objective increases without bound

[Shalom Abate]

Greetings everyone,

I've been using PICOS + CVXOPT to run some SDPs for my research in quantum information (NPA hierarchy for the CHSHp game, if anyone's curios). I have a relatively simple SDP (which I know is bounded and know the optimal value for) which I've been migrating from YALMIP (yes I was calling matlab from python) to PICOS. After getting stuck with problems in my conventional solver (Mosek) I've been trying to use cvxopt to solve the SDP. However, the solver keeps on iterating and seemingly un-boundedly increases in objective and then terminates with either division by 0 or some arithmetic error because the numbers get too big.

Here's the SDP i'm running, and the CVXOPT verbose output (nicely printed via picos):

---------------------

optimization problem  (SDP):

45 variables, 41 affine constraints, 45 vars in 1 SD cones

X : (9, 9), symmetric

maximize 0.25*( X[1,5] + X[2,6] + X[1,7] + X[2,8] + X[3,5] + X[4,6] + X[3,8] + X[4,7] )

such that

(Positivity) : X ≽ |0|

(Psi Valid State) : X[0,0] = 1.0

(Valid Probabilities (P)) : X[1,0] + X[2,0] = X[0,0]

(Valid Probabilities (Q)) : X[5,0] + X[6,0] = X[0,0]

(Valid Probabilities (P)) : X[1,1] + X[2,1] = X[0,1]

(Valid Probabilities (Q)) : X[5,1] + X[6,1] = X[0,1]

(Valid Probabilities (P)) : X[1,2] + X[2,2] = X[0,2]

(Valid Probabilities (Q)) : X[5,2] + X[6,2] = X[0,2]

(Valid Probabilities (P)) : X[1,3] + X[2,3] = X[0,3]

(Valid Probabilities (Q)) : X[5,3] + X[6,3] = X[0,3]

(Valid Probabilities (P)) : X[1,4] + X[2,4] = X[0,4]

(Valid Probabilities (Q)) : X[5,4] + X[6,4] = X[0,4]

(Valid Probabilities (P)) : X[1,5] + X[2,5] = X[0,5]

(Valid Probabilities (Q)) : X[5,5] + X[6,5] = X[0,5]

(Valid Probabilities (P)) : X[1,6] + X[2,6] = X[0,6]

(Valid Probabilities (Q)) : X[5,6] + X[6,6] = X[0,6]

(Valid Probabilities (P)) : X[1,7] + X[2,7] = X[0,7]

(Valid Probabilities (Q)) : X[5,7] + X[6,7] = X[0,7]

(Valid Probabilities (P)) : X[1,8] + X[2,8] = X[0,8]

(Valid Probabilities (Q)) : X[5,8] + X[6,8] = X[0,8]

(Valid Probabilities (P)) : X[3,0] + X[4,0] = X[0,0]

(Valid Probabilities (Q)) : X[7,0] + X[8,0] = X[0,0]

(Valid Probabilities (P)) : X[3,1] + X[4,1] = X[0,1]

(Valid Probabilities (Q)) : X[7,1] + X[8,1] = X[0,1]

(Valid Probabilities (P)) : X[3,2] + X[4,2] = X[0,2]

(Valid Probabilities (Q)) : X[7,2] + X[8,2] = X[0,2]

(Valid Probabilities (P)) : X[3,3] + X[4,3] = X[0,3]

(Valid Probabilities (Q)) : X[7,3] + X[8,3] = X[0,3]

(Valid Probabilities (P)) : X[3,4] + X[4,4] = X[0,4]

(Valid Probabilities (Q)) : X[7,4] + X[8,4] = X[0,4]

(Valid Probabilities (P)) : X[3,5] + X[4,5] = X[0,5]

(Valid Probabilities (Q)) : X[7,5] + X[8,5] = X[0,5]

(Valid Probabilities (P)) : X[3,6] + X[4,6] = X[0,6]

(Valid Probabilities (Q)) : X[7,6] + X[8,6] = X[0,6]

(Valid Probabilities (P)) : X[3,7] + X[4,7] = X[0,7]

(Valid Probabilities (Q)) : X[7,7] + X[8,7] = X[0,7]

(Valid Probabilities (P)) : X[3,8] + X[4,8] = X[0,8]

(Valid Probabilities (Q)) : X[7,8] + X[8,8] = X[0,8]

(Off Diagonals) : X[1,2] = 0

(Off Diagonals) : X[5,6] = 0

(Off Diagonals) : X[3,4] = 0

(Off Diagonals) : X[7,8] = 0

---------------------

0 secs: Solving Program...

--------------------------

  cvxopt CONELP solver

--------------------------

     pcost       dcost       gap    pres   dres   k/t

 0:  1.4066e+00  1.1562e+00  7e+01  2e+01  4e+00  1e+00

 1: -2.1177e+01  7.3801e+01  4e+05  9e+02  2e+03  3e+02

 2:  1.1882e+00  1.1340e+02  2e+02  2e+01  2e+02  1e+02

 3:  2.1549e+00  2.5880e+02  3e+02  3e+01  3e+02  3e+02

 4:  2.7772e+00  7.1004e+02  1e+03  4e+01  1e+03  7e+02

 5: -5.2925e-01  4.3954e+03  2e+04  1e+02  5e+03  4e+03

 6:  1.0399e+00  6.3467e+03  3e+04  1e+02  7e+03  6e+03

 7: -1.4600e+00  8.3463e+03  7e+04  2e+02  1e+04  8e+03

 8:  2.0511e+00  5.9767e+02  9e+03  5e+01  3e+03  3e+02

 9:  1.9139e+00  6.8854e+02  1e+05  2e+02  2e+04  3e+03

10:  4.2164e+00  4.8022e+03  2e+05  2e+02  3e+04  8e+03

... (truncated for space)...

1445:  2.1771e+13  3.4827e+22  1e+32  1e+20  6e+23  2e+22

1446:  2.1771e+13  3.4827e+22  1e+32  1e+20  5e+23  2e+22

1447:  2.1771e+13  3.4827e+22  1e+32  1e+20  5e+23  2e+22

Traceback (most recent call last):

  File "generate_picos.py", line 245, in <module>

    run_mod_p(1,2)

  File "generate_picos.py", line 240, in run_mod_p

    sol = modp.solve(solver='cvxopt', verbose=1)

  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/picos/problem.py", line 5004, in solve

    primals, duals, obj, sol = self._cvxopt_solve()

  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/picos/problem.py", line 5258, in _cvxopt_solve

    self.cvxoptVars['b'])

  File "/Users/shalomabate/Library/Python/2.7/lib/python/site-packages/cvxopt-1.1.9+1.gd8bd930-py2.7-macosx-10.6-intel.egg/cvxopt/coneprog.py", line 1395, in conelp

    misc.update_scaling(W, lmbda, ds, dz)

  File "/Users/shalomabate/Library/Python/2.7/lib/python/site-packages/cvxopt-1.1.9+1.gd8bd930-py2.7-macosx-10.6-intel.egg/cvxopt/misc.py", line 628, in update_scaling

    a = 1.0 / math.sqrt(lmbda[ind+i])

ZeroDivisionError: float division by zero

Any ideas what's going on here? Is there a known issue with CVXOPT going off the rails with a bounded SDP? I've been successfully running this same SDP with YALMIP on Matlab via mosek.

[Joachim Dahl]

Being affiliated with MOSEK I suggest you send problematic instances to support@mosek, rather than changing solver.  There are many ways to formulate a problem, including how to scale it, and you probably have chosen a bad one if MOSEK struggles with it.  Or you may, in fact, have an ill-posed problem.

--
You received this message because you are subscribed to the Google Groups "CVXOPT" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cvxopt+unsubscribe@googlegroups.com.
To post to this group, send email to cvx...@googlegroups.com.
Visit this group at https://groups.google.com/group/cvxopt.
For more options, visit https://groups.google.com/d/optout.

[Dima Pasechnik]

Shalom,

Could you post the code to generate your SDP instance somewhere?

Otherwise it's clearly not easy to say what's going wrong - could be the combo of your OS/hardware...

Thanks,

Dima

[Shalom Abate]

Hi Dima,

Would a .CBF suffice? My code for generating the SDP is pretty long.

Best,

Shalom

[Shalom Abate]

I've attached my SDP in CBF format. When I load this with PICOS and solve it with CVXOPT it goes on forever and terminates with division by 0 error.

But I can solve it fine with MOSEK and get the same result as always.

Best,

Shalom

[Dima Pasechnik]

cvxopt has a feature to use DSDP instead of the builtin conic solver. Did you try it?

[Shalom Abate]

I did not know that. I am calling 'CVXOPT' as a solver from PICOS, I'm not directly using 'CVXOPT' to create my model. Is there a way I can tell PICOS to make cvxopt use 'DSDP' instead? Some kind of option I can set?

[Martin]

It looks like it's a numerical issue with the default KKT solver in CVXOPT. Using another KKT solver seems to work:

c = prob.cvxoptVars['c']
G = prob.cvxoptVars['Gs'][0]
h = prob.cvxoptVars['hs'][0]
A = prob.cvxoptVars['A']
b = prob.cvxoptVars['b']
dims = {'l':0, 'q':[], 's':[9]}
cvxopt.solvers.conelp(c,G,h,dims,A,b,kktsolver='ldl2')

     pcost       dcost       gap    pres   dres   k/t

 0: -5.0000e-01 -5.0000e-01  2e+01  3e+00  4e+00  1e+00

 1: -6.3559e-01 -5.4390e-01  7e-01  2e-01  3e-01  2e-01

 2: -8.5233e-01 -8.4859e-01  2e-02  6e-03  7e-03  6e-03

 3: -8.5354e-01 -8.5350e-01  2e-04  6e-05  7e-05  6e-05

 4: -8.5355e-01 -8.5355e-01  2e-06  6e-07  7e-07  6e-07

 5: -8.5355e-01 -8.5355e-01  2e-08  6e-09  7e-09  6e-09

Optimal solution found.

cvxopt.solvers.conelp(c,G,h,dims,A,b,kktsolver='ldl',options={'kktreg':1e-9})

     pcost       dcost       gap    pres   dres   k/t

 0: -5.0000e-01 -5.0000e-01  2e+01  3e+00  4e+00  1e+00

 1: -6.3559e-01 -5.4390e-01  7e-01  2e-01  3e-01  2e-01

 2: -8.5233e-01 -8.4859e-01  2e-02  6e-03  7e-03  6e-03

 3: -8.5354e-01 -8.5350e-01  2e-04  6e-05  7e-05  6e-05

 4: -8.5355e-01 -8.5355e-01  2e-06  6e-07  7e-07  6e-07

 5: -8.5355e-01 -8.5355e-01  2e-08  6e-09  7e-09  6e-09

Optimal solution found.

I am not sure if it is possible to override the default CVXOPT KKT solver from PICOS.

[Shalom Abate]

Oh wow that's really interesting. I mean I'm glad to see it working.

So that means I would have to model the SDP using CVXOPT and then use another solver? I can't like pass an option down from PICOS or something...

Or if I really wanted to be hacky, I could just find the CVXOPT source and change the default solver? Is that hard to do?

[Shalom Abate]

I suppose I could also just dump the .cbf from picos and just run it in CVXOPT like you did here.

[Shalom Abate]

Alright, I just forced picos to set kktsolver='ldl2' in the source, and all is well now.

Thank you!

[Martin]

I think that the real issue is that your problem instance is numerically close to violating the "full rank condition" required by CVXOPT (i.e., either Rank(A) < p or Rank([G; A]) < n). In fact, the LDL KKT-solver throws an exception because of this. Thus, it would probably be better to go back to your model to check if you have one or more redundant constraints. Alternatively, I would recommend that you use the LDL KKT-solver with regularization, since this tends to work fairly well even when the full rank condition is violated:

cvxopt.solvers.conelp(c,G,h,dims,A,b,kktsolver='ldl',options={'kktreg':1e-9})

[Joachim Dahl]

Then using MOSEK is not an option either.  It's probably worthwhile to figure out where linear dependencies come from, and get rid of them. Not for the sake of MOSEK, but in order to have a better problem formulation and to be able solve problems faster and more reliably.

--

[Shalom Abate]

Thanks for the suggestion guys. I will go back and take a look to see if there's any redundant constraints, although I feel like that's something PICOS should detect and reject.

To unsubscribe from this group and stop receiving emails from it, send an email to cvxopt+un...@googlegroups.com.

[Joachim Dahl]

Linear dependencies cannot in general be removed reliably (that's an NP-hard problem).  It's really up to the user to formulate problems without redudant constrants, and with somewhat sensible scaling.  If PICOS introduces linear dependencies in a systematic way, then that's very bad,  but I doubt that's the case. It could also be that your data are so badly scaled that ADA' is numerically rank-deficient,  but that's just as a bad, and not something the solver can always fix automatically.

MOSEK employs a number of heuristics to detect and remove *some* dependencies, but there are no guarantees all dependencies are detected,  and when not all dependencies are removed then usually the solver runs into numerical problems.  It's also a very time-consuming step for large models, so often users prefer to disable the dependency-checker. 

To unsubscribe from this group and stop receiving emails from it, send an email to cvxopt+unsubscribe@googlegroups.com.