SCS Fails on Quad Program from OSQP Benchmark (LASSO)

[RoyiAvital]

Specifications

• OS: Windows 10
• SCS Version: 2.1.4 (Python PIP)
• Compiler: Python

Description

I have a pretty simple Quad Program which SCS fails to solve (Returns nan) while Gurobi / OSQP manages to solver it.
I am using CVXPY to utilize SCS

How to reproduce

Python code to reproduce:

# General tools
import numpy as np
import scipy as sp

# Sparse
import scipy.sparse as spa

# Solvers
import cvxpy

# Code from OSQP Benchmark
class LassoExample(object):
    '''
    Lasso QP example
    '''
    def __init__(self, n, seed=1):
        # Set random seed
        np.random.seed(seed)

        self.n = int(n)               # Number of features
        self.m = int(self.n * 10)    # Number of data-points

        self.Ad = spa.random(self.m, self.n, density=0.15,
                             data_rvs=np.random.randn)
        self.x_true = np.multiply((np.random.rand(self.n) >
                                   0.5).astype(float),
                                  np.random.randn(self.n)) / np.sqrt(self.n)
        self.bd = self.Ad.dot(self.x_true) + np.random.randn(self.m)
        self.lambda_max = np.linalg.norm(self.Ad.T.dot(self.bd), np.inf)
        self.lambda_param = (1./5.) * self.lambda_max

        self.qp_problem = self._generate_qp_problem()

    def _generate_qp_problem(self):
        # Construct the problem
        #       minimize	y' * y + lambda * 1' * t
        #       subject to  y = Ax - b
        #                   -t <= x <= t
        P = spa.block_diag((spa.csc_matrix((self.n, self.n)), 2*spa.eye(self.m), spa.csc_matrix((self.n, self.n))), format='csc')
        q = np.append(np.zeros(self.m + self.n), self.lambda_param * np.ones(self.n))
        In = spa.eye(self.n)
        Onm = spa.csc_matrix((self.n, self.m))
        A = spa.vstack([spa.hstack([self.Ad, -spa.eye(self.m),
                                    spa.csc_matrix((self.m, self.n))]),
                        spa.hstack([In, Onm, -In]),
                        spa.hstack([In, Onm, In])]).tocsc()
        l = np.hstack([self.bd, -np.inf * np.ones(self.n), np.zeros(self.n)])
        u = np.hstack([self.bd, np.zeros(self.n), np.inf * np.ones(self.n)])

        return P, q, A, l, u, A.shape[0], A.shape[1]

lassoQp = LassoExample(3);
P, q, A, l, u, m, n = lassoQp._generate_qp_problem();


x = cvxpy.Variable(lassoQp.n + lassoQp.m + lassoQp.n);

cvxObj  = cvxpy.Minimize(cvxpy.quad_form(x, P) + (x @ q));
cvxCons = [l <= A @ x, A @ x <= u];

cvxPrblm    = cvxpy.Problem(cvxObj, cvxCons);
cvxRes      = cvxPrblm.solve(solver = cvxpy.SCS, verbose = True);

I am also attaching as a .txt file SCSFail.txt.

Additional information

It is a LASSO problem converted into QP.
Actually SCS fails other problems from OSQP Benchmark such as: Huber Fitting, SVM.

Running the same problems with OSQP or Gurobi yields a valid solution.

Output

===============================================================================
                                     CVXPY                                     
                                    v1.1.15                                    
===============================================================================
(CVXPY) Aug 24 03:29:23 PM: Your problem has 36 variables, 2 constraints, and 0 parameters.
(CVXPY) Aug 24 03:29:23 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Aug 24 03:29:23 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Aug 24 03:29:23 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
                                  Compilation                                  
-------------------------------------------------------------------------------
(CVXPY) Aug 24 03:29:23 PM: Compiling problem (target solver=SCS).
(CVXPY) Aug 24 03:29:23 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS
(CVXPY) Aug 24 03:29:23 PM: Applying reduction Dcp2Cone
(CVXPY) Aug 24 03:29:23 PM: Applying reduction CvxAttr2Constr
(CVXPY) Aug 24 03:29:23 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Aug 24 03:29:23 PM: Applying reduction SCS
(CVXPY) Aug 24 03:29:23 PM: Finished problem compilation (took 1.600e-02 seconds).
-------------------------------------------------------------------------------
                                Numerical solver                               
-------------------------------------------------------------------------------
(CVXPY) Aug 24 03:29:23 PM: Invoking solver SCS  to obtain a solution.
----------------------------------------------------------------------------
	SCS v2.1.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 144
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 = 37, constraints m = 104
Cones:	linear vars: 72
	soc vars: 32, soc blks: 1
