Cvxpy: problem translated wrongly for different solvers?

[aliirmak]

I have this code snippet taken from here with minimal changes: https://cs.stanford.edu/~tachim/optimization_code.html
this is a svm classifier.

import cvxpy
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(5)

x1 = np.random.normal(2, 1, (2, 40))
x2 = np.random.normal(-2, 1, (2, 40))

w = cvxpy.Variable(2); b = cvxpy.Variable()
obj = 0
for i in range(40):
    obj += cvxpy.pos(1 - (w @ x1[:, i] + b))
    obj += cvxpy.pos(1 + (w @ x2[:, i] + b))
cvxpy.Problem(cvxpy.Minimize(obj), []).solve(solver=cvxpy.ECOS)

x = np.arange(-5, 5)
y = -(w.value[0] * x + b.value) / w.value[1]

plt.figure()
plt.plot(x, y, color='red')
plt.scatter(x1[0, :], x1[1, :], color='blue')
plt.scatter(x2[0, :], x2[1, :], color='green')

When obj=0.0*cvxpy.norm2(w0)**2 instead of obj=0, some solvers reach to a different solution. Some solvers give the same solution either way. Is that intentional?
version=1.0.25

[akshayka]

By different values, I assume you mean that the optimal objective is always close to zero, but not exactly zero. That's expected, because solvers return solutions by running numerical iterative algorithms, and some solvers are more accurate than others.

If my answer doesn't match what you mean, it would be helpful if you provided the code that we can run (without modification) that produces the exact output that you find confusing.

[aliirmak]

I understand the fact that solvers may have different precision. This is not a precision issue.

What I mean is, you use the same solver for the _almost_ same problem and you get a different result.
See below the two code pieces given below. They have the same solver.
The first one has obj=0, and the second one has obj=0.0*cvxpy.norm2(w0)**2 instead of obj=0.
To me, they are equivalent and _almost_ same problem. But they produce different (but close) results. CVXOPT generates the same result for both code pieces.

==1==

import cvxpy
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(5)

x1 = np.random.normal(2, 1, (2, 40))
x2 = np.random.normal(-2, 1, (2, 40))

w = cvxpy.Variable(2); b = cvxpy.Variable()
obj = 0
for i in range(40):
    obj += cvxpy.pos(1 - (w @ x1[:, i] + b))
    obj += cvxpy.pos(1 + (w @ x2[:, i] + b))
cvxpy.Problem(cvxpy.Minimize(obj), []).solve(solver=cvxpy.SCS)

x = np.arange(-5, 5)
y = -(w.value[0] * x + b.value) / w.value[1]

plt.figure()
plt.plot(x, y, color='red')
plt.scatter(x1[0, :], x1[1, :], color='blue')
plt.scatter(x2[0, :], x2[1, :], color='green')

==2==

import cvxpy
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(5)

x1 = np.random.normal(2, 1, (2, 40))
x2 = np.random.normal(-2, 1, (2, 40))

w = cvxpy.Variable(2); b = cvxpy.Variable()
obj = 0.0*cvxpy.norm2(w0)**2
for i in range(40):
    obj += cvxpy.pos(1 - (w @ x1[:, i] + b))
    obj += cvxpy.pos(1 + (w @ x2[:, i] + b))
cvxpy.Problem(cvxpy.Minimize(obj), []).solve(solver=cvxpy.SCS)

x = np.arange(-5, 5)
y = -(w.value[0] * x + b.value) / w.value[1]

plt.figure()
plt.plot(x, y, color='red')
plt.scatter(x1[0, :], x1[1, :], color='blue')
plt.scatter(x2[0, :], x2[1, :], color='green')

[akshayka]

Got it, thank you for clarifying.

CVXPY will translate the expression 0.0*cvxpy.norm2(w0)**2 by introducing a variable and a (trivial) conic constraint. Hence, when SCS is used, the two code snippets end up solving different, but equivalent, convex cone problems. Both problems have the same optimal value, but because they are specified differently, they give (just slightly) different answers, due to numerical issues. It's also totally possible for the two solves to give different optimal values for the variables.

In the future, one might imagine CVXPY could eliminate redundant expressions, like 0.0*cvxpy.norm2(w0)**2. These kinds of optimizations or pre-solves are currently out-of-scope for CVXPY.

[rileyjmurray]

[Once I finished writing this I saw Akshay's comment, but I'll post anyway.]

@aliirmak the difference is that in the first case, cvxpy compiles the problem into an LP, and in the second case, cvxpy compiles the problem into a second order cone program. If you use verbose=True, you'll see that

Lin-sys: sparse-indirect, nnz in A = 400, CG tol ~ 1/iter^(2.00)
eps = 1.00e-04, alpha = 1.50, max_iters = 2500, normalize = 1, scale = 1.00
Variables n = 83, constraints m = 160
Cones:  linear vars: 160

-- note the only cones are "linear vars". Compare to the second case, where you''ll find

Lin-sys: sparse-indirect, nnz in A = 406, CG tol ~ 1/iter^(2.00)
eps = 1.00e-04, alpha = 1.50, max_iters = 2500, normalize = 1, scale = 1.00
Variables n = 85, constraints m = 166
Cones:  linear vars: 160
    soc vars: 6, soc blks: 2

--- with two "soc blks", and total length of second-order-cone constraints equal to 6.

The number of nonzeros in the constraint matrix A is the same in both of these situations. SCS produces different results, because even trivially satisfied constraints (such as the second order cone constraints introduced in case 2) can affect the trajectory of an optimization algorithm.

It is unlikely that you would see this behavior with a commercial solver such as MOSEK or GUROBI, because commercial solvers usually include a presolve phase to eliminate trivial constraints. Open source projects (including rewriting systems such as cvxpy, or solvers like SCS) rarely have presolve functionality.

Interestingly, the design of cvxpy version "1.x" is such that users could contribute presolve methods, which cvxpy could integrate in a cross-solver way. Someone would have to step up and volunteer for such a thing, though.

[aliirmak]

Thanks for the excellent responses! I am closing the issue. Is this already asked for a feature?