Levenberg–Marquardt variant for underdetermined least-square problems?

[stevengj]

What kind of problems is it mostly used for? Please describe.

Underdetermined nonlinear least-square problems $\min_x \Vert f(x) \Vert_2$ for $f: \mathbb{R}^n \mapsto \mathbb{R}^m$, with $n > m$ (more parameters than equations), trying to evolve the initial guess to a "nearby" solution, possibly with constraints (most simply, bound constraints on $x$).

Describe the algorithm you’d like

As described in this discourse post, it is a straightforward modification of Levenberg–Marquardt (L–M) for underdetermined cases.

Instead of solving a regularized/damped least-squares problem $(J^T J + \lambda I) \delta = -J^T f(x)$ for each step $x \to x + \delta$ as in ordinary L–M, we solve a regularized minimum-norm problem $(JJ^T + \lambda I) z = - f(x) \implies \delta = J^T z$ at each step (finding the smallest step length to solve the linearized equations as $\lambda \to 0^+$). This keeps the matrix $JJ^T + \lambda I$ small ($m \times m$), and it is well-conditioned if you choose e.g. $\lambda = \alpha \Vert J \Vert_F^2$ for a sufficient dimensionless damping parameter $\alpha > 0$ as recommended in the paper (we chose $\alpha = 0.1/n$ IIRC, but the basic idea here is similar to L–M so you can do something analogous to whatever damping you have now) — so you can directly use Cholesky to solve it. Or you could use LQ factorization of $J$ (= QR of $J^T$) or other methods if you want.

Other implementations to know about

Right now our code for the paper is in a LeastSquaresOptim.jl fork at Inverse-Design-of-Resonances/LeastSquaresOptim/levenberg_marquardt.jl at main · aristeidis-karalis/Inverse-Design-of-Resonances · GitHub

References

We described this in the reference Chen et al. (2026), appendix B, but I wouldn't be surprised if other authors have stumbled across the same trick. (At the time of publication, however, we couldn't find any examples.)

[stevengj]

cc @aristeidis-karalis and @mochen4

[oscardssmith]

I believe our Levenberg–Marquardt implementation uses a QR decomposition for better numerical stability. Would that work here? We generally try to avoid forming JJ'.

[stevengj]

First, note that standard L–M uses J'J ($J^T J$), not JJ' (not $JJ^T$), or more precisely it solves the Tikhonov-regularized least-square problem $\min_x (\Vert J\delta - f \Vert_2^2 + \lambda \Vert \delta \Vert_2^2)$ by one way or another. If the problem is underdetermined (so that the first term can be made exactly zero), this will find the minimum-norm solution, but then the magnitude of $\lambda$ becomes irrelevant so effectively there is no step-size regularization. Even then, how you solve this least-squares problem is crucial.

The problem is that if you use QR with regularization ($\lambda > 0$) then you end up doing QR on the $(m+n) \times n$ matrix [ J ; sqrt(λ) * I ], which is still $O(n^2)$ for $n > m$ (the underdetermined case) and much bigger than necessary if $n \gg m$. I think the solver you are using is in FastLevenbergMarquardt.jl here? It uses a fancier solver from MINPACK which only finds the dense QR factorization of $J$, and then does Givens rotations of the diagonal $\sqrt{\lambda} I_n$ part (or more generally some other diagonal matrix $D$) to take advantage of the sparsity, but the algorithm looks like it still has $O(n^2)$ cost, which is wasteful for large $n \gg m$.

It's also possible to avoid this with an SVD of $J$ instead of QR (since there is a simple modification of the pseudo-inverse to included the Tikhonov regularization), but as I said in my linked post this seems like overkill. (And for the ordinary over-determined case, I don't entirely see the point of not using Cholesky, since the $\lambda > 0$ regularization should normally make the problem well-conditioned? At least if you use $D = I$)

However, I want to make sure that they key point is not lost here:

Our suggestion (which worked rather well in our application) was not simply to change the solver in the underdetermined case, but to change the minimization problem: do a regularized minimum-norm step directly, instead of regularized least-squares. Just like the least-squares problem, you should be able solve the regularized minimum-norm problem in a variety of ways, either by direct Cholesky factorization of $JJ^T + \lambda I$ (assuming it is well conditioned) or by SVD or by QR-based algorithms. In particular, you would do QR of $J^T$ , or equivalently LQ of $J$, and even including the regularization it is $O(m^2 + mn)$ storage rather than $O(n^2)$. And probably there is a similar Givens trick to exploit the sparsity of $\lambda I$ … indeed, it's near the same problem just with a re-interpretation of the terms, so maybe you can re-use the same code.

… Though since the whole point of the $\lambda > 0$ regularization is to make it well-conditioned, I'm again not sure I quite see the point of avoiding the simpler and faster Cholesky solver.

[JamesWrigley]

I believe our Levenberg–Marquardt implementation uses a QR decomposition for better numerical stability. Would that work here? We generally try to avoid forming JJ'.

I think QR is the default, though it is possible to change the linear solver for LM: https://docs.sciml.ai/CurveFit/stable/tutorials/linsolve_choice/