Optimal objective from a simple QP problem:
Answered% set up QP problem
QPproblem.Q = [8, 1; 1, 8];
QPproblem.c = [3, -4]';
QPproblem.A = [1, 1; 1, -1];
QPproblem.b = [5, 0]';
QPproblem.lb = [0, 0]';
QPproblem.ub = [inf, inf]';
QPproblem.x0 = [0, 1]'; % starting point
QPproblem.modelsense = 'min';
QPproblem.sense = ['<'; '='];
The optimal value of the objective from gurobi (version 9.01) via the gurobi.m matlab interface is:
-0.0138889
while the value constructed from c'*x + 0.5*x'*F*x:
-0.0208333
while the value constructed from c'*x + x'*F*x :
-0.0138889
It would appear as though c'*x + 0.5*x'*F*x is the correct objective because that is consistent with the dual optimality conditions. Please can you check the implementation?
-
Official comment
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
-
-
Hi Ronan,
Could you clarify your question? Gurobi optimizes with respect to the objective \(x^\top Q x + c^\top x\), with no scaling of the \( Q \) matrix.
A lot of theory for quadratic programming assumes that the objective function is of the form \( \frac{1}{2} x^\top D x + a^\top x \). The problem solved by Gurobi fits this form when we rewrite the objective as \(\frac{1}{2} x^\top Q' x + c^\top x\), where \( Q' = 2Q \).
Thanks,
Eli
0 -
-
Hi Eli,
most theory in quadratic optimisation assumes that the objective function is \( \frac{1}{2} x^\top Q x + c^\top x \) because the gradient thereof is \( Q x + c \) and it seemes like the latter is part of the optimality condition that is satisfied by the gurobi solutions, rather than \( \frac{1}{2} Q x + c \), so there seems to be some inconsistency. Please can you check the source code and let us know exactly what the optimality conditions are, that are being tested for in a gurobi QP problem?
Regards,
Ronan
0 -
-
-
Hi Ronan,
I think this might just be a matter of writing the problem in the correct way to apply the optimality conditions. We can write your problem as:
$$\begin{alignat}{2} \min_x\ && \frac{1}{2} x^\top &(2Q) x + c^\top x \\ \textrm{s.t.}\ && Ax &\leq b \\ && x &\geq 0. \end{alignat}$$
I'm going to replace your equality constraint with the single inequality constraint \( -x_1 + x_2 \leq 0\); this doesn't affect the optimal solution at all, and it's a little easier to concisely write the optimality conditions when the constraint senses are all the same.
Gurobi uses nonpositive dual values for \( \leq \) constraints in minimization problems (see here). So, the Lagrangian function is
$$ \mathcal{L}(x, \lambda) = \frac{1}{2} x^\top (2Q) x + c^\top x - \lambda^\top(Ax - b),$$
where \( \lambda \leq 0 \) are the Lagrange multipliers. Then
$$ \nabla_x \mathcal{L} (x, \lambda) = (2Q)x + c - A^\top \lambda. $$
The first-order optimality conditions for this problem (which are also sufficient, since \( Q \) is positive definite) are:
$$\begin{alignat}{2} \nabla \mathcal{L}_x (x^*, \lambda^*) &= 0 \qquad && \textrm{(stationarity)} \\ Ax^* \leq b \textrm{ and } x^* &\geq 0 && \textrm{(primal feasibility)} \\ \lambda^* &\leq 0 && \textrm{(dual feasibility)} \\ (\lambda^*)^\top(Ax^* - b) &= 0 &&\textrm{(complementarity)}. \end{alignat}$$
We can write and solve this problem in the MATLAB API as follows:
% Define data
c = [3; -4];
b = [5; 0];
Q = sparse([8, 1; 1, 8]);
A = sparse([1, 1; -1, 1]);
% Build model
qp.Q = Q;
qp.obj = c;
qp.A = A;
qp.rhs = b;
qp.lb = [0, 0];
qp.ub = [inf, inf];
qp.modelsense = 'min';
qp.sense = ['<'; '<'];
% Solve
result = gurobi(qp);It remains to check each of the optimality conditions for \( (x^*, \lambda^*) \), where \( x^* \) is the optimal primal solution and \( \lambda^* \) is the optimal Lagrange multipliers (dual values):
% Get primal/dual values
x = result.x;
lam = result.pi;
% Check optimality conditions
disp("Check 2*Q*x + c - A'*lam = 0 (stationarity):");
disp(2*Q*x + c - A'*lam);
disp('Check A*x - b <= 0 (primal feasibility):');
disp(A*x - b);
disp('Check x >= 0 (primal feasibility):');
disp(x);
disp('Check lam <= 0 (dual feasibility):');
disp(lam);
disp("Check lam'*(A*x - b) = 0 (complementarity):");
disp(lam'*(A*x - b))The output shows us that all of the optimality conditions are satisfied:
Check 2*Q*x + c - A'*lam = 0 (stationarity):
1.0e-09 *
0.2836
0.0774
Check A*x - b <= 0 (primal feasibility):
-4.9444
-0.0000
Check x >= 0 (primal feasibility):
0.0278
0.0278
Check lam <= 0 (dual feasibility):
-0.0000
-3.5000
Check lam'*(A*x - b) = 0 (complementarity):
1.2953e-13Does this help? Thanks!
Eli
0 -
-
Hi Eli,
that helped a lot. For two out of these 3 problems, gurobi is giving different solutions to cplex and pdco. Some of them are extreme cases where either c is zero or F is zero. What solutions and objectives do you get?
I need to be sure because I am writing the QP equivalent of this:
https://github.com/opencobra/cobratoolbox/blob/develop/src/base/solvers/solveCobraLP.m
Regards,
Ronan
QPproblem2.Q = [1, 0, 0; 0, 1, 0; 0, 0, 1];
QPproblem2.modelsense = 'min';
QPproblem2.c = -1*[0, 0, 0]';
QPproblem2.A = [1, -1, -1 ; 0, 0, 1];
QPproblem2.rhs = [0, 5]';
QPproblem2.lb = [0, -inf, 0]';
QPproblem2.ub = [inf, inf, inf]';
QPproblem2.csense = ['='; '='];
QPproblem3.Q = [1, 0, 0; 0, 1, 0; 0, 0, 1];
QPproblem3.modelsense = 'min';
QPproblem3.c = -1*[1, 1, 1]';
QPproblem3.A = [1, -1, 0 ; 0, 1, -1];
QPproblem3.rhs = [0, 0]';
QPproblem3.lb = [0, 0, 0]';
QPproblem3.ub = [inf, inf, inf]';
QPproblem3.sense = ['='; '='];
QPproblem4.c = [200; 400];
QPproblem4.A = [1 / 40, 1 / 60; 1 / 50, 1 / 50];
QPproblem4.b = [1; 1];
QPproblem4.lb = [0; 0];
QPproblem4.ub = [1; 1];
QPproblem4.modelsense = 'max';
QPproblem4.sense = ['L'; 'L'];
QPproblem4.Q = zeros(2,2);0 -
-
-
Hi Ronan,
In case these are snippets of code you're running, there are a few small errors:
- The model struct should have a field "obj" instead of a field "c"
- For QP2, the "csense" field should be named "sense"
- For QP4, 'L' is not a valid constraint sense (I suppose these should be \( \leq \) constraints, denoted with '<')
Just to make sure Gurobi reads your problem correctly, you could use
gurobi_write(QPproblemX, 'qp.lp')
to write an LP model file for the problem, then visually inspect it for errors.
Could you try solving the problems again after fixing these issues? This is what I get:
- For QP2, the optimal solution is \( x^* = (2.5, -2.5, 5) \), with optimal objective value \( z^* = 37.5 \)
- For QP3, the optimal solution is \( x^* = (0.5, 0.5, 0.5) \) with \( z^* = -0.75 \)
- For QP4, the optimal solution is \( x^* = (0, 0) \) with \( z^* = 0 \). This problem is strange, because \( x \in [0,1]^2 \), but the two constraints are \( \frac{1}{40} x_1 + \frac{1}{60} x_2 \leq 0 \) and \( \frac{1}{50} x_1 + \frac{1}{50} x_2 \leq 0 \). Thus, \( (0, 0) \) is the only feasible solution.
Thanks,
Eli
0 -
-
Hi Eli,
with the code below, I can recover your solutions. In the next post, I have a more complicated example, where there are Ax \ge b constraints.
Regards,
Ronan
%tests if Gurobi correctly solves a set of QP problems
solvers = prepareTest('requiredSolvers',{'gurobi'});
fprintf('Checking gurobi solver ...\n')clear params
params.OutputFlag = 0;tol= 1e-6;
c = [3; -4];
b = [5; 0];
Q = sparse([8, 1; 1, 8]);
A = sparse([1, 1; -1, 1]);% Build model
qp.Q = Q;
qp.obj = c;
qp.A = A;
qp.rhs = b;
qp.lb = [0, 0];
qp.ub = [inf, inf];
qp.modelsense = 'min';
qp.sense = ['<'; '<'];% Solve
result = gurobi(qp,params);% Get primal/dual values
x = result.x;
lam = result.pi;fprintf('%s\n','Problem qp')
% Check optimality conditions
disp('Check 2*Q*x + c - A''*lam = 0 (stationarity):');
disp(2*Q*x + c - A'*lam);
assert(norm(2*Q*x + c - A'*lam,inf)<1e-8)disp('Check A*x - b <= 0 (primal feasibility):');
disp(A*x - b);
assert(all((A*x - b)<=0))disp('Check x >= 0 (primal feasibility):');
disp(x);
assert(all(x>0))disp('Check lam <= 0 (dual feasibility):');
disp(lam);
assert(all(lam<=0))disp('Check lam''*(A*x - b) = 0 (complementarity):');
disp(lam'*(A*x - b))
assert(norm(lam'*(A*x - b),inf)<1e-8)QPproblem2.Q = sparse([1, 0, 0; 0, 1, 0; 0, 0, 1]);
QPproblem2.modelsense = 'min';
QPproblem2.obj = -1*[0, 0, 0]';
QPproblem2.A = sparse([1, -1, -1 ; 0, 0, 1]);
QPproblem2.rhs = [0, 5]';
QPproblem2.lb = [0, -inf, 0]';
QPproblem2.ub = [inf, inf, inf]';
QPproblem2.sense = ['='; '='];result = gurobi(QPproblem2,params);
fprintf('%s\n%s%g\n','QPproblem2','Optimal objective is: ', result.objval)
fprintf('%s\n','optimal primal is: ')
disp(result.x)
assert(abs(result.objval - 37.5)<tol)
assert(all((result.x - [2.5;-2.5;5])<tol))
fprintf('\n')QPproblem3.Q = sparse([1, 0, 0; 0, 1, 0; 0, 0, 1]);
QPproblem3.modelsense = 'min';
QPproblem3.obj = -1*[1, 1, 1]';
QPproblem3.A = sparse([1, -1, 0 ; 0, 1, -1]);
QPproblem3.rhs = [0, 0]';
QPproblem3.lb = [0, 0, 0]';
QPproblem3.ub = [inf, inf, inf]';
QPproblem3.sense = ['='; '='];result = gurobi(QPproblem3,params);
fprintf('%s\n%s%g\n','QPproblem3','Optimal objective is: ', result.objval)
fprintf('%s\n','optimal primal is: ')
disp(result.x)
assert(abs(result.objval + 0.75)<tol)
assert(all((result.x - [0.5;0.5;0.5])<tol))
fprintf('\n')QPproblem4.obj = [200; 400];
QPproblem4.A = sparse([1 / 40, 1 / 60; 1 / 50, 1 / 50]);
QPproblem4.b = [1; 1];
QPproblem4.lb = [0; 0];
QPproblem4.ub = [1; 1];
QPproblem4.modelsense = 'max';
QPproblem4.sense = ['<'; '<'];
QPproblem4.Q = sparse(2,2);result = gurobi(QPproblem4,params);
fprintf('%s\n%s%g\n','QPproblem4','Optimal objective is: ', result.objval)
fprintf('%s\n','optimal primal is: ')
disp(result.x)
assert(abs(result.objval + 0)<tol)
assert(all((result.x - [0;0])<tol))
fprintf('\n')fprintf('Done.\n')
0 -
-
Hi Eli,
this problem:
I solve like this:
resultgurobi = gurobi(gurobiQP,params);
Then test for stationarity like this:
res = (2*gurobiQP.Q*resultgurobi.x + gurobiQP.obj) - gurobiQP.A'*resultgurobi.pi - resultgurobi.rc;
disp(norm(res,inf))Does one have to change the sign of resultgurobi.pi depending on what gurobiQP.sense is?
Regards,
Ronan
0 -
-
-
Hi Ronan,
If Gurobi solves the QP using the barrier algorithm, there is no crossover phase and Gurobi will return an interior point solution. Could you try setting Method=1 to solve the problem using dual simplex, then check these conditions again? This way, you should end up with a boundary point solution.
Eli
0 -
-
Hi Eli,
yes, I will do that.
In the meantime, please can you express the full test of optimality conditions, in the general case, where reduced cost is shown and where there may be <, >, and = in the model.sense array?
Different solvers have different ways to pass out an optimal result, esp with respect to the sign of the dual variables. However, what I need is to be able to see that all solvers are returning the same result, where there is a unique optimum primal and dual vector to a given problem.
Regards,
Ronan
0 -
-
-
Hi Ronan,
The MATLAB code that you posted for the stationarity test looks correct to me (I believe it should work irrespective of the constraint senses). I don't know how other solvers handle the dual variable signs, though.
Thanks,
Eli
0 -
-
Hi Eli,
I am getting issues with optimality conditions apparently not being satisfied by possibly atypical QP problems.
https://groups.google.com/g/cobra-toolbox/c/Ikj7DE8hyzo
Please can you email ronan.mt.fleming@gmail.com so I can send you a problematic example?
Regards,
Ronan
0 -
-
-
I'll open a ticket for this in our online support portal.
0 -
Post is closed for comments.
Comments
13 comments