Large-scale least squares with ModelingToolkit.jl?

I’m starting to dig deep into ModelingToolkit.jl with the intention of using it to model a version of the SLAM problem. Ultimately, SLAM problems end up looking like large-scale nonlinear least squares problems. Instead of working with a single scalar objective function, large-scale least squares solvers exploit the extra structure in the problem by working with the residual function and its Jacobian directly.

ModelingToolkit.jl provides OptimizationSystem, but it is designed to wrap a single scalar objective function, and doesn’t provide a method for generate_jacobian.

NonlinearSystem seems like almost the right thing to use to build up the residual function by components (including having a generate_jacobian method) but requires the system be expressed as a vector of Equations, rather than simple expressions.

My current thought is to just make use of NonlinearSystem for now, prepending my residual expressions with 0 ~ to turn them into equations.

A few questions, though:

  1. Immediate term: will this hack work? I think so, but I haven’t tried yet.
  2. Slightly longer term, what is the right way to make things less hacky for this use case? A little generalization of NonlinearSystem? A new concrete sub-type of AbstractTimeIndependentSystem?
  3. Even longer term, where might large-scale least squares solvers for problems of this sort find their way into the ecosystem? Once I’ve got my residual function modeled, I’m perfectly comfortable cobbling together the solver bits I require from elsewhere in the julia ecosystem and beyond, but it would be nice to button things up with a more complete story fully contemplated by the reach of SciML.

I opened an issue on github about this, figuring it would ultimately lead to some opportunity to bring relevant extensions to ModelingToolkit.jl and/or packages in its larger ecosystem, but figured it would probably be best to cross post here, since at this point it is more about questions and discussions than specific software changes.

2 Likes

Well the gradient is the one-row Jacobian of a scalar equation.

Residual equations are the same thing, and sometimes can be a little bit more naturally written as an equation! Indeed, NonlinearSystem automatically moves them to the right and drops the zeros internally.

If there are solvers which specialize on this, then it would make sense to do. Right now there are no solvers in the SciML ecosystem to target. See https://scimlbase.sciml.ai/dev/ . So it really needs to happen the other way: first there needs to be a set of specific least squares solvers and a LeastSquaresProblem, and then a symbolic form would make sense.

1 Like