# How to solve a large quadratic programming problem

**URL:** <https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optimization, nonlinear\
**Created:** [March 5, 2024, 3:20pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200 "2024-03-05T15:20:31Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 6, 2024, 9:30am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/21 "2024-03-06T09:30:08Z")

</div>

> [@aris](#):
>
> I was looking for the most appropriate way to solve the following minimization problem in Julia

I’m curious as to the origin of the problem. Matrix `K` seems to be tall, i.e. few rows and many columns. So this indicates an underdetermined problem with much fewer equations than unknowns.

Is the problem really to find f from:

Kf = g

with the added condition that f\ge 0?

---

<div class="post-metadata">

**Author:** ![aris](https://avatars.discourse-cdn.com/v4/letter/a/258eb7/32.png) [@aris](https://discourse.julialang.org/u/aris)\
**Post date:** [March 6, 2024, 10:17am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/22 "2024-03-06T10:17:32Z")

</div>

@geoffroyleconte and @odow, your RipQP implementations work very well indeed. Thanks a ton for writing these down!

I suppose “solution” credit should go to @dpo for being the first to recommend using RipQP.

That’ll do for now, but I do plan to benchmark the other options recommended in this thread, whenever I manage to find some time to study how to use the different packages. My goal is to write a package similar to RegularizationTools.jl to automatically deal with these problems, so having multiple solvers would be beneficial.

> [@stevengj](#):
>
> I’m not sure why the condition number should matter here because there is a Tikhonov regularization (the \alphaα\alpha term), no?

Precisely, this is a textbook case of Tikhonov regularization.  
The condition number in the K example above has also been reduced **a lot** by SVD truncation during the preprocessing of the data.

> [@gdalle](#):
>
> you would be able to differentiate _through_ the optimization problem itself

This method sounds very interesting, thank you for the recommendation. I will try to go throught the documentation and see how that could work.

> [@BLI](#):
>
> I’m curious as to the origin of the problem.

The problem arises in the processing of NMR relaxation data. Your suggestion makes sense, it is effectively just Kf=g, but this problem is very ill posed, and slight perturbations in g (experimental noise in this case) would lead in dramatically different f results, if we just solve the problem using the aforementioned “naive” way.  
This is why we use Tikhonov regularization and introduce this \alpha term, to make the results more consistent and the resulting f distribution more “smooth”. Larger \alpha values result in smoother distributions which are less affected by the noise in g. Of course, if you increase \alpha too much you begin to lose a lot of information, so you need to strike a balance between ignoring noise and retaining actual information.

> [@aris](#):
>
> For anyone interested in more details for any reason, [here is a paper](https://doi.org/10.1016/j.pnmrs.2011.07.002)

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 6, 2024, 10:18am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/23 "2024-03-06T10:18:54Z")

</div>

> [@aris](#):
>
> I will try to go throught the documentation and see how that could work.

Here’s an example that is very similar to yours:

> **[Auto-tuning Hyperparameters · DiffOpt.jl](https://jump.dev/DiffOpt.jl/stable/examples/autotuning-ridge/)**
>
> Documentation for DiffOpt.jl.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 6, 2024, 10:31am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/24 "2024-03-06T10:31:59Z")

</div>

> [@aris](#):
>
> This is why we use Tikhonov regularization and introduce this \alphaα\alpha term, to make the results more consistent and the resulting fff distribution more “smooth”.

Is there then a chance that for the optimal solution, Kf=g is _not_ satisfied?

---

<div class="post-metadata">

**Author:** ![aris](https://avatars.discourse-cdn.com/v4/letter/a/258eb7/32.png) [@aris](https://discourse.julialang.org/u/aris)\
**Post date:** [March 6, 2024, 10:46am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/25 "2024-03-06T10:46:39Z")

</div>

I would imagine with data like these, equality would never be satisfied.  
You want to strike a balance between minimizing both terms ||Kf-g||^2 and ||\alpha f||^2. The L-curve method illustrates this goal nicely. If you’re interested in more information, [here is a link](https://www.youtube.com/playlist?list=PLMn2aW3wpAtMnYBVsPHF5n1Z0uMBvrVt8) to a YouTube playlist with some lectures on the subject, which I found very helpful.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 6, 2024, 10:56am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/26 "2024-03-06T10:56:32Z")

</div>

No its the other way around. Since f is underdetermined (more variables than lines), we are actually in high dimensionsional settings.

The space of f’s that satify exactlu Kf = g is large. Among them, you want to choose the one witht he smallest L2 norm.

Actually here whatever \alpha you should obtain the same result.

You shoudl actually implement the problem with Kf=g as _equality constraints_ and the loss being just `norm2(f)` : this form of the QP is equivalent in your case (K wide), but will probably be much faster. No need for \alpha.

---

<div class="post-metadata">

**Author:** ![aris](https://avatars.discourse-cdn.com/v4/letter/a/258eb7/32.png) [@aris](https://discourse.julialang.org/u/aris)\
**Post date:** [March 6, 2024, 11:15am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/27 "2024-03-06T11:15:34Z")

</div>

||Kf-g||^2 and ||f||^2 seem to be the things most people try to minimize.  
I don’t see Kf=g as a constraint being used (maybe I’m not looking hard enough).  
Perhaps there’s something I’m missing here, or I’m getting something wrong in the methodology (not unlikely).  
[https://www.sciencedirect.com/science/article/pii/S0377042700004143](https://www.sciencedirect.com/science/article/pii/S0377042700004143)

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [March 6, 2024, 11:25am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/28 "2024-03-06T11:25:33Z")

</div>

> [@aris](#):
>
> I don’t see Kf=gKf=gKf=g as a constraint being used (maybe I’m not looking hard enough).

Well, consider the sizes of things: size the matrix K is sized (n,d) with d \> n, the linear system of equations Kf = g will always have a solution (in fact, the space of solutions is d-n-dimensional.

So there is a way the optimizer will make the first part of the loss zero, namely \lVert Kf - g\rVert^{2} = 0. Among these solutions, the optimizer simply has to pick the one with smaller \lVert f \rVert, **regardless of the value of \alpha**.

Overwise said, Ridge regression does not handle high dimensional cases as well as Lasso regression would.

* * *

Edit: I forgot about the positivity constraint. Depending on the context, sometimes the positivity constraint is satisfied out of the box and you can just “check” for it at the end of the non-constrainted optimization. Otherwise, what I just said does not hold anymore and you have to use the full QP and find the right hyperparameter \alpha. But since the time is orders of magnitude smaller, I would still make this check for luck at the begining of my implementation 🙂

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 6, 2024, 3:04pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/29 "2024-03-06T15:04:01Z")

</div>

So – in principle one could just say:

```julia
f0 = K\g

```

which produces the `f` with minimum norm. If `f0 >= 0`, then the solution is found (perhaps poorly conditioned…).

What if `f0 >= 0` is not satisfied? Could one do:

```julia
M = nullspace(K)
f = f0 + M*z

```

and search for z such that `f >= 0`? This would reduce (slightly) the size of the problem?

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [March 6, 2024, 3:30pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/30 "2024-03-06T15:30:21Z")

</div>

It’s not because K is underdetermined that Kf = g is consistent. But if indeed it is, you should be able to remove the r variable from my formulation of the problem. You can also remove \alpha. RipQP will continue to apply. But if there’s a chance that Kf = g is inconsistent, or nearly inconsistent, the formulation with r is more robust.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [March 6, 2024, 4:05pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/31 "2024-03-06T16:05:27Z")

</div>

> [@lrnv](#):
>
> So there is a way the optimizer will make the first part of the loss zero, namely \lVert Kf - g\rVert^{2} = 0. Among these solutions, the optimizer simply has to pick the one with smaller \lVert f \rVert, **regardless of the value of \alpha**.

I’m wondering if I missed a turn in this discussion and am maybe misinterpreting. But are you talking about finding a minimum norm solution or finding the solution with Tikhonov regularization? In the former case \alpha is not involved. But with Tikhonov regularization, even neglecting nonnegativity constraints, you won’t generally get an exact solution to Kf = g. In fact as \alpha becomes large, you can expect the Tikhonov solution f to approach zero. Given that the problem is ill posed and K is ill conditioned, seeking an exact solution seems likely to give bad results, even if it is an exact solution of minimum norm.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 6, 2024, 4:15pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/32 "2024-03-06T16:15:53Z")

</div>

That’s what I meant when I asked if there is a chance that Kf=g is not satisfied. And I agree with you.

---

<div class="post-metadata">

**Author:** ![jondeuce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jondeuce/32/16378_2.png) [@jondeuce](https://discourse.julialang.org/u/jondeuce)\
**Post date:** [March 7, 2024, 7:37am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/33 "2024-03-07T07:37:29Z")

</div>

In case it is useful, I solve regularized small-scale nonnegative least-squares problems in my package [DECAES.jl](https://github.com/jondeuce/DECAES.jl/tree/171587a40423d9498d0210f56b336f1e2ab540f) using [this `NNLS` submodule](https://github.com/jondeuce/DECAES.jl/blob/171587a40423d9498d0210f56b336f1e2ab540f3/src/NNLS.jl#L151) which was forked from [NonNegLeastSquares.jl](https://github.com/jondeuce/NonNegLeastSquares.jl/blob/a122bf7acb498efcaf140b719133691e7c4cd03d/src/nnls.jl) some years ago.

To use this for your problem, you would pass `A = [K; sqrt(alpha) * I]` and `b = [g; zeros(n)]` where `size(K) = (m, n)`. Given that your `K` has many more columns than rows this is probably inefficient for your use case, though.

---

<div class="post-metadata">

**Author:** ![jondeuce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jondeuce/32/16378_2.png) [@jondeuce](https://discourse.julialang.org/u/jondeuce)\
**Post date:** [March 7, 2024, 7:59am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/34 "2024-03-07T07:59:13Z")

</div>

Also, slight aside, but I saw [an interesting twitter thread](https://x.com/docmilanfar/status/1763454859842552023?s=46&t=kxqKc2T4uClVZZghcawUAQ) about regularization when you have many more columns than rows. They [reference a paper](https://dl.acm.org/doi/abs/10.5555/3455716.3455885) showing that the optimal regularization parameter can sometimes be zero (or even negative!) due to implicit regularization in the least-norm solution `f = K' * ((K * K') \ g)`.

---

<div class="post-metadata">

**Author:** ![aris](https://avatars.discourse-cdn.com/v4/letter/a/258eb7/32.png) [@aris](https://discourse.julialang.org/u/aris)\
**Post date:** [March 7, 2024, 8:46am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/35 "2024-03-07T08:46:12Z")

</div>

Thanks for the suggestion @jondeuce.  
Interestingly enough, we are solving a very similar problem here.  
However, if I understand correctly you are interested in multiple T\_2 distributions (one for each point in the k-space I guess?) while I’m looking for T\_1-T\_2 correlation maps, which is the 2D analogue without any magnetic field gradients involved.

I installed DECAES.jl and tried `DECAES.nnls(A,b)` but it seems the size of b as defined in your comment will not work here, any insights into this?  
Also, would the nonneg\_lsq() solver from NonNegLeastSquares.jl work in a similar way here?

The paper you cite for the negative regularization parameter is very interesting, I’ll definitely have a read at some point.

---

<div class="post-metadata">

**Author:** ![jondeuce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jondeuce/32/16378_2.png) [@jondeuce](https://discourse.julialang.org/u/jondeuce)\
**Post date:** [March 7, 2024, 8:52am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/36 "2024-03-07T08:52:21Z")

</div>

Ahh cool, yeah probably DECAES as a whole is related but not totally relevant to you.

I edited my previous comment to make the sizes more clear; if `size(K) = (m, n)` then the padded arrays should be `size(A) = (m+n, n)` and `size(b) = (m+n,)`.

---

<div class="post-metadata">

**Author:** ![aris](https://avatars.discourse-cdn.com/v4/letter/a/258eb7/32.png) [@aris](https://discourse.julialang.org/u/aris)\
**Post date:** [March 7, 2024, 9:58am UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/37 "2024-03-07T09:58:36Z")

</div>

Works indeed! Although it is a lot slower than RipQP unfortunately.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 7, 2024, 1:14pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/38 "2024-03-07T13:14:08Z")

</div>

> [@jondeuce](#):
>
> To use this for your problem, you would pass `A = [K; sqrt(alpha) * I]`

The problem is that this is a huge matrix when `K` is “wide”, so you sacrifice a lot of efficiency by working in terms of A.

---

<div class="post-metadata">

**Author:** ![jondeuce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jondeuce/32/16378_2.png) [@jondeuce](https://discourse.julialang.org/u/jondeuce)\
**Post date:** [March 7, 2024, 6:58pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/39 "2024-03-07T18:58:50Z")

</div>

Yeah I figured this would probably be prohibitively costly. Also that code is single-threaded and doesn’t use BLAS, so would need to be tweaked to be competitive for large problems in any case.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [March 9, 2024, 6:10pm UTC](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200/40 "2024-03-09T18:10:58Z")

</div>

A sparse factorization will not care about the \sqrt{\alpha} I block, and it improves conditioning. It’s just that your vectors are longer, so a very small price to pay.

[Previous page](https://discourse.julialang.org/t/how-to-solve-a-large-quadratic-programming-problem/111200.md?page=1)