GUROBI failing SOC dual-feasibility check

[rileyjmurray]

@SteveDiamond mentioned that GUROBI failed a dual variable complementarity check on his local machine. I've verified that now on my machine. The failing test is socp_1. Although the failure happens at a complementarity check, there is a larger issue that the dual value

dual_var = [np.array([0.28656026]), np.array([2.86560261, 2.22062537, 1.22062538])]

is very far from dual-feasible. (The dual-feasibility issue wasn't reported because we don't check for dual domains with ExpCone or SOC constraints.)

Paging contributors to the GUROBI interface who might be able to help: @rostyboost @rmcgibbo @rluce @gglockner.

Version

• OS: macOS 10.13.6
• CVXPY Version: 1.1.0a3 (master; slightly ahead of public 1.1.0a3 release).
• GUROBI Version: 9.0.1 (find via gurobipy.gurobi.version()).

[rluce]

Thanks for pinging, I'll check what's happening here.

[rileyjmurray]

Thanks @rluce!

I realize that CVXPY doesn't have a decent guide for maintaining solver interfaces. I just wrote a long comment about dual variables on another thread. You can find that comment here.

[rluce]

@rileyjmurray The post on the dual values was in intersting read, thanks. In this context I would like to understand what exactly SOC.dual_value actually represents. If I query it, I see that a nested array is returned, like for the socp_1 test,

>>> constraints[0].dual_value
[array([2.86560261]), array([[ 2.22062537],
       [ 1.22062538],
       [-1.33812376]])]

From what I see, these values are the dual values of some auxiliary linear constraints related to the SOC, but what is their interpretation?

Ultimately the reason for which the test socp_1 (and _2, _3) fail is related to these dual values: Upon construction of the gurobipy problem a different convention for ordering the constraints is used than later to retrieve the dual values.

To be specific, the first ordering is (1) linear equality constraints, (2) linear inequality constraints, and (3) SOC constraints expressed by (3a) one quadratic constraint and (3b) some auxiliary linear constraints. The (raw) vector of dual variables, solution["y"] in GUROBI.solve_via_data has this ordering, too, with an explicit dual value for the quadratic constraint.

Later in GUROBO.invert, it appears though that the dual value for the quadratic constraint (in the Gurobi representation) is not needed at all: The query

                leq_dual = utilities.get_dual_values(
                    solution['ineq_dual'],
                    utilities.extract_dual_value,
                    inverse_data[GUROBI.NEQ_CONSTR])

invokes the generic utilities.extract_dual_value function to extract the dual for the (cvxpy) SOC constraint, and only wants to access the dual values of the auxiliary constraints from (3b) above. Hence the dual value for the quadratic constraint (3a) is in the way and confuses the extraction logic.

The simplest fix appears to be chaning the order of the gurobipy constraints:

index 9ee2267a..519ada10 100644
--- a/cvxpy/reductions/solvers/conic_solvers/gurobi_conif.py
+++ b/cvxpy/reductions/solvers/conic_solvers/gurobi_conif.py
@@ -193,8 +193,7 @@ class GUROBI(SCS):
             variables += new_vars
             soc_start += constr_len
 
-        gur_constrs = eq_constrs + ineq_constrs + \
-            soc_constrs + new_leq_constrs
+        gur_constrs = eq_constrs + ineq_constrs + new_leq_constrs + soc_constrs
         model.update()
 
         # Set verbosity and other parameters

With this change in place, there are no failing tests anymore.

But since no tests are failing after this change, it means that there are no tests that actually ever use the dual value for the (gurobipy-) quadratic constraint. Do you know when/how this value is needed? If it is never needed, we should probably skip computing them, because this can be an expensive operation.

[rmcgibbo]

I believe that in our internal branch, we instruct Gurobi not to compute these duals and similarly ignore them in cvxpy, and nothing else in cvxpy seems to mind.

[rileyjmurray]

Thanks for that detailed investigation @rluce! It sounds like a good place to start is change the order of the gurobipy constraints as you suggest; this will get the SOCP tests passing. Because of how CVXPY passes the problem to gurobipy, it's possible that gurobipy QP dual variables are not needed (i.e. it's possible we're recovering all desired dual variables from constraints 3(b)). If that is indeed the case, then we can change CVXPY default settings to deactivate QP dual variable computation. I will look into this a little further over the next couple weeks.

For your question about exactly what SOC.dual_value represents, here's a walkthough:

• Given a t, a vector x, we create a conic representation of con = (cp.norm(x,2) <= t) by calling con = cp.SOC(t,x).
• CVXPY internally tracks primal expressions t and x by sticking them in a list called args (accessible by con.args). A list of numeric primal values can generically be recovered by prim_vals = [arg.value for arg in con.args].
• CVXPY stores dual variables in a way that matches the list structure of args: dual_vals = con.dual_values.
• Dual feasibility for this constraint can be checked by np.norm(dual_vals[1]) <= dual_vals[0].
• Complementarity can be computed by cs = prim_vals[0]*dual_vals[0] + prim_vals[1] @ dual_vals[1], and in exact arithmetic we would have cs == 0.0.

More generally, complementary slackness with a conic constraint (SOC, PSD, ExpCone) can be checked by

prim_vals = [arg.value for arg in con.args]
dual_vals = con.dual_values
comp_slack = cp.scalar_product(prim_vals, dual_vals).value
# comp_slack <= eps is necessary for an eps-accurate solution

We don't have a clean, generic function for checking conic dual-feasibility yet. I'll probably add such a function relatively soon.

[rluce]

The above change is part of the PR #955.

[rileyjmurray]

Resolved by @rluce's #955.