Setup time: 1.03e-03s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0|      nan       nan       nan       nan -nan(ind)       nan  1.01e-03 
   100|      nan       nan       nan       nan -nan(ind)       nan  4.29e-03 
   200|      nan       nan       nan       nan -nan(ind)       nan  7.17e-03 
   300|      nan       nan       nan       nan -nan(ind)       nan  8.96e-03 
   400|      nan       nan       nan       nan -nan(ind)       nan  9.89e-03 
   500|      nan       nan       nan       nan -nan(ind)       nan  1.13e-02 
   600|      nan       nan       nan       nan -nan(ind)       nan  1.37e-02 
   700|      nan       nan       nan       nan -nan(ind)       nan  1.50e-02 
   800|      nan       nan       nan       nan -nan(ind)       nan  1.64e-02 
   900|      nan       nan       nan       nan -nan(ind)       nan  1.76e-02 
  1000|      nan       nan       nan       nan -nan(ind)       nan  1.88e-02 
  1100|      nan       nan       nan       nan -nan(ind)       nan  1.99e-02 
  1200|      nan       nan       nan       nan -nan(ind)       nan  2.16e-02 
  1300|      nan       nan       nan       nan -nan(ind)       nan  2.35e-02 
  1400|      nan       nan       nan       nan -nan(ind)       nan  2.49e-02 
  1500|      nan       nan       nan       nan -nan(ind)       nan  2.60e-02 
  1600|      nan       nan       nan       nan -nan(ind)       nan  2.73e-02 
  1700|      nan       nan       nan       nan -nan(ind)       nan  2.87e-02 
  1800|      nan       nan       nan       nan -nan(ind)       nan  2.97e-02 
  1900|      nan       nan       nan       nan -nan(ind)       nan  3.08e-02 
  2000|      nan       nan       nan       nan -nan(ind)       nan  3.17e-02 
  2100|      nan       nan       nan       nan -nan(ind)       nan  3.31e-02 
  2200|      nan       nan       nan       nan -nan(ind)       nan  3.48e-02 
  2300|      nan       nan       nan       nan -nan(ind)       nan  3.65e-02 
  2400|      nan       nan       nan       nan -nan(ind)       nan  3.75e-02 
  2500|      nan       nan       nan       nan -nan(ind)       nan  3.85e-02 
  2600|      nan       nan       nan       nan -nan(ind)       nan  3.99e-02 
  2700|      nan       nan       nan       nan -nan(ind)       nan  4.18e-02 
  2800|      nan       nan       nan       nan -nan(ind)       nan  4.36e-02 
  2900|      nan       nan       nan       nan -nan(ind)       nan  4.54e-02 
  3000|      nan       nan       nan       nan -nan(ind)       nan  4.68e-02 
  3100|      nan       nan       nan       nan -nan(ind)       nan  4.79e-02 
  3200|      nan       nan       nan       nan -nan(ind)       nan  4.99e-02 
  3300|      nan       nan       nan       nan -nan(ind)       nan  5.10e-02 
  3400|      nan       nan       nan       nan -nan(ind)       nan  5.20e-02 
  3500|      nan       nan       nan       nan -nan(ind)       nan  5.32e-02 
  3600|      nan       nan       nan       nan -nan(ind)       nan  5.42e-02 
  3700|      nan       nan       nan       nan -nan(ind)       nan  5.51e-02 
  3800|      nan       nan       nan       nan -nan(ind)       nan  5.62e-02 
  3900|      nan       nan       nan       nan -nan(ind)       nan  5.71e-02 
  4000|      nan       nan       nan       nan -nan(ind)       nan  5.82e-02 
  4100|      nan       nan       nan       nan -nan(ind)       nan  5.92e-02 
  4200|      nan       nan       nan       nan -nan(ind)       nan  6.02e-02 
  4300|      nan       nan       nan       nan -nan(ind)       nan  6.14e-02 
  4400|      nan       nan       nan       nan -nan(ind)       nan  6.23e-02 
  4500|      nan       nan       nan       nan -nan(ind)       nan  6.33e-02 
  4600|      nan       nan       nan       nan -nan(ind)       nan  6.42e-02 
  4700|      nan       nan       nan       nan -nan(ind)       nan  6.52e-02 
  4800|      nan       nan       nan       nan -nan(ind)       nan  6.63e-02 
  4900|      nan       nan       nan       nan -nan(ind)       nan  6.72e-02 
  5000|      nan       nan       nan       nan -nan(ind)       nan  6.81e-02 
----------------------------------------------------------------------------
Status: Unbounded/Inaccurate
Hit max_iters, solution may be inaccurate, returning best found solution.
Timing: Solve time: 6.81e-02s
	Lin-sys: nnz in L factor: 303, avg solve time: 6.89e-07s
	Cones: avg projection time: 1.12e-07s
	Acceleration: avg step time: 9.24e-06s
----------------------------------------------------------------------------
Certificate of dual infeasibility:
dist(s, K) = 0.0000e+00
|Ax + s|_2 * |c|_2 = -nan(ind)
c'x = nan
============================================================================
-------------------------------------------------------------------------------
                                    Summary                                    
-------------------------------------------------------------------------------
(CVXPY) Aug 24 03:29:23 PM: Problem status: unbounded_inaccurate
(CVXPY) Aug 24 03:29:23 PM: Optimal value: -inf
(CVXPY) Aug 24 03:29:23 PM: Compilation took 1.600e-02 seconds
(CVXPY) Aug 24 03:29:23 PM: Solver (including time spent in interface) took 7.002e-02 seconds
<Path>\site-packages\cvxpy\problems\problem.py:1278: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.

[bodono]

Thanks for posting this in a way that made it easy to reproduce!

Turns out that the data has inf values in several places, which propagate to nans inside SCS (and also breaks ECOS it turns out), so it seems like an issue with CVXPY. To see this simply modify these lines and inspect the data object here:

cvxPrblm    = cvxpy.Problem(cvxObj, cvxCons);
data = cvxPrblm.get_problem_data(cvxpy.SCS);
import pdb;pdb.set_trace()

In SCS 3.0 we can support inf values in the new 'box' cone, which is the 'natural' way to handle them. I hacked SCS 3.0 into your script as follows:

# General tools
import numpy as np
import scipy as sp

# Sparse
import scipy.sparse as spa

# Solvers
import cvxpy
import scs

# Code from OSQP Benchmark
class LassoExample(object):
    '''
    Lasso QP example
    '''
    def __init__(self, n, seed=1):
        # Set random seed
        np.random.seed(seed)

        self.n = int(n)               # Number of features
        self.m = int(self.n * 10)    # Number of data-points

        self.Ad = spa.random(self.m, self.n, density=0.15,
                             data_rvs=np.random.randn)
        self.x_true = np.multiply((np.random.rand(self.n) >
                                   0.5).astype(float),
                                  np.random.randn(self.n)) / np.sqrt(self.n)
        self.bd = self.Ad.dot(self.x_true) + np.random.randn(self.m)
        self.lambda_max = np.linalg.norm(self.Ad.T.dot(self.bd), np.inf)
        self.lambda_param = (1./5.) * self.lambda_max

        self.qp_problem = self._generate_qp_problem()

    def _generate_qp_problem(self):
        # Construct the problem
        #       minimize    y' * y + lambda * 1' * t
        #       subject to  y = Ax - b
        #                   -t <= x <= t
        P = spa.block_diag((spa.csc_matrix((self.n, self.n)), 2*spa.eye(self.m),
spa.csc_matrix((self.n, self.n))), format='csc')
        q = np.append(np.zeros(self.m + self.n), self.lambda_param *
np.ones(self.n))
        In = spa.eye(self.n)
        Onm = spa.csc_matrix((self.n, self.m))
        A = spa.vstack([spa.hstack([self.Ad, -spa.eye(self.m),
                                    spa.csc_matrix((self.m, self.n))]),
                        spa.hstack([In, Onm, -In]),
                        spa.hstack([In, Onm, In])]).tocsc()
        l = np.hstack([self.bd, -np.inf * np.ones(self.n), np.zeros(self.n)])
        u = np.hstack([self.bd, np.zeros(self.n), np.inf * np.ones(self.n)])

        return P, q, A, l, u, A.shape[0], A.shape[1]

lassoQp = LassoExample(3);
P, q, A, l, u, m, n = lassoQp._generate_qp_problem();

# Hack out the equality constraints
idxs = (u - l < 1e-6)
idxs &= (u < 1e20)
idxs &= (l > -1e20)
o_idxs = np.array(range(A.shape[0])) # no need +1 for box cone
o_idxs = np.hstack((o_idxs[idxs], o_idxs[~idxs]))
inv_perm = np.argsort(o_idxs)

A_scs = spa.vstack((A[idxs, :], np.zeros((1, n)), -A[~idxs, :]))
b_scs = np.hstack((u[idxs], 1, np.zeros(m - np.sum(idxs))))

data = dict(P=spa.csc_matrix(P), c=q,
            A=spa.csc_matrix(A_scs), b=b_scs)
cone = dict(z=np.int(np.sum(idxs)), bl=l[~idxs].tolist(), bu=u[~idxs].tolist())

scs.solve(data, cone, verbose=True)

With result:

=> python tmp.py
------------------------------------------------------------------
	       SCS v3.0.0 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 36, constraints m: 37
cones: 	  z: primal zero / dual free vars: 30
	  b: box cone vars: 7
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, warm_start: 0
lin-sys:  sparse-direct
	  nnz(A): 56, nnz(P): 30
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 9.39e+00  5.44e-01  6.36e+01  6.62e+01  1.00e-01  9.72e-05
    50| 2.20e-06  3.49e-07  4.90e-07  2.14e+01  1.00e-01  1.80e-04
------------------------------------------------------------------
status:  solved
timings: total: 4.09e-04s = setup: 2.25e-04s + solve: 1.84e-04s
	 lin-sys: 3.64e-05s, cones: 8.49e-06s, accel: 1.31e-06s
------------------------------------------------------------------
cones: dist(s, K) = 0.00e+00, dist(y, K*) = 0.00e+00
comp slack: s'y/|s||y| = 0.00e+00, gap: |x'Px+c'x+b'y| = 4.90e-07
pri res: |Ax+s-b| = 2.20e-06, dua res: |Px+A'y+c| = 3.49e-07
------------------------------------------------------------------
objective = 21.406531
------------------------------------------------------------------

[RoyiAvital]

Indeed it fails on ECOS as well as I posted in embotech/ecos#200.