# Solving a Convex Quadratic Program Using Gurobi & Mosek
**URL:** https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214
**Category:** Optimization (Mathematical)
**Tags:** optimization, mosek, gurobi, convex-optimization
**Created:** [August 28, 2021, 10:39am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214 "2021-08-28T10:39:02Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 28, 2021, 10:39am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/1 "2021-08-28T10:39:02Z")
I would like to solve the following Quadratic Program / Quadratic Programming:
\begin{aligned} \arg \min\_{x} \quad & \frac{1}{2} {x}^{T} P x + {x}^{T} q \\ \text{subject to} \quad & l \leq A x \leq u \end{aligned}
The matrix P is a Positive Definite matrix. It is a sparse matrix of size `40_000 x 40_000`.
I tried solving it using `Convex.jl` but it failed with `Out of Memory` (See [https://github.com/jump-dev/Convex.jl/issues/254](https://github.com/jump-dev/Convex.jl/issues/254)) at the problem formulation phase (I have `32 GB` so it might be memory is not properly handled).
I would like to solve it with `Gurobi.jl` / `Mosek.jl` using the direct interface.
Does anyone have a sample code for this?
Whta about formulation using `JuMP.jl`? Will it also have large overhead?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 28, 2021, 11:01am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/2 "2021-08-28T11:01:58Z")
JuMP should work fine for this.
The documentation has a range of examples
[https://jump.dev/JuMP.jl/stable/tutorials/Quadratic%20programs/portfolio/](https://jump.dev/JuMP.jl/stable/tutorials/Quadratic%20programs/portfolio/)
Or
[https://jump.dev/JuMP.jl/stable/manual/objective/#Set-a-quadratic-objective](https://jump.dev/JuMP.jl/stable/manual/objective/#Set-a-quadratic-objective)
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 28, 2021, 11:15am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/3 "2021-08-28T11:15:54Z")
I tried following the documentation and produced this:
```julia
using JuMP;
jumpModel = Model(Gurobi.Optimizer);
@variable(jumpModel, vT[1:numElements]);
@objective(jumpModel, Min, 0.5 * (vT' * mP * vT) + (vT' * vQ));
@constraint(jumpModel, vL .<= mA * vT);
@constraint(jumpModel, mA * vT .<= vU);
optimize!(jumpModel);
```
But I get errors: `Expression contains an invalid NaN constant. This could be produced by `Inf - Inf`.`.
I have `Inf` in the some of the elements of `vU` (Namely no constraint). It works with `Convex.jl` for small matrices.
This why I approached the forum.
Could it be that `JuMP` expect to `Inf` / `-Inf` in the constraints vectors (`vL` and `vU` above)?
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 28, 2021, 12:58pm UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/4 "2021-08-28T12:58:51Z")
The solution I had, for using `JuMP.jl` is discarding those volatile inequalities.
So the function I created, for reference, is:
```julia
using JuMP;
using Gurobi;
function SolveQpJump(mP, vQ, mA, vL, vU; hOptFun = nothing)
if (isnothing(hOptFun))
hOptFun = Gurobi.Optimizer;
end
jumpModel = Model(hOptFun);
@variable(jumpModel, vT[1:numElements]);
@objective(jumpModel, Min, 0.5 * (vT' * mP * vT) + (vT' * vQ));
vV = vL .≠ -Inf;
@constraint(jumpModel, vL[vV] .<= mA[vV, :] * vT);
vV .= vU .≠ Inf;
@constraint(jumpModel, mA[vV, :] * vT .<= vU[vV]);
optimize!(jumpModel);
return value.(vT);
end
```
I still wonder what’s the proper way to define `Inf` in `JuMP.jl`.
---
### Author:  [@mtanneau](https://discourse.julialang.org/u/mtanneau)
#### Post date: [August 28, 2021, 2:58pm UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/5 "2021-08-28T14:58:04Z")
When you give JuMP a constraint `A*x .>= vU`, my understanding is that it will try to reformulate it in a canonical form f(x) \in Set. If this involves moving the constant term to the f(x) term, and said constant is infinite… you’re out of luck.
For reference, I think [this issue](https://github.com/jump-dev/JuMP.jl/issues/2634) is related to the error you’re seeing.
I would recommend you create an additional variable s = Ax, and write your constraints as Ax = s, l \leq s \leq u (which, BTW, is more solver-friendly as it avoids duplicating the constraint matrix).
You can set lower/upper bound with [`set_lower_bound`](https://jump.dev/JuMP.jl/stable/reference/variables/#JuMP.set_lower_bound) and [`set_upper_bound`](https://jump.dev/JuMP.jl/stable/reference/variables/#JuMP.set_upper_bound), and these shouldn’t error if passed infinite values.
```julia
using JuMP;
numConstraints = size(mA, 1)
jumpModel = Model(Gurobi.Optimizer);
@variable(jumpModel, vT[1:numElements]);
@variable(jumpModel, s[1:numConstraints]); # <--- new
@objective(jumpModel, Min, 0.5 * (vT' * mP * vT) + (vT' * vQ));
# the following lines have been updated
@constraint(jumpModel, mA * vT .== s);
set_lower_bound.(s, vL)
set_upper_bound.(s, vU)
# Ready to optimize!
optimize!(jumpModel)
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 28, 2021, 9:38pm UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/6 "2021-08-28T21:38:11Z")
I would even go
```nohighlight
jumpModel = Model(Gurobi.Optimizer)
@variable(jumpModel, vL[i] <= s[i=1:numConstraints] <= vU[i])
```
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 29, 2021, 3:56am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/7 "2021-08-29T03:56:22Z")
> [@odow](#):
>
> ```Julia
> jumpModel = Model(Gurobi.Optimizer)
> @variable(jumpModel, vL[i] <= s[i=1:numConstraints] <= vU[i])
>
> ```
Does it support `Inf` in this formulation?
By the way, what’s the reason `Convex.jl` can’t formulate the problem with the same efficiency as `JuMP`?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 29, 2021, 6:12am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/8 "2021-08-29T06:12:13Z")
For Gurobi it does. Some solvers will have issues, which is the reason for the issue linked above.
JuMP and Convex.jl are very different, so it’s not surprising the have different performance characteristics.
---
### Author:  [@mtanneau](https://discourse.julialang.org/u/mtanneau)
#### Post date: [August 29, 2021, 7:43pm UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/9 "2021-08-29T19:43:23Z")
Convex.jl converts everything to conic form using the Disciplined Convex Programming (DCP) rules.
When you give it a quadratic objective, it will check for convexity then reformulate it as a second-order cone constraint. This implies computing a factorization of the objective matrix, which can be costly.
cc @ericphanson: I strongly suspect part of the reason is [this line](https://github.com/jump-dev/Convex.jl/blob/6371b5e3563c8998c7dc36afcd4597dd06cc14f0/src/atoms/second_order_cone/quadform.jl#L5), which will end up computing a dense cholesky factorization.
When you model the problem with JuMP, no such transformation occurs: the quadratic objective is passed as-is to Gurobi, who will be responsible for convexity checks and internal reformulations (if any).
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [August 29, 2021, 8:30pm UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/10 "2021-08-29T20:30:18Z")
Ah, could be. Convex has had memory problems long before that line but that does sound like it could be an issue. Perhaps we can provide a `PositiveSemidefinite` wrapper users can use to promise their matrix is PSD and skip the check.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 30, 2021, 4:40am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/11 "2021-08-30T04:40:48Z")
I am aware about DCP and the idea behind (Long time `CVX` user).
I thought that `JuMP` does some transformations as well in order to align all solvers.
By the way, would letting the user add some knowledge about the problem be helpful?
Namely, will there be any benefit telling `JuMP` the problem is Convex Quadratic Program? Can it pass it to the Optimizer?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 30, 2021, 5:37am UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/12 "2021-08-30T05:37:02Z")
> I thought that `JuMP` does some transformations as well in order to align all solvers.
Only if the constraint is not natively supported by the solver. In this case, Gurobi supports arbitrary quadratic objectives (including non-convex), so there is no need to analyze or reformulate it.
> Namely, will there be any benefit telling `JuMP` the problem is Convex Quadratic Program?
In this case, no. In other cases (e.g., conic solvers like SCS), you may be better to formulate this as minimizing a single epigraph variable with a second-order cone constraint.
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [September 14, 2021, 10:18pm UTC](https://discourse.julialang.org/t/solving-a-convex-quadratic-program-using-gurobi-mosek/67214/13 "2021-09-14T22:18:25Z")
> [@mtanneau](#):
>
> This implies computing a factorization of the objective matrix, which can be costly.
> cc @ericphanson: I strongly suspect part of the reason is [this line](https://github.com/jump-dev/Convex.jl/blob/6371b5e3563c8998c7dc36afcd4597dd06cc14f0/src/atoms/second_order_cone/quadform.jl#L5), which will end up computing a dense cholesky factorization
Just to update the thread, @mtanneau kindly updated Convex to use a sparse factorization when possible ([Use sparse symmetric factorization when possible by mtanneau · Pull Request #457 · jump-dev/Convex.jl · GitHub](https://github.com/jump-dev/Convex.jl/pull/457)) and I added an `assume_psd=false` argument to allow users to opt-out of the check altogether. Both are in Convex v0.14.15 which should be out shortly ([New version: Convex v0.14.15 by JuliaRegistrator · Pull Request #44875 · JuliaRegistries/General · GitHub](https://github.com/JuliaRegistries/General/pull/44875)).