bodono
commented
3 years ago 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
Specifications
2.1.4(Python PIP)Description
I have a pretty simple Quad Program which SCS fails to solve (Returns
nan) while Gurobi / OSQP manages to solver it. I am usingCVXPYto utilizeSCSHow to reproduce
Python code to reproduce:
I am also attaching as a
.txtfile 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