Describe the bug
Incorrect solution returned when using Gurobi solver
To Reproduce
import cvxpy as cp
import numpy as np
S = cp.Variable((2,2))
v = cp.Variable()
D = cp.Variable((2,2))
X = np.array([[1, 1], [1, 1]])
M = np.array([[1, 2/3], [2/3, 0.5]])
k = np.array([1, 2])
l = np.array([1, 0.5])
K = np.tile(k, (2,1))
L = np.tile(l, (2,1))
E = np.ones((2,2)) @ D @ np.diag(l)
constraints = [S >= 0, D >= 0, v >= 0]
constraints += [S <= np.ones((2,2)) @ S @ np.diag(l)]
constraints += [cp.sum(S) <= 1]
Q = S @ np.diag(k)
qs = 1 - cp.sum(cp.multiply(Q, L))
q0 = cp.multiply(qs, np.ones((2,1)))
dP = Q @ M.T
dp0 = q0 / 2
dP -= cp.multiply(cp.sum(dp0), L)
w = cp.sum(dp0)
G = dP @ np.diag(k)
constraints += [v >= G - D + E]
obj = v - cp.sum(cp.multiply(cp.multiply(S, K), X)) + w
lp_prob = cp.Problem(cp.Minimize(obj), constraints)
Expected behavior
Correct solution: -1.0
Outputs using other solvers:
In [32]: lp_prob.solve(solver=cp.ECOS)
Out[32]: -1.0000000000021085
In [33]: lp_prob.solve(solver=cp.OSQP)
Out[33]: -0.9999999999999996
In [34]: lp_prob.solve(solver=cp.MOSEK)
Out[34]: -1.0000000004011982
In [35]: lp_prob.solve(solver=cp.CVXOPT)
Out[35]: -0.9999999947383706
Output
Output using Gurobi solver
In [31]: lp_prob.solve(solver=cp.GUROBI)
Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored
Out[31]: -1.3333333333333335
In [36]: lp_prob.solve(solver=cp.GUROBI, verbose = True)
Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored
Parameter OutputFlag unchanged
Value: 1 Min: 0 Max: 1 Default: 1
Changed value of parameter QCPDual to 1
Prev: 0 Min: 0 Max: 1 Default: 0
Optimize a model with 18 rows, 9 columns and 41 nonzeros
Coefficient statistics:
Matrix range [5e-01, 3e+00]
Objective range [1e+00, 3e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 1e+00]
Presolve removed 14 rows and 5 columns
Presolve time: 0.00s
Presolved: 4 rows, 4 columns, 14 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 -3.0000000e+00 3.333333e-01 0.000000e+00 0s
1 -2.3333333e+00 0.000000e+00 0.000000e+00 0s
Solved in 1 iterations and 0.00 seconds
Optimal objective -2.333333333e+00
Out[36]: -1.3333333333333335
Version
Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored
I believe this is the issue. @rluce any Gurobi setting to prevent that?
@bstellato from what I can tell, there are no coefficients in the problem on the order of < 1e-13. In fact, when compiling that problem with ECOS, the smallest nonzero piece of problem data (in absolute value) is equal to 0.5. Specifically, the matrix term of ECOS problem data has
data['G'].data =
array([-1. , -1. , -1. , -1. , -1. ,
-1. , -1. , 1. , 2. , 1. ,
2.33333333, 1. , -1. , -1. , 1. ,
1. , 2. , 1. , 2.33333333, -1. ,
0.5 , -0.5 , 1. , 2.33333333, 1. ,
3. , 1. , -1. , -0.5 , 0.5 ,
1. , 1. , 2.33333333, 1. , 3. ,
-1. , 1. , -1. , 1. , -1. ,
-0.5 , 0.5 , -1. , 0.5 , -0.5 ])
with 45 entries, all nonzero. Compare to GUROBI's output,
Optimize a model with 18 rows, 9 columns and 41 nonzeros
Coefficient statistics:
Matrix range [5e-01, 3e+00]
Objective range [1e+00, 3e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 1e+00]
which states that the overall model has 41 nonzeros. Even assuming the number of nonzeros here refers to nnz of the constraint matrix, this is still missing 4 terms, relative to ECOS.
I would test more, but calling lp_prob.get_problem_data(solver='GUROBI') fails, since I don't have GUROBI installed on my machine.
Here is the output of the get_problem_data for Gurobi:
In [13]: lp_prob.get_problem_data(solver='GUROBI')[0]['F'].todense()
Out[13]:
matrix([[ 0. , -1. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , -1. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , -1. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
-1. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , -1. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , -1. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , -1. ],
[-1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , -1. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0.5 , -0.5 ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , -0.5 , 0.5 ,
0. , 0. , 0. , 0. ],
[ 0. , 1. , 1. , 1. , 1. ,
0. , 0. , 0. , 0. ],
[-1. , 2. , 1. , 2.33333333, 1. ,
0. , 1. , 0. , 0. ],
[-1. , 1. , 2. , 1. , 2.33333333,
1. , 0. , 0. , 0. ],
[-1. , 2.33333333, 1. , 3. , 1. ,
0. , 0. , -0.5 , 0.5 ],
[-1. , 1. , 2.33333333, 1. , 3. ,
0. , 0. , 0.5 , -0.5 ]])
In [14]: lp_prob.get_problem_data(solver='GUROBI')[0]['G']
Out[14]:
array([-0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
1., 1., 1., 1., 1.])
In [15]: lp_prob.get_problem_data(solver='GUROBI')[0]['q']
Out[15]: array([ 1., -2., -2., -3., -3., 0., 0., 0., 0.])
They are exactly the same as the output for OSQP:
In [19]: lp_prob.get_problem_data(solver='OSQP')[0]['F'].todense()
Out[19]:
matrix([[ 0. , -1. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , -1. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , -1. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
-1. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , -1. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , -1. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , -1. ],
[-1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , -1. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0.5 , -0.5 ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , -0.5 , 0.5 ,
0. , 0. , 0. , 0. ],
[ 0. , 1. , 1. , 1. , 1. ,
0. , 0. , 0. , 0. ],
[-1. , 2. , 1. , 2.33333333, 1. ,
0. , 1. , 0. , 0. ],
[-1. , 1. , 2. , 1. , 2.33333333,
1. , 0. , 0. , 0. ],
[-1. , 2.33333333, 1. , 3. , 1. ,
0. , 0. , -0.5 , 0.5 ],
[-1. , 1. , 2.33333333, 1. , 3. ,
0. , 0. , 0.5 , -0.5 ]])
In [20]: lp_prob.get_problem_data(solver='OSQP')[0]['G']
Out[20]:
array([-0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0., -0.,
1., 1., 1., 1., 1.])
In [21]: lp_prob.get_problem_data(solver='OSQP')[0]['q']
Out[21]: array([ 1., -2., -2., -3., -3., 0., 0., 0., 0.])
@rizalzaf Something isn't quite right -- which Gurobi version are you using?
@rluce I use Gurobi 8.1.1
I don't think it's Gurobi's problem, though.
I have a Julia implementation of the same code using Gurobi via JuMP, and it works fine
@rizalzaf I'm currently looking into what's happening on the gurobipy side. Could you actually try using the newer Gurobi version 9.0? While I could reproduce the problem with 8.1.1, the problem doesn't show with 9.0.
I see. I will let you know if I have a chance to update to v 9.0.
@rluce I just updated to Gurobi 9.0.0 and it solves the problem. Thanks.
Btw, coincidentally this issue is #900 :D
@rluce Why Gurobi keeps outputting the warning Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored even in version 9.0.0?
Most of the constraints in the problem are just saying that the variables need to be greater than 0.
Here is the complete picture. All coefficents of the linear constraint matrix being smaller than 1e-13 in absolute value are considered "noise", and automatically removed in order increase the sparsity of the matrix. The warning is issued in this case to alert about this, because most often such tiny coefficients (or exact zeros) arise from cancellations upon data assembly that should be cleaned up before passing over to the solver. In the example above, the matrix F has four structural zeros, so the warning is triggered no matter which Gurobi version is used.
In Gurobi version 8.1.1 there is a bug in gurobipy.Model._v811_addMConstrs, that may remove some nonzeros if the matrix being passed has structural zeros. This is the reason for the wrong answer this issue is originally about. Should there ever be another bugfix release v8.1.X, this bug will be fixed. A workaround in the interface would be:
diff --git a/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py b/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py
index 36ddd1a8..4cefad5f 100644
--- a/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py
+++ b/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py
@@ -128,6 +128,7 @@ class GUROBI(QpSolver):
elif hasattr(model, '_v811_addMConstrs'):
# We can pass all of A @ x == b at once, API only for Gurobi
# v811
+ A.eliminate_zeros() # Work around bug in gurobipy v811
sense = np.repeat(grb.GRB.EQUAL, A.shape[0])
model._v811_addMConstrs(A, sense, b)
else:
@@ -150,6 +151,7 @@ class GUROBI(QpSolver):
elif hasattr(model, '_v811_addMConstrs'):
# We can pass all of F @ x <= g at once, API only for Gurobi
# v811.
+ F.eliminate_zeros() # Work around bug in gurobipy v811
sense = np.repeat(grb.GRB.LESS_EQUAL, F.shape[0])
model._v811_addMConstrs(F, sense, g)
else:
Or, if this said warning is too annoying, one could also cleanup the linear constraint matrices unconditionally:
diff --git a/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py b/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py
index 36ddd1a8..4aa079d5 100644
--- a/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py
+++ b/cvxpy/reductions/solvers/qp_solvers/gurobi_qpif.py
@@ -92,8 +92,10 @@ class GUROBI(QpSolver):
P = data[s.P]
q = data[s.Q]
A = data[s.A].tocsr() # Convert A matrix to csr format
+ A.eliminate_zeros()
b = data[s.B]
F = data[s.F].tocsr() # Convert F matrix to csr format
+ F.eliminate_zeros()
g = data[s.G]
n = data['n_var']
Anyone with strong opinions on these two options?
@rluce Anyway to suppress the warning? Setting OutputFlag = 0 does not remove the warning.
The warning isn't suppressed because the Gurobi parameters are set only before the optimization is started, but after the Model object is populated. I can post a PR to address this.
Most helpful comment
@rluce I just updated to Gurobi 9.0.0 and it solves the problem. Thanks.
Btw, coincidentally this issue is #900 :D