cvxgrp / scs

Splitting Conic Solver
MIT License
553 stars 136 forks source link

SCS returns wrong result #156

Closed astghikhakobyan closed 3 years ago

astghikhakobyan commented 3 years ago

Hello,

I am using SCS for solving an SDP problem in cvxpy. However, the result I get is incorrect when I compare it to other solvers such as MOSEK, SEDUMI, and SDPT3. All three solvers return somewhat similar results, while the SCS one is completely different. Similarly, when I solve the dual version of my problem, in the case of MOSEK the optimal values are very close, while in the case of SCS completely different. Please see below my python code.

import cvxpy as cp import numpy as np

ny = 2 p = 3 theta = 0.0001 alpha = 0.95 r = 0.2

Z = cp.Variable((ny,ny), PSD=True) z = cp.Variable((1,1)) tau = cp.Variable((1,1)) lambda1 = cp.Variable() eps = cp.Variable((1,1)) gamma = cp.Variable((ny,1)) Gamma = cp.Variable((ny,ny), symmetric=True) mu_hat = cp.Parameter(ny) Sigma_hat = cp.Parameter(ny+1) y = cp.Parameter(ny) mu_hat1 = cp.Parameter(1) Sigma_hat1 = cp.Parameter(ny+1) y1 = cp.Parameter(1)

Sigma = cp.Expression() Sigma = cp.bmat([[Sigma_hat[0], Sigma_hat[1]],[Sigma_hat[1], Sigma_hat[2]]]) Sigma1 = cp.Expression() Sigma1 = cp.bmat([[Sigma_hat1[0], Sigma_hat1[1]],[Sigma_hat1[1], Sigma_hat1[2]]])

constraints = [ cp.bmat([[lambda1 np.eye(ny) - Gamma, gamma + cp.vstack(lambda1 mu_hat)], [(gamma + cp.vstack(lambda1 mu_hat)).T, eps] ]) >> 0, cp.bmat([[lambda1 np.eye(ny) - Gamma, lambda1 Sigma], [(lambda1 Sigma).T, Z] ]) >> 0, cp.bmat([[np.eye(ny) + Gamma, gamma - cp.vstack(y)], [(gamma - cp.vstack(y)).T, tau + z] ]) >> 0, cp.bmat([[Gamma, gamma], [gamma.T, tau] ]) >> 0, eps >= 0, lambda1 >= 0 ] obj = cp.Minimize(- y1 + r r + z + 1/(1 - alpha) (tau + eps + cp.trace(Z) + lambda1 (theta theta - mu_hat1 - cp.trace(Sigma1))) ) prob = cp.Problem(obj, constraints)

In = np.array([7.09364830858073, 5.36952559622012, 0.305480348185421, 0.015403254279958, 0.639418852045523, 8.54820909414561, 6.36110970560502]) mu_hat.value = In[5:7] Sigma_hat.value = In[2:5] y.value = In[0:2] mu_hat1.value = np.array([In[5:7].T @ In[5:7]]) Sigma_hat1.value = np.array([In[2]In[2]+In[3]In[3], In[2]In[3]+In[3]In[4], In[3]In[3]+In[4]In[4]]) y1.value = np.array([In[0:2].T @ In[0:2]]) prob.solve(solver = cp.SCS, verbose = False) print("Optimal Value (SCS): ", prob.value, "solver time: ", prob.solver_stats.solve_time) prob.solve(solver = cp.MOSEK, verbose = False) print("Optimal Value (MOSEK): ", prob.value, "solver time: ", prob.solver_stats.solve_time)

In the cvxpy repository I was advised to increase the solver precision, but that didn't really help.

bodono commented 3 years ago

When I run your code (with verbose=True) I get:

----------------------------------------------------------------------------
    SCS v2.1.2 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 40
eps = 1.00e-04, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 12, constraints m = 33
Cones:  linear vars: 2
    sd vars: 31, sd blks: 5
Setup time: 2.67e-04s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 3.41e+19  2.50e+19  9.87e-01  5.89e+20  8.98e+22  4.83e+22  7.68e-04
   100| 3.66e-03  5.93e-03  4.13e-02  7.94e+01  8.63e+01  9.05e-14  3.40e-03
   200| 3.71e-04  2.91e-04  8.44e-04  7.68e+01  7.67e+01  3.39e-13  5.14e-03
   300| 8.22e-04  9.25e-04  4.90e-03  7.70e+01  7.78e+01  4.75e-14  6.67e-03
   400| 5.32e-04  7.56e-04  1.64e-03  7.69e+01  7.67e+01  2.68e-14  8.17e-03
   500| 2.56e-04  6.32e-04  2.32e-03  7.62e+01  7.66e+01  2.39e-13  9.69e-03
   600| 2.69e+18  5.12e+17  9.36e-01 -5.08e+19 -1.69e+18  6.37e+20  1.12e-02
   700| 4.53e-05  3.57e-05  7.74e-04  7.94e+01  7.96e+01  1.03e-13  1.28e-02
   760| 9.86e-06  2.97e-06  1.79e-07  7.91e+01  7.91e+01  1.30e-13  1.37e-02
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.37e-02s
    Lin-sys: nnz in L factor: 87, avg solve time: 3.41e-07s
    Cones: avg projection time: 1.31e-05s
    Acceleration: avg step time: 3.52e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 2.1514e-09, dist(y, K*) = 2.2195e-09, s'y/|s||y| = -1.2370e-11
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 9.8637e-06
dual res:   |A'y + c|_2 / (1 + |c|_2) = 2.9723e-06
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 1.7878e-07
----------------------------------------------------------------------------
c'x = 79.1368, -b'y = 79.1369
============================================================================
Optimal Value (SCS):  0.025195336215986686 solver time:  13.736438

I don't have a MOSEK license so I can't run the rest, but it looks to me like it's solving the problem fine. I can run it with CVXOPT, which gets:

     pcost       dcost       gap    pres   dres   k/t
 0:  9.4227e+02  7.3368e+02  1e+04  2e+00  3e-01  1e+00
 1:  1.0471e+02  4.5939e+01  3e+04  1e+00  2e-01  1e+02
 2:  8.2810e+02  1.6318e+03  5e+04  1e+00  2e-01  9e+02
 3:  1.1637e+03  1.7584e+03  2e+04  5e-01  8e-02  7e+02
 4: -3.4335e+02  1.3772e+02  1e+04  3e-01  5e-02  5e+02
 5: -3.4805e+01  4.9283e+01  2e+03  6e-02  1e-02  9e+01
 6:  6.2995e+01  8.0625e+01  3e+02  1e-02  2e-03  2e+01
 7:  7.6330e+01  7.8702e+01  4e+01  1e-03  2e-04  3e+00
 8:  7.8529e+01  7.9122e+01  1e+01  4e-04  6e-05  6e-01
 9:  7.9010e+01  7.9098e+01  2e+00  5e-05  9e-06  9e-02
10:  7.9124e+01  7.9136e+01  2e-01  7e-06  1e-06  1e-02
11:  7.9132e+01  7.9137e+01  9e-02  3e-06  4e-07  6e-03
12:  7.9135e+01  7.9137e+01  2e-02  6e-07  9e-08  2e-03
13:  7.9131e+01  7.9132e+01  1e-02  2e-07  4e-08  1e-03
14:  7.9129e+01  7.9130e+01  8e-03  1e-07  2e-08  6e-04
15:  7.9128e+01  7.9129e+01  2e-02  2e-07  3e-08  1e-03
16:  7.9128e+01  7.9128e+01  1e-02  1e-07  2e-08  6e-04
17:  7.9125e+01  7.9125e+01  4e-03  3e-08  4e-09  3e-04
18:  7.9124e+01  7.9124e+01  5e-03  2e-08  3e-09  3e-04
19:  7.9123e+01  7.9123e+01  2e-03  5e-09  7e-10  1e-04
20:  7.9122e+01  7.9123e+01  2e-03  1e-07  1e-07  1e-04
21:  7.9122e+01  7.9123e+01  5e-05  1e-07  3e-08  3e-06
Terminated (singular KKT matrix).

Unfortunately CVXOPT fails as it runs into numerical difficulties, but the objective value it is converging to is roughly 79.122, and SCS is getting 79.1368, so they are very close and looks like SCS is solving it correctly. What's the issue?

bodono commented 3 years ago

In reduced eps to 1e-6 and increased max_iters, and SCS solves it as

 11800| 3.09e-06  9.38e-06  1.62e-05  7.91e+01  7.91e+01  2.22e-13  2.00e-01
 11900| 6.20e-07  8.89e-07  3.82e-06  7.91e+01  7.91e+01  6.70e-14  2.02e-01
 12000| 6.74e-06  7.22e-06  5.26e-06  7.91e+01  7.91e+01  1.00e-13  2.03e-01
 12100| 1.82e-07  1.99e-07  7.32e-07  7.91e+01  7.91e+01  8.46e-14  2.05e-01
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 2.05e-01s
    Lin-sys: nnz in L factor: 87, avg solve time: 2.97e-07s
    Cones: avg projection time: 1.25e-05s
    Acceleration: avg step time: 3.19e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 1.3223e-07, dist(y, K*) = 2.3510e-09, s'y/|s||y| = -2.2286e-12
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 1.8165e-07
dual res:   |A'y + c|_2 / (1 + |c|_2) = 1.9936e-07
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 7.3195e-07
----------------------------------------------------------------------------
c'x = 79.1226, -b'y = 79.1227
============================================================================

Looks to me like it's working just fine, even at the higher precision.

astghikhakobyan commented 3 years ago

You are right, the solution becomes closer to the optimal for eps = 1e-6, but for example, when solving the same problem for eps=1e-7, I get

17100| 2.27e-06  4.05e-06  8.30e-07  7.91e+01  7.91e+01  1.42e-13  2.33e-01 
 17200| 3.24e-07  1.46e-07  1.92e-07  7.91e+01  7.91e+01  2.40e-13  2.34e-01 
 17300| 2.48e-06  2.14e-06  3.10e-06  7.91e+01  7.91e+01  1.02e-13  2.36e-01 
 17400| 9.56e-08  5.32e-08  1.16e-07  7.91e+01  7.91e+01  5.37e-14  2.37e-01 
 17440| 3.31e-08  4.35e-08  2.16e-08  7.91e+01  7.91e+01  3.63e-14  2.37e-01 
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 2.37e-01s
    Lin-sys: nnz in L factor: 87, avg solve time: 2.33e-07s
    Cones: avg projection time: 9.55e-06s
    Acceleration: avg step time: 2.45e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 9.8424e-08, dist(y, K*) = 1.9700e-09, s'y/|s||y| = 1.0673e-11
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 3.3093e-08
dual res:   |A'y + c|_2 / (1 + |c|_2) = 4.3497e-08
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 2.1630e-08
----------------------------------------------------------------------------
c'x = 79.1216, -b'y = 79.1216
============================================================================

with prob.value = 0.00992025968734822 and solver_time = 237.341735. Meanwhile MOSEK returns the following:

19  5.3e-11  1.7e-07  1.3e-15  6.83e-01   7.912130088e+01   7.912129381e+01   1.5e-13  0.02  
20  2.1e-11  2.2e-07  5.3e-16  4.25e-01   7.912126443e+01   7.912125853e+01   6.1e-14  0.02  
21  3.9e-12  5.5e-06  1.9e-17  8.39e-01   7.912123626e+01   7.912123494e+01   1.1e-14  0.02  
22  6.9e-13  1.8e-05  8.3e-18  8.02e-01   7.912122816e+01   7.912122783e+01   2.0e-15  0.02  
Optimizer terminated. Time: 0.02    
`

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 7.9121228159e+01    nrm: 2e+03    Viol.  con: 6e-08    var: 2e-11    barvar: 0e+00  
  Dual.    obj: 7.9121227834e+01    nrm: 4e+04    Viol.  con: 0e+00    var: 0e+00    barvar: 1e-10  

with prob.value = 0.009576704754323373 and solver_time = 0.02266097068786621. (Note that other solvers such as SEDUMI and SDPT3 return solutions very close to the MOSEK solutions in very little iterations)

This gap is more critical when I solve the problem for different parameters, e.g.

In = np.array([9.70592781760616, 6.49677170169942, 0.254647324863969, 0.233294850086384, 0.524667637294570, 1.47366805213582, 7.69765550241156])

Then the solution returned by SCS in 5 milion iterations looks like the following:

4999300| 1.03e-06  5.75e-07  7.00e-06  8.13e+01  8.13e+01  5.27e-14  5.07e+01 
4999400| 6.47e-07  4.48e-07  1.68e-06  8.13e+01  8.13e+01  5.27e-14  5.07e+01 
4999500| 9.82e-07  4.52e-07  2.40e-06  8.13e+01  8.13e+01  5.33e-14  5.07e+01 
4999600| 8.87e-07  4.26e-07  1.10e-06  8.13e+01  8.13e+01  1.59e-13  5.07e+01 
4999700| 6.63e-07  6.62e-07  4.34e-07  8.13e+01  8.13e+01  1.59e-13  5.07e+01 
4999800| 8.84e-07  4.23e-07  1.95e-06  8.13e+01  8.13e+01  5.33e-14  5.07e+01 
4999900| 9.54e-07  5.21e-07  7.02e-06  8.13e+01  8.13e+01  5.33e-14  5.07e+01 
5000000| 8.06e-07  5.03e-07  2.56e-06  8.13e+01  8.13e+01  5.27e-14  5.07e+01 
----------------------------------------------------------------------------
Status: Solved/Inaccurate
Hit max_iters, solution may be inaccurate, returning best found solution.
Timing: Solve time: 5.07e+01s
    Lin-sys: nnz in L factor: 87, avg solve time: 2.16e-07s
    Cones: avg projection time: 8.66e-06s
    Acceleration: avg step time: 1.65e-08s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 0.0000e+00, dist(y, K*) = 2.4399e-09, s'y/|s||y| = 6.1715e-11
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 6.0611e-07
dual res:   |A'y + c|_2 / (1 + |c|_2) = 5.2853e-07
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 3.6760e-07
----------------------------------------------------------------------------
c'x = 81.2564, -b'y = 81.2564

While MOSEK solves in only 28 iterations:

21  1.4e-10  9.1e-09  7.9e-14  6.57e-01   8.120551330e+01   8.120465588e+01   9.0e-13  0.02  
22  5.3e-11  2.2e-10  3.0e-14  3.18e-01   8.119968437e+01   8.119878374e+01   3.4e-13  0.02  
23  1.4e-11  4.9e-08  4.9e-15  7.02e-01   8.119531511e+01   8.119497079e+01   8.9e-14  0.02  
24  6.2e-12  7.7e-08  2.1e-15  4.37e-01   8.119294233e+01   8.119266220e+01   3.0e-14  0.02  
25  6.2e-12  5.8e-09  6.3e-17  7.82e-01   8.119152657e+01   8.119143969e+01   7.0e-15  0.02  
26  9.5e-12  1.9e-07  1.3e-16  6.98e-01   8.119100142e+01   8.119096348e+01   2.0e-15  0.02  
27  7.3e-12  2.2e-07  3.8e-17  9.56e-01   8.119077906e+01   8.119077639e+01   1.3e-16  0.02  
28  2.9e-12  6.7e-04  1.1e-17  9.86e-01   8.119076317e+01   8.119076292e+01   1.2e-17  0.02  
Optimizer terminated. Time: 0.02    

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 8.1190763170e+01    nrm: 1e+03    Viol.  con: 4e-07    var: 1e-12    barvar: 0e+00  
  Dual.    obj: 8.1190762930e+01    nrm: 1e+06    Viol.  con: 0e+00    var: 4e-12    barvar: 1e-09  

And still, the gap is too big between the two solutions.

bodono commented 3 years ago

Firstly, it's meaningless to compare the number of iterations. Mosek / sedumi / sdpt3 are all second-order solvers, SCS is a first-order solver. Secondly I don't see how these results are that different, note that mosek is returning solutions with constraint violations of about 6e-8 in the first example and 4e-7 in the second, which is about the same as SCS in both. The prob.value is from cvxpy, not what the solvers are returning. You can compare what the solvers return directly, for the first problem SCS gets 79.1216 and mosek gets 79.1212, both with high accuracy (O(1e-8)) so it's not clear which is closer to the true value. Similarly for the second case SCS gets 81.2564 and mosek gets 81.1907, again both with high accuracy (O(1e-7)). If you want high accuracy and mosek is working for your problem, then you should use that. SCS is a first-order solver so is better suited for moderate accuracy solutions for very large problems.