# Solving a quadratic program inside a JuMP model solve

**URL:** <https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [May 3, 2018, 8:31pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690 "2018-05-03T20:31:32Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [May 3, 2018, 8:31pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/1 "2018-05-03T20:31:32Z")

</div>

I’m trying to solve a two-step optimization problem (to fit a complex econometric model) that requires solving a [fairly straight forward] quadratic program (as a function of the parameters of the model) at every iteration.

Is there any issue with initializing a new JuMP model to solve the quadratic program at each step? Would a different approach (perhaps a direct quad-prog solver or something) be preferable?

A simple example of what I want to do (for intuition) might be something like this:

```julia
# Quad Prog Function
function solveQadProg( H, f, ceq, beq)
    
    m = Model(solver=IpoptSolver(print_level=0, tol=1e-12))
    
    @variables m begin
        b[i=1:length(f)] >= 0
    end
    
    @objective(m, Min, sum(H[i]*b[i]^2 + f[i]*b[i] for i=1:length(f) ))
    
    @constraint(m, sum( b[i] * ceq[i] for i=1:length(ceq) ) == beq)
    
    status = JuMP.solve(m)
    
    b_min = getvalue(b)
    
    return b_min

end

### Main Model #user

M_main = Model(solver=IpoptSolver())

@variables M_main begin
    f_var[i=1:N]
    beq_var
end

## This would require solveQadProg to be registered properly, etc. and wouldn't work as is, but it gets the idea across
@expression(M_main, b_quadprog[1:N], @solveQadProgm, [H_data, f_var, ceq_data, beq]) 

@objective(M_main, sum(someData[i] - b_quadprog[i] for i=1:N)^2 / N)

status = JuMP,solve(M_main)
 

```

---

<div class="post-metadata">

**Author:** ![joehuchette](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joehuchette/32/32_2.png) [@joehuchette](https://discourse.julialang.org/u/joehuchette)\
**Post date:** [May 3, 2018, 8:47pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/2 "2018-05-03T20:47:27Z")

</div>

> Is there any issue with initializing a new JuMP model to solve the quadratic program at each step?

If the model generation time is a significant bottleneck in your code, then possibly yes. Otherwise, no. I’d suggest doing the easy thing (rebuilding the model from scratch each time) first, and then trying something fancier if it’s too slow.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [May 3, 2018, 10:04pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/3 "2018-05-03T22:04:58Z")

</div>

See also [https://github.com/JuliaOpt/JuMP.jl/issues/1241](https://github.com/JuliaOpt/JuMP.jl/issues/1241).

I’m working on a package currently called [SimpleQP](https://github.com/tkoolen/SimpleQP.jl), built on top of MathOptInterface.jl (the new backend for future versions of JuMP), which allows you to set up this kind of ‘parameterized QP’ and efficiently solve instances of it (with zero allocation). See [https://github.com/JuliaOpt/JuMP.jl/issues/1241#issuecomment-385806395](https://github.com/JuliaOpt/JuMP.jl/issues/1241#issuecomment-385806395) for a demo.

Note that the package is all of 11 days old, has zero documentation, has only been tested with the OSQP QP solver, and is still likely to change significantly.

---

<div class="post-metadata">

**Author:** ![miles.lubin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/miles.lubin/32/279_2.png) [@miles.lubin](https://discourse.julialang.org/u/miles.lubin)\
**Post date:** [May 4, 2018, 1:34am UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/4 "2018-05-04T01:34:40Z")

</div>

You could consider embedding the inner QP’s KKT optimality conditions directly as extra constraints and variables in the outer problem. This is a pretty standard approach for bilevel problems.

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [May 4, 2018, 1:49am UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/5 "2018-05-04T01:49:36Z")

</div>

Yeah, I thought of that and tried implementing it in a previous version of my model. It had a hard time finding feasible solutions, which discouraged me but I’ve made some simplifications to the model since then that might help. [update - it is having trouble finding feasible solutions… :X]

If this doesn’t work, do you think it makes any sense to use a QP solver in the form of a function, the solution to which is passed in to the main level in the form of an expression? I know JuMP isn’t really designed for this kind of thing, but I’d prefer to work in it rather going back to matlab, and I could provide a manual gradient to the QP…

Thanks for the thought!

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [May 7, 2018, 6:56pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/6 "2018-05-07T18:56:28Z")

</div>

Wondering if anyone has any thoughts on this being a good/bad approach? It seems like using the KKT optimality conditions as additional constraints makes it very difficult for the solver to find a feasible point…

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [May 7, 2018, 7:04pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/7 "2018-05-07T19:04:18Z")

</div>

> [@shoshievass](#):
>
> If this doesn’t work, do you think it makes any sense to use a QP solver in the form of a function, the solution to which is passed in to the main level in the form of an expression?

In my experience JuMP probably won’t be the limiting factor here: it’s more likely to be the startup time of your solver. I recommend taking a look at the wrapper code for your solver and figuring out if there is a way to keep a persistent one that you just reset every time (even this may not work well).

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [June 12, 2018, 4:15pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/8 "2018-06-12T16:15:37Z")

</div>

I’ve circled back to trying to get the KKT optimality conditions approach to work. I think the issue I’m hitting is this: my quadratic program requires that the solution be non-negative (e.g. b[t] \>= 0 for all t). When the boundary of 0 is not hit, the KKT model works perfectly. However, when the boundary is hit, the solver craps out and tells me there is no feasible solution.

@miles.lubin, @joehuchette, and @ExpandingMan - do you guys know what might be going on?  
My notebook is [here](https://github.com/shoshievass/QuadProgForAuctions/blob/master/QuadProg_KKT_Example/QuadProg_KKT.ipynb) and the data behind it is uploaded to [github as well](https://github.com/shoshievass/QuadProgForAuctions/tree/master/QuadProg_KKT_Example)

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [June 12, 2018, 4:47pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/9 "2018-06-12T16:47:03Z")

</div>

Your links are broken for me.

I only have the vaguest idea of what you are trying to do, but my first instinct would be to ask whether you can be sure that your Q matrix (i.e. the matrix of coefficients of the quadratic term in your objective function) is guaranteed to be positive definite. If it isn’t, your problem will not be convex. Most solvers should complain specifically that they got a non-positive-definite matrix rather than reporting the problem as infeasible in that case, so that’s not necessarily a great candidate for what’s going on. Can you fix your links?

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [June 12, 2018, 4:50pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/10 "2018-06-12T16:50:39Z")

</div>

Sorry, I just realized that github repo was private. Does [this](https://github.com/shoshievass/QuadProgForAuctions/blob/master/QuadProg_KKT_Example/QuadProg_KKT.ipynb) work?

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [June 12, 2018, 4:55pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/11 "2018-06-12T16:55:49Z")

</div>

I can prove that my Q is positive definite, actually. The math is [here](https://github.com/shoshievass/QuadProgForAuctions/blob/master/problem_doc.pdf)

---

<div class="post-metadata">

**Author:** ![ExpandingMan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/expandingman/32/866_2.png) [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)\
**Post date:** [June 12, 2018, 5:33pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/12 "2018-06-12T17:33:28Z")

</div>

> [@shoshievass](#):
>
> I can prove that my Q is positive definite, actually.

Ah, ok good. I’ve looked at this briefly, I’m afraid I can’t offer anything very helpful without going through all of this thoroughly, but if it were me doing this, I’d have far more confidence that the JuMP solution is correct than I would in myself having manually computed and entered the KKT conditions without any errors.

Can you try a different solver and see what error you get?

I’m sure you have been wrestling this for a while now, but since it isn’t that much code it might be worth entering in the KKT code from scratch just to be extra sure that it’s right. One of the big advantages of Julia and JuMP is that you could have done, e.g.

```julia
@objective(m, Min,
           -γ*q_b⋅(b .- α*c) - (γ/2)*sum(σ_sq .* (b .- α*c).^2)
)     

```

which might be helpful for making things a little bit clearer as you try to spot any possible errors (granted I have no idea how good the Jupyter notebook unicode input plugin is, if it exists at all).

Also, why are you allowing `b` and `omega` to go very slightly negative in the KKT formulation? It might not matter, but I wouldn’t do that.

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [June 12, 2018, 5:46pm UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/13 "2018-06-12T17:46:19Z")

</div>

Thanks for this! I’ve tried a bunch of solvers (gurobi among them) and the results are very similar. I’m very confident that the JuMP model is correct, and I’ve double checked the KKT conditions. What’s weird is that I keep getting infeasibility errors… I’m honestly not sure why that’s happening. Any insight would be greatly appreciated.

Re: the slightly negative variable constraints - this was my attempt to make it a bit easier for the solver to find something feasible but it didn’t work ☹

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 13, 2018, 1:12am UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/14 "2018-06-13T01:12:48Z")

</div>

> the results are very similar

Does Gurobi prove infeasibility or does it crash with numeric error?

Provided that there is a feasible solution, you are likely hitting numerical error. This comment is relevant:

> [@Ipopt infeasibility on travis-ci](https://discourse.julialang.org/t/ipopt-infeasibility-on-travis-ci/11612/2):
>
> This could be for any number of reasons. I don’t have any good way of troubleshooting this. You may want to try providing a starting solution using the start keyword in @variable. (See the docs [Nonlinear Modeling — JuMP -- Julia for Mathematical Optimization 0.18 documentation](http://www.juliaopt.org/JuMP.jl/0.18/nlp.html)) I would also try changing the tolerances in Ipopt. If you have any random inputs make sure you set the same seed. Your problem might also be poorly scaled. It looks like your objective value is on the order of…

---

<div class="post-metadata">

**Author:** ![Juser](https://avatars.discourse-cdn.com/v4/letter/j/34f0e0/32.png) [@Juser](https://discourse.julialang.org/u/Juser)\
**Post date:** [June 13, 2018, 3:06am UTC](https://discourse.julialang.org/t/solving-a-quadratic-program-inside-a-jump-model-solve/10690/15 "2018-06-13T03:06:37Z")

</div>

It could be something with numeric error or tolerances as @odow suggested. Alternatively, you can solve the interior optimization problem when doing the function eval for the econometric problem and then either compute the analytic gradient of the econometric objective function by hand (since the argmax of your interior function should be easy to differentiate because of the envelope theorem). If other parts of your function are better handled by Automatic Differentiation, you can do the differentiation manually for your argmax and overwrite/overload that function when the inputs are duals. See, for example, [this post](https://discourse.julialang.org/t/automatic-differentiation-slow-slower-than-finite-differences/11209/14).