SCS does not fully satisfy constraints?

[h-vetinari]

Hi!

Not sure if this issue should also go in https://github.com/bodono/scs-python/ and/or https://github.com/cvxgrp/cvxpy, but I think it's best to start here.

I'm working from python, and managed to upgrade everything to the most current versions (which is not easy on Windows+Python3.6, but possible thanks to the wheels from https://www.lfd.uci.edu/~gohlke/pythonlibs/)

My issue is that the constraints I pass to SCS (admittedly through cvxpy) are not respected by the solver. I'm not talking about > vs >=, but actual violations, and not on the order of machine precision.

I tried to pare down the following example as much as possible - it is just a toy example, of course. In this toy example, I can add a penalty term for the constraints and get things to work, but in my real-life application, I have to do a parameter-sweep of 3 separate parameters to hopefully get one result that satisfies my constraints. This is not a very satisfying (or scalable) situation, of course.

So the desired outcome would be clearly: fully respect constraints - e.g. x >=0 must rule out any element having x<0.

import numpy as np # 1.15.0rc2+mkl
import pandas as pd # 0.23.3
import cvxpy as cvx # 1.0.6, scs: 2.0.2, ecos: 2.0.4

scale = 10 ** 3
nrow = 100
ncol = 10
A = pd.DataFrame(np.random.uniform(0, scale, (nrow, ncol)))

# normally, the target p is given; here we try to generate an example
tmp = pd.DataFrame(np.random.uniform(scale, scale*2, (nrow, ncol)))
weight = pd.Series(np.random.uniform(10, 50, (ncol,)))
p = tmp.dot(weight)

# to interpret x as a weight, we normalize the columns
# (also helps making problem feasible; plus other reasons)
A_nrm = A.apply(lambda x: x.divide(x.sum()), axis = 0)

# parameters needed to tune for ACTUALLY achieving the constraints
augm = 1000
maxit = 1000
safety = 0.1
inflate = 1
A_nrm = A_nrm.apply(lambda x: x * inflate, axis = 0)

x = cvx.Variable(len(A_nrm.columns))
objective = cvx.Minimize(cvx.sum(cvx.abs(A_nrm.values * x - p.values))
                         ## with this penalty term, constraints work out
                         # + augm ** 2 * cvx.sum(cvx.abs(cvx.neg(x)))
                         )
constraints = [x >= 0 + safety, A_nrm.values * x <= p.values]
prob = cvx.Problem(objective, constraints)
prob.solve(solver = 'SCS', max_iters = maxit)

# from here on, just "reporting"
if x.value is None:
    print(f'Status: {prob.status}')
x = pd.Series(x.value, index = A_nrm.columns)
if (x < 0).any():
    print('x<0')
    tmp = x.loc[x<0]
    tmp_rel = (tmp/x.abs().max()).apply(lambda x: f'{x:.2e}')
    print(pd.concat([tmp, tmp_rel], axis=1, keys=['error', 'rel2max']))
q = A_nrm.dot(x)
z = p - q
if (z < 0).any():
    print('A*x > p')
    tmp = z.loc[z<0]
    tmp_rel = (tmp/z.abs().max()).apply(lambda x: f'{x:.2e}')
    print(pd.concat([tmp, tmp_rel], axis=1, keys=['error', 'rel2max']))
rsv = q.sum()/p.sum()
print(f'Resolved {rsv*100:.2f}%')

Output would be something like:

x<0
       error    rel2max
3 -12.995466  -2.24e-06
A*x > p
       error    rel2max
25 -0.497438  -1.49e-06
28 -0.204003  -6.10e-07
54 -0.010426  -3.12e-08
55 -0.054533  -1.63e-07
80 -0.179385  -5.37e-07
87 -0.911080  -2.73e-06
Resolved 77.30%

Some further notes:

• the ECOS solver does enforce the constraints, but often fails; I could solve more situations overall with SCS + the parameter sweep I mentioned; Plus ECOS had yielded weird results (may have improved with the new version I downloaded to open this issue).
• upgrading from SCS=1.2.7 (the last version currently on https://www.lfd.uci.edu/~gohlke/pythonlibs/) to 2.0.2, reduced the above (relative) errors by about 2 orders of magnitude.

[bodono]

Thanks for reporting such a detailed problem. The issue here is that SCS is a first-order solver and as such is good at getting moderate accuracy quickly, but possibly taking a long time to get high accuracy. In theory it will converge to any accuracy you desire, so to change the default accuracy (which I think is 1e-5) just pass eps=1e-8 or whatever to cvxpy when calling prob.solve. Also if you pass verbose=True then you will get more print out about what SCS and ECOS are actually doing.

However, you might not get the constraints satisfied exactly with any solver, since all of them will have some small tolerance. Your best bet would be to round the solution, or just specify a slight tolerance in your problem specification if this is critical to your use-case, for example rather than x>=0 have x>=1e-5 and have eps=1e-9 for example.

[h-vetinari]

@bodono Thanks for taking the time to answer!

In theory it will converge to any accuracy you desire

I found that the increase in iterations did not change anything, but I'll very happily try if I can get there just by playing with eps. I could certainly live with eps=1e-10.

However, I've played around with it a bit now, and didn't immediately see a change. E.g., for eps=1e-12, I still get

x<0
      error    rel2max
8 -0.053904  -7.23e-09
A*x > p
       error    rel2max
31 -0.001641  -5.31e-09
36 -0.004082  -1.32e-08
66 -0.006870  -2.22e-08
72 -0.003382  -1.09e-08
84 -0.003510  -1.14e-08
Resolved 71.69%

while using a reasonable number for maxit. Although, for something like maxit = 200000, I at least get close(r) to it.

A*x > p
       error    rel2max
11 -0.000001  -3.97e-12
14 -0.000002  -5.60e-12
18 -0.000010  -3.21e-11
19 -0.000002  -5.06e-12
53 -0.000007  -2.29e-11
55 -0.000005  -1.64e-11
61 -0.000001  -3.69e-12
80 -0.000002  -6.98e-12

[bodono]

SCS will terminate if it either hits max_iters or the tolerance is below eps, so you might need to increase max_iters and decrease eps to get the accuracy you want. Though I would suggest trying a small tolerance, e.g., x>=1e-4 or something if you really need those constraints to hold (this is especially true for a first-order solver like SCS, but is somewhat true for any solver really).

[bodono]

Also setting verbose=True will allow you to see why SCS has decided to terminate and eye-ball the progress toward the end.