qpsolvers / qpbenchmark

Benchmark for quadratic programming solvers available in Python
Apache License 2.0
110 stars 10 forks source link

Need to flip the sign of dual solution from Gurobi solver #112

Closed 563925743 closed 7 months ago

563925743 commented 7 months ago

Hi,

Thanks for the cool implementation!

The solution of dual variables output by Gurobi is negative. You need to flip its sign and then use it to compute the dual residual/ gap, etc.

I implemented the dual residual/ gap by myself. In my experiment, I can replicate the results from your codes by not flipping the sign, i.e. plug in the negative dual variables to compute the dual residual/ gap.

stephane-caron commented 7 months ago

Thank you for your report. Do you have a minimal example QP where the output of Gurobi is wrong?

We have unit tests to check dual multipliers with all solvers in https://github.com/qpsolvers/qpsolvers/blob/main/tests/test_solve_problem.py but the more the better. Feel also free to open a PR at https://github.com/qpsolvers/qpsolvers/pulls to propose your change.

563925743 commented 7 months ago

Thank you for your report. Do you have a minimal example QP where the output of Gurobi is wrong?

We have unit tests to check dual multipliers with all solvers in https://github.com/qpsolvers/qpsolvers/blob/main/tests/test_solve_problem.py but the more the better. Feel also free to open a PR at https://github.com/qpsolvers/qpsolvers/pulls to propose your change.

I found this issue when testing random QP only with inequality constraints:

import numpy as np
import gurobipy as gpy

n = 1000
m = 1000

P = np.random.rand(n, n)
q = (np.random.rand(n)).reshape(-1,1)
Q = P.T @ P + 1e-4*np.eye(n)

G = np.random.rand(m, n)
z0 = np.ones(n)
s0 = np.ones(m)
h = (G @ z0 + s0).reshape(-1,1)
#%%

m = gpy.Model()
m.Params.LogToConsole = 0

zhat = m.addMVar((n, 1), lb=-gpy.GRB.INFINITY, ub=gpy.GRB.INFINITY)
obj = 0.5 * zhat.transpose() @ Q @ zhat + q.T @ zhat
m.setObjective(obj)
m.addConstr(G @ zhat <= h)

m.optimize()

x_star = np.zeros(n)
for j in range(n):
    x_star[j] = zhat[j].x[0]

lamb = - np.array(m.Pi)
dual_res = np.linalg.norm((Q @ x_star + q.reshape(-1) + G.T @ lamb),np.inf)
dual_gap = np.abs(x_star.dot(Q @ x_star) + q.T @ x_star + h.T @ lamb)[0]

lamb_incorrect = np.array(m.Pi)
dual_res_incorrect = np.linalg.norm((Q @ x_star + q.reshape(-1) + G.T @ lamb_incorrect),np.inf)
dual_gap_incorrect = np.abs(x_star.dot(Q @ x_star) + q.T @ x_star + h.T @ lamb_incorrect)[0]

The dual variables of ineq constraints lamb output by Gurobi by m.Pi is actually negative. So it's sign needs to be flipped to positive before computing the dual residual and duality gap. If I don't flip the sign I will have dual_res_incorrect and dual_gap_incorrect to be the same as the output from your code.

stephane-caron commented 7 months ago

Thank you! :+1: The bug is now fixed with qpsolvers v4.3.1 and in the benchmark as of v2.2.1.