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.

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:

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, and send me the updated code (link to a gist is best).

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.

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))

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

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:

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.

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

I should have said 2 or 0, given the power of 2, but this is the context: Hessian matrix - Wikipedia.
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).

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.

The largest is 10000.

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

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:

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.

I’ve implemented 3D rotation matrices (see the implementation here) and bounded variables (as we discussed before), 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.

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

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).

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…