# A new solver for robustified non-linear least squares problems, with non-Euclidean variables

**URL:** <https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085>\
**Category:** Optimization (Mathematical)\
**Created:** [March 14, 2023, 8:30pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085 "2023-03-14T20:30:40Z")\
**Posts on this page:** 7\
**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:** [March 14, 2023, 8:30pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/1 "2023-03-14T20:30:40Z")

</div>

I have written a new solver for robustified non-linear least squares problems, with non-Euclidean variables. It can compute auto-differentiated Jacobians when necessary. I have called the package [NLLSsolver](https://github.com/ojwoodford/NLLSsolver.jl). Please feel free to [clone it from GitHub](https://github.com/ojwoodford/NLLSsolver.jl) and try it out. Feedback is welcomed.

**Background:**  
Until recently I prototyped in MATLAB, and wrote production code in C++. I have been following Julia for a while, as a potential replacement for both. Recently I decided to try it out, and since I tend to do a lot of non-linear least squares optimizations of geometry problems in my work, my first goal was to optimize [bundle adjustment (BA) problems](https://grail.cs.washington.edu/projects/bal/).

I first had a look at existing Julia optimizer frameworks, [JuMP](https://github.com/jump-dev/JuMP.jl) and [NLPModels](https://github.com/JuliaSmoothOptimizers/NLPModels.jl). I found [an implementation of BA as an NLPModel](https://github.com/JuliaSmoothOptimizers/BundleAdjustmentModels.jl). However, I found that the model did not support robustification (necessary for downweighting outlier or erroneous measurements), did not seem able to model non-Euclidean variables such as 3D rotation matrices, and also had hand-coded derivatives (I’m not sure if auto-differentiation (AD) was or is now available). I needed to be able to do all 3 of these things. In addition, I was put off JuMP by the [assessment of the authors of that BA implementation](https://hal.science/hal-03155907/document), who made it sound even less suitable.

This led me to decide to implement my own optimization framework for robustified non-linear least squares problems. It’s been a very interesting first major project in Julia for me, and now I’m ready to share it publicly. It’s still in the early stages, and there isn’t much documentation (just the readme, and file comments), but there are a couple of examples, BA being one of them. I tested it against the NLPModels BA optimization on an unrobustified problem, and it was competitive (actually slightly faster, reaching the same cost), despite having auto-differentiated derivatives instead of hand-coded ones. There’s also still huge scope for reducing allocations.

Many thanks to everyone who has answered my questions on this forum since I started learning Julia. Hopefully I can now give a bit back to the community.

---

<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:** [April 14, 2023, 1:57pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/2 "2023-04-14T13:57:55Z")

</div>

You should now be able to install it from the package registry.

---

<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:** [April 14, 2023, 3:39pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/3 "2023-04-14T15:39:12Z")

</div>

this sounds very interesting!

I just saw rotations (SO(3) in your announcement, to you also have other non-Euclidean variables?

There is a bunch of algorithms for Optimization on Manifolds (like SO(3)) in `Manopt.jl` ( [Home · Manopt.jl](https://manoptjl.org/stable/), disclaimer – I am a developer). Maybe a solver close to yours might be [Levenberg–Marquardt · Manopt.jl](https://manoptjl.org/stable/solvers/LevenbergMarquardt/) ?

Concerning the derivative, you can employ Euclidean AD and “convert” the result to a Riemannian one. See for example Nicolas Boumals book [An introduction to optimization on smooth manifolds](https://www.nicolasboumal.net/book/index.html) section 4.7 for gradients (of embedded manifold) and I think 5.5 for the Hessian?  
The gradient conversion is implemented in [`Manifolds.jl`](https://juliamanifolds.github.io/Manifolds.jl/stable/) and [`ManifoldsDiff.jl`](https://juliamanifolds.github.io/ManifoldDiff.jl/stable/) (again disclaimer, I am among the developers). The Hessian is still future work.  
With that you can compute the Euclidean (Matrix-)Gradient and turn it into a Riemannian (Rotations-)Gradient.

---

<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:** [April 14, 2023, 7:52pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/4 "2023-04-14T19:52:12Z")

</div>

> [@kellertuer](#):
>
> I just saw rotations (SO(3) in your announcement, to you also have other non-Euclidean variables?

There are lower bounded, and upper & lower bounded variables, but I’m not sure how well they work in practice.

> [@kellertuer](#):
>
> There is a bunch of algorithms for Optimization on Manifolds (like SO(3)) in `Manopt.jl`

Do they support robustified least squares cost functions?

---

<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:** [April 14, 2023, 8:18pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/5 "2023-04-14T20:18:42Z")

</div>

> [@user664303](#):
>
> There are lower bounded, and upper & lower bounded variables, but I’m not sure how well they work in practice.

Just bounds would not make Euclidean variables non-Euclidean – that would just be constrained (but-still-Euclidean) Optimisation (like Augmented Lagrangian)

> [@user664303](#):
>
> Do they support robustified least squares cost functions?

Can you phrase this a bit more mathematical? What does robustified mean?

For example I know the robust PCA, which instead of the distance quared (Euclidean: norm-squared) minimised the distance (norm). But robust per se is a bit vague mathematically – so what are you referring to when you speak about robustified?  
Sure like for the PCA using the norm - is no longer least-squares – so I am unsure what your robust refers to.

In general all algorithms work with a cost (any cost `f(M,p)` where p is a point on the manifold `M`. Without a gradient for example `particle_swarm`. If you have a gradient, for example `gradient_descent` or `quasi_Newton`. If you even have a(n approximation of) Hessian, `trust_regions`.

---

<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:** [April 14, 2023, 9:51pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/6 "2023-04-14T21:51:55Z")

</div>

> [@kellertuer](#):
>
> Just bounds would not make Euclidean variables non-Euclidean – that would just be constrained (but-still-Euclidean) Optimisation (like Augmented Lagrangian)

It depends how the “bounds” are implemented. I’m using unconstrained, non-Euclidean variables, as opposed to constrained Euclidean variables. My solver doesn’t support constraints.

Robustified non-linear least squares is explained well in the [Ceres solver documentation](http://ceres-solver.org/nnls_tutorial.html).

---

<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:** [April 14, 2023, 10:30pm UTC](https://discourse.julialang.org/t/a-new-solver-for-robustified-non-linear-least-squares-problems-with-non-euclidean-variables/96085/7 "2023-04-14T22:30:19Z")

</div>

Thanks I will check the ceres docs. Ah I see – well with a Chain rule, we could handle that as well (to have the gradient), m but we do not have a specific solver for that. Might be interesting to give that a look.

And sure, the tradeoff/choice is usually either constraint Euclidean _or_ unconstraint Riemannian (non-Euclidean). `Manopt.jl` even mixes both and has two solvers (augmented Lagrangian and exact penalty method) that can do _constrained Riemannian optimisation_ (or in your words constrained non-Euclidean).