# Comparing non-linear least squares solvers

**URL:** <https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752>\
**Category:** Optimization (Mathematical)\
**Tags:** nonlinear, nonlinear-optimizati\
**Created:** [October 9, 2023, 1:37pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752 "2023-10-09T13:37:45Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 9, 2023, 1:37pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/1 "2023-10-09T13:37:45Z")

</div>

I’m creating this thread as a resource for people looking to solve non-linear least squares problems in Julia, so they know what tools are available, what features they have and how well they perform.

UPDATE: This comparison has been added to [juliapackagecomparisons.github.io](https://juliapackagecomparisons.github.io/pages/nonlinear_solvers/#nonlinear_least_squares_solvers)

I’ll start with solvers and features. This list is simply what I found and was able to test, and my understanding of the features.

| | [Ipopt](https://jump.dev/) | [JSO](https://github.com/JuliaSmoothOptimizers) | [NLLSsolver.jl](https://github.com/ojwoodford/NLLSsolver.jl) | [LeastSquaresOptim.jl](https://github.com/matthieugomez/LeastSquaresOptim.jl) |
| --- | --- | --- | --- | --- |
| Registered package(s) | ✅ | ✅ | ✅ | ✅ |
| Uses JuMP model definition | ✅ | ✅ | ❌ | ❌ |
| Bound constraints | ✅ | ✅ | ❌ | ✅ |
| Equality constraints | ✅ | ❌ | ❌ | ❌ |
| Non-linear constraints | ✅ | ❌ | ❌ | ❌ |
| Robustified cost functions | ❌ | ❌ | ✅ | ❌ |
| Non-Euclidean variables | ❌ | ❌ | ✅ | ❌ |
| Dense auto-differentiation | ✅ | ✅ | ✅ | ✅ |
| Supports sparsity | ✅ | ✅ | ✅ | ✅ |
| Sparse auto-differentiation | ✅ | ✅ | ✅ | ❌ |

The following performance comparisons are evaluated on an Apple M1 Pro CPU, using [this script](https://gist.github.com/ojwoodford/789e85197b18dcddb349e1f695bffc31).

Optimization times on small, dense problems of varying sizes:

 ![Screenshot 2023-10-12 at 15.28.44](https://global.discourse-cdn.com/julialang/original/3X/5/6/56ce0eb1f6ece60c9f8bf0a22797dc8e13c99545.png)

All the solvers used auto-differentiation to compute the derivatives. Ipopt got stuck in a local minimum for the final (largest) problem, but all other solvers found the global optimum.

Optimization times on medium sized, sparse problems of varying sizes:

 ![Screenshot 2023-10-12 at 15.28.26](https://global.discourse-cdn.com/julialang/original/3X/7/5/752376d5fd3f92bb84b591ca3483b554050712ad.png)

All the solvers except LeastSquaresOptim used auto-differentiation to compute the derivatives. The latter used hand-computed derivatives, since auto-differentiation of sparse systems is not supported. All solvers reached the global optimum.

If you believe other solvers should be in this comparison, please first add them to [the benchmark](https://gist.github.com/ojwoodford/789e85197b18dcddb349e1f695bffc31), and send me the updated code (link to a gist is best).

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 9, 2023, 1:48pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/2 "2023-10-09T13:48:16Z")

</div>

In terms of performance it would be nice to benchmark all solvers over a range of problem types and sizes.

For now I’ve benchmarked a few solvers on a multi-dimensional Rosenbrock problem, which leads to a sparse, bi-diagonal linear system. I’ve also added a similar, dense problem.

I’ve now added the graphs to the first post in this thread. I’ll aim to keep this graphs updated as more solvers are added.

---

<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:** [October 9, 2023, 2:21pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/3 "2023-10-09T14:21:02Z")

</div>

I guess the other option to add to the table is JuMP directly with Ipopt. Something like:

```julia
using JuMP, Ipopt
function main(f::Function, data::Vector, y::Vector)
    model = Model(Ipopt.Optimizer)
    @variable(model, x)
    @variable(model, residuals[1:length(y)])
    @constraint(model, residuals .== f(data, x) .- y)
    @objective(model, Min, sum(residuals.^2))
    optimize!(model)
    return value.(x)
end

```

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 9, 2023, 4:33pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/4 "2023-10-09T16:33:40Z")

</div>

> [@odow](#):
>
> ```julia
> @variable(model, residuals[1:length(y)])
> @constraint(model, residuals .== f(data, x) .- y)
> @objective(model, Min, sum(residuals.^2))
> 
> ```

I’m definitely interested to try Ipopt. Why do you add residuals as variables with a constraint, instead of just doing this?:  
`@objective(model, Min, sum((f(data, x) .- y) .^ 2))`

---

<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:** [October 9, 2023, 4:43pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/5 "2023-10-09T16:43:04Z")

</div>

See also the thread: [Suggestions needed for bound-constrained least squares solver](https://discourse.julialang.org/t/suggestions-needed-for-bound-constrained-least-squares-solver/35611)

The [StructuredOptimization.jl](https://github.com/kul-forbes/StructuredOptimization.jl) package is worth considering here too.

---

<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:** [October 9, 2023, 7:51pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/6 "2023-10-09T19:51:20Z")

</div>

So that the Hessian is diagonal. Its usually more efficient

---

<div class="post-metadata">

**Author:** ![mateuszbaran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mateuszbaran/32/221842_2.png) [@mateuszbaran](https://discourse.julialang.org/u/mateuszbaran)\
**Post date:** [October 9, 2023, 8:55pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/7 "2023-10-09T20:55:33Z")

</div>

Manopt.jl also has a non-linear least square solver: [Levenberg–Marquardt · Manopt.jl](https://manoptjl.org/stable/solvers/LevenbergMarquardt/) (though it’s not the main focus of the package). It supports manifold constraints only and can be used with AD and sparse Jacobians.

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 10, 2023, 8:19am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/8 "2023-10-10T08:19:32Z")

</div>

> [@odow](#):
>
> So that the Hessian is diagonal.

I don’t understand this. The hessian block for `x` will (in general) not be diagonal.

In any case I haven’t been able to successfully create a sum of squares cost function for Rosenbrock. Whatever I do generates a compilation error. Here’s one such implementation, along the lines of your approach:

```julia
model = Model(Ipopt.Optimizer)
x0 = [(n ≤ 16 && i % 2 == 0) ? 1.0 : -1.2 for i = 1:n]
@variable(model, x[i = 1:n], start = x0[i])
@NLexpression(model, FA[k = 1:(n - 1)], 10 * (x[k + 1] - x[k]^2))
@expression(model, FB[k = 1:(n - 1)], 1 - x[k])
@variable(model, residuals[1:2*(n-1)])
@constraint(model, residuals .== [FA; FB])
@objective(model, Min, sum(residuals .^ 2))

```

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [October 10, 2023, 10:29am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/9 "2023-10-10T10:29:52Z")

</div>

> [@user664303](#):
>
> I don’t understand this. The hessian block for `x` will (in general) not be diagonal.

The Hessian of the objective function with respect to all problem variables ([`x; residuals]`) will be diagonal with either 1 or 0s.

Changing

> [@user664303](#):
>
> `@constraint(model, residuals .== [FA; FB])`

to

```julia
@NLconstraint(model, [k=1:2*(n-1)], residuals[k] == [FA; FB][k])

```

should do it.

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 10, 2023, 10:56am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/10 "2023-10-10T10:56:36Z")

</div>

> [@jd-foster](#):
>
> The Hessian of the objective function with respect to all problem variables ([`x; residuals]`) will be diagonal with either 1 or 0s.

I’m very confused. A simple by-hand derivation of the Hessian of the cost function I define above won’t be either diagonal, or 1 or 0. Note that the variables `x` are not integers; they’re real values.

> [@jd-foster](#):
>
> ```julia
> @NLconstraint(model, [k=1:2*(n-1)], residuals[k] == [FA; FB][k])
> 
> ```

Thanks. This worked (ran). However, it’s incredibly slow/fails to converge. Is there a way to set the maximum optimization time in JuMP?

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [October 10, 2023, 11:04am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/11 "2023-10-10T11:04:42Z")

</div>

> [@user664303](#):
>
> either diagonal, or 1 or 0

I should have said 2 or 0, given the power of 2, but this is the context: [Hessian matrix - Wikipedia](https://en.wikipedia.org/wiki/Hessian_matrix#Use_in_optimization).  
If f(x\_1, r\_1, r\_2) = r\_1^2 + r\_2^2 then H(x\_1, r\_1, r\_2) = \left(\begin{array}{ccc}0 & 0 & 0 \\0 & 2 & 0 \\0 & 0 & 2\end{array}\right).

> [@user664303](#):
>
> However, it’s incredibly slow/fails to converge.

What is `n` in your tests?

> [@user664303](#):
>
> Is there a way to set the maximum optimization time in JuMP?

`JuMP.set_time_limit_sec(model::Model, limit::Real)`

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 10, 2023, 11:29am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/12 "2023-10-10T11:29:59Z")

</div>

> [@jd-foster](#):
>
> I should have said 2 or 0, given the power of 2, but this is the context: [Hessian matrix - Wikipedia](https://en.wikipedia.org/wiki/Hessian_matrix#Use_in_optimization).

Ah yes. Thanks. But then all the non-linearity is in the equality constraints. Do these not appear in some Hessian of the dual? My memory of Lagrangian multipliers is hazy, but I thought they too had a Hessian? If that information isn’t present in the system then it’s going to converge very slowly. No wonder it performs badly (see the graph below).

Is this sensible, @odow? What would it look like without the equality constraints; i.e. a standard, unconstrained least squares formulation? I can’t get anything written like that to compile.

> [@jd-foster](#):
>
> What is `n` in your tests?

The largest is 10000.

> [@jd-foster](#):
>
> JuMP.set\_time\_limit\_sec(model::Model, limit::Real)

Thanks. Here’s a plot with the time limited to 30s:

 ![Screenshot 2023-10-10 at 12.24.02](https://global.discourse-cdn.com/julialang/original/3X/c/4/c423940bb5d23ee1e1fc660bcc1c9cada8aa1e47.png)

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [October 10, 2023, 11:37am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/13 "2023-10-10T11:37:12Z")

</div>

> [@user664303](#):
>
> 10000 2000

Note it looks like you’re simply terminating at 2000 iterations? Based on line 23 in your comparison code:

```julia
function NLLStest(sizes)
    # First run the optimzation to compile everything
    options = NLLSOptions(maxiters = 2000)

```

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 10, 2023, 11:43am UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/14 "2023-10-10T11:43:54Z")

</div>

> [@jd-foster](#):
>
> Note it looks like you’re simply terminating at 2000 iterations?

Yes, NLLSsolver does indeed terminate after 2000 iterations currently. This happens only for the largest test case, (n=10000). The other solvers terminate after far fewer iterations (and at a higher cost), but after 30 seconds of runtime. Currently NLLSsolver can’t terminate based on runtime, but I will implement that. I’m also planning to add visualizations of the final cost.

In the meantime here’s a graph of tests for which all solvers reached the optimal solution:

 ![Screenshot 2023-10-10 at 13.14.38](https://global.discourse-cdn.com/julialang/original/3X/5/0/5073b472e26d1385ee98f61f407ce86585d60be2.png)

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [October 10, 2023, 1:25pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/15 "2023-10-10T13:25:52Z")

</div>

That is an interesting comparison, as @MateuszBaran mentioned, `Manopt.jl` can do least squares on arbitrary manifolds (those from `Manifolds.jl`), so I would be interested how (and which). non-Euclidean variables are realised in NLLSsolver.jl? Maybe we can join forces a bit in that non-Euclidean matter.

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 10, 2023, 1:34pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/16 "2023-10-10T13:34:51Z")

</div>

> [@kellertuer](#):
>
> I would be interested how (and which). non-Euclidean variables are realised in NLLSsolver.jl

I’ve implemented 3D rotation matrices ([see the implementation here](https://github.com/ojwoodford/NLLSsolver.jl/blob/ced2471c9225b176de4c800eca34fb903e1bb3f6/src/visualgeometry/geometry.jl#L56)) and bounded variables (as [we discussed before](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/4)), but it’s easy to add your own. It should be very straightforward to use variables from `Manifolds.jl` in `NLLSsolver.jl`. You just need to overload the `NLLSsolver.update` method for the particular variable type.

I’m happy to discuss it further. Please message me direct.

> [@kellertuer](#):
>
> `Manopt.jl` can do least squares on arbitrary manifolds (those from `Manifolds.jl`)

I’m slowly adding more solvers. I may try to add that one. But it’s a Euclidean problem.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [October 10, 2023, 1:38pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/17 "2023-10-10T13:38:07Z")

</div>

Yes, I remember. We could check whether you could be based on `ManifoldsBase.jl` (as well) in the update, then it might work on manifold more generally (but I would have to check how you update).

---

<div class="post-metadata">

**Author:** ![TheLateKronos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thelatekronos/32/12824_2.png) [@TheLateKronos](https://discourse.julialang.org/u/TheLateKronos)\
**Post date:** [October 10, 2023, 7:46pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/18 "2023-10-10T19:46:08Z")

</div>

This is a great comparison! Would you concider making a PR to add it to [JuliaPackageComparisons](https://github.com/JuliaPackageComparisons/JuliaPackageComparisons.github.io)?

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [October 10, 2023, 8:10pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/19 "2023-10-10T20:10:44Z")

</div>

Just one note on terminology: it is my understanding that Ipopt is the solver, while JuMP is the software that allows you to write your optimisation problem in nice mathematical, solver independent terms and from there it builds the actual problem to pass to the solver…

---

<div class="post-metadata">

**Author:** ![user664303](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/user664303/32/37843_2.png) [@user664303](https://discourse.julialang.org/u/user664303)\
**Post date:** [October 10, 2023, 10:02pm UTC](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752/20 "2023-10-10T22:02:29Z")

</div>

> [@TheLateKronos](#):
>
> This is a great comparison! Would you concider making a PR to add it to [JuliaPackageComparisons](https://github.com/JuliaPackageComparisons/JuliaPackageComparisons.github.io)?

Thanks. I’ll aim for that. I’d like to broaden the tests and generate a few more statistics first.

[Next page](https://discourse.julialang.org/t/comparing-non-linear-least-squares-solvers/104752.md?page=2)