JuliaCon 2025

Solving nonlinear least squares problems lives in a very interesting part of the ecosystem. Naively, many would think that nonlinear least squares is very similar to nonlinear optimization, since after all it's solving a problem to minimize ||f(u)||_2 where ||.||_2 is the 2-norm (i.e. sum of squares). You can solve such a problem using a standard nonlinear optimization tool, like Optimization.jl or Optim.jl, but it turns out this is sub-optimal because it does not specialize on the special properties involved in the evaluation of the 2-norm.

In the Julia ecosystem, we have traditionally had the package LsqFit.jl that fits into this part of this ecosystem with a Levenberg–Marquardt method. However, its implementation, along with many other implementations (such as those from SciPy) are not numerically stable, and thus can have convergence issues for many problems.

This problem has been solved by integrating numerically stable and high-performance methods into NonlinearSolve.jl. The reason these methods are part of NonlinearSolve.jl and not an optimization package like Optimization.jl is because nonlinear least squares problems are actually more numerically similar to solving f(u)=0 than it is to doing an optimization. The reason for this is subtle. Newton's method is equivalent to the process of "linearize my function f at u, solve where the zero is to get a new u, and repeat". If you change that process to "linearize my function f at u, solve the least squares problem for the linear f to get a new u, and repeat", then you now have the generalization of Newton's method to nonlinear least squares problems, which is known as Gauss-Newton. The naive way of doing this is to simply change the Jacobian J in Newton to the square form J'J, that's all that's required! The downside of this form is that is the numerically unstable aspect of the aforementioned algorithms because this squares the condition number. Instead, if you replace the LU factorization with a QR factorization, the natural result that a QR factorization gives the solution to a linear least squares problem is precisely the solution. In other words, solving nonlinear least squares problems is the same as solving nonlinear systems of equations except you now allow for non-square Jacobians and you use QR factorizations instead of LU factorizations. Everything from Levenburg-Marquardt and beyond then simply falls out of this formulation, it's all numerically stable, and gives generalizations to line searches and trust regions to help ensure global convergence.

While that result is beautiful in its own right, it gives an odd aspect to the Julia ecosystem: the fastest and most robust nonlinear least squares solvers are in a package called NonlinearSolve.jl, which is not a package that mentions nonlinear optimization or nonlinear least squares in its name. So what do we do about discoverability, and how do we share this more widely with the community? The end of this is meant to be a discussion, since we really don't know at this point. Should we make a separate package NonlinearLeastSquares.jl that simply re-exports NonlinearSolve.jl and documents the nonlinear least squares tooling separately? Do we use a name like CurveFit that more matches the SciPy naming, but doing the same thing? Are there other ways to lead people towards this tooling? We don't know, but hopefully this talk will help people throughout the Julia ecosystem better understand how we got into the current situation to better brainstorm how to solve the last steps of the discoverability problem.