Convergence to zero at certain problem sizes about scs HOT 6 CLOSED

Thanks for flagging this and giving such a detailed description. This is a very strange issue indeed. I added the following lines to your script:

data = prob.get_problem_data(cvxpy.SCS)    # get SCS problem data                                       

params = {}  # no additional params required                                                            
sol = scs.solve(data, data['dims'], **params)    # solve in SCS                                         

from genRandomConeProb import *     # cone projection routines                                          
# do the following because genRandomConeProb doesn't like missing cones                                 
data['dims']['ep'] = 0                                                                                  
data['dims']['ed'] = 0                                                                                  
data['dims']['p'] = {}                                                                                  

print('Duality gap:  {0}'.format(np.dot(data['c'], sol['x']) + np.dot(data['b'], sol['y'])))            
print('Primal residual:  {0}'.format(np.linalg.norm(data['A'].dot(sol['x']) + sol['s'] - data['b'])))   
print('Dual residual:  {0}'.format(np.linalg.norm(data['A'].T.dot(sol['y']) + data['c'])))              

print('Primal cone error:  {0}'.format(np.linalg.norm(sol['s'] - proj_cone(sol['s'], data['dims']))))   
print('Dual cone error:  {0}'.format(np.linalg.norm(sol['y'] - proj_dual_cone(sol['y'], data['dims']))))

CVXPY has a nice feature that you can pull out the data that is getting passed to the solver directly, in the canonical form it accepts. This allows us to talk directly to SCS, and check in python whether the solution it returns is actually optimal or not, given the data it received.

The genRandomConeProb file contains the routines required to project onto the primal and dual cones, and lives in the scs/examples/python directory, so you'll need to grab that to run this.

Running with those extra lines, on an example that breaks, I get:

Duality gap:  -2.0763261198597505e-05
Primal residual:  8.58649469424906e-05
Dual residual:  6.655642390908208e-05
Primal cone error:  4.664535778086328e-16
Dual cone error:  1.5716250041522723e-16

This is basically saying that the solution SCS is finding is close to optimal, to those tolerances, since the KKT conditions say that having those five quantities be zero is necessary and sufficient for global optimality.

Clearly the solution is not optimal for the problem you specified however. But this suggests that the problem is likely in cvxpy, rather than SCS, in either the problem transformation or transforming the solution back. @SteveDiamond would you mind taking a look at this to see if my hunch is right? I'm happy to keep digging into SCS if you think the problem is there.

from scs.

I can also reproduce on python 2.7.11. I thought it might have been an integer division issue.

from scs.

Thanks for taking a look. In the output from SCS, there does appear to be some difference happening in the problem generation/transformation. In the 256^2 case, there are relatively more zeros in the A matrix than in the 257^2 case or for that matter the 128^2 case (a power-of-two that exhibits no issue). At 512^2, another problem case, the same discrepancy exists. I don't recall seeing a discrepancy in the "tree-form" of the problem from CVXPY, but I would need to take another look at that. @SteveDiamond can probably track down where this is happening much faster than I can, but I will try to take a closer look soon.

from scs.

This is the most bizarre bug I've seen in CVXPY or the solvers. I have no idea what the problem is. I'm not convinced it's CVXPY, since I couldn't find anything in the canonicalization related to powers of two. Here's a simpler example that fails for the same problem sizes:

import cvxpy
SHAPE = (256, 256)
X = cvxpy.Variable(*SHAPE)
reg = cvxpy.norm(X, 'nuc')  # regularizer (trace norm)
prob = cvxpy.Problem(cvxpy.Minimize(reg), [X[0,0] >= 1])
prob.solve(solver=cvxpy.SCS, verbose=True)
print('CVXPY status:  {0}'.format(prob.status))
print('objective(X):  {0}'.format(prob.objective.value))

from scs.

That minimal example is solved for all sizes up to 256, but returns infeasible for size 256. Examining the tolerances of the certificate SCS returns, with eps = 1e-9, I get

b'y:  -1.0
Dual residual:  9.25926002537e-14
Dual cone error:  0.0

which is extremely high tolerance and clearly a real certificate that the data coming into SCS is infeasible. I dug into @jbruer hint about NNZs in A for sizes 1 up to 300, and got the following plot:

The problem at size 256 has a noticeable deviation from the trend and is the only one that fails. Note that this plot was generated without SCS, it's straight out of CVXPY:

import cvxpy                                                        
import numpy as np                                                  
import matplotlib.pyplot as plt                                     

nnzs = []                                                           
for i in range(1, 300):                                             
  SHAPE = (i, i)                                                    
  X = cvxpy.Variable(*SHAPE)                                        
  reg = cvxpy.norm(X, 'nuc')  # regularizer (trace norm)            
  prob = cvxpy.Problem(cvxpy.Minimize(reg), [X[0,0] >= 1])          

  data = prob.get_problem_data(cvxpy.SCS)    # get SCS problem data 
  nnzs.append(len(data['A'].nonzero()[0]))                          

plt.plot(nnzs)                                                      
plt.show()                                                          

I really think now that this is something in CVXPY, I have no idea what it could be though.

from scs.

Note that if you ask for the problem data using CVXOPT as the solver, you also get a similar blip in nnz(A). The actual numbers are different due to the implementation of the SDP.

I've done some more testing and have possibly narrowed down where the issue occurs. Given that this is outside of SCS, I will make an issue over at the CVXPY page with those details. Thanks, @bodono, for looking into this on your end.

from scs.