Add underdetermined Levenberg-Marquardt support (minimum-norm mode)

[ChrisRackauckas-Claude]

PR: Add underdetermined Levenberg-Marquardt support (minimum-norm mode)

Summary

This PR implements support for underdetermined nonlinear least-squares problems (n > m, more unknowns than equations) in the Levenberg-Marquardt algorithm, addressing issue #851.

Background

The current L-M implementation in DampedNewtonDescent assumes overdetermined systems (m ≥ n) and forms JᵀJ + λD. For underdetermined systems (n > m), this is inefficient because:

• JᵀJ is n×n (large)
• We want the minimum-norm solution among infinitely many solutions

Algorithm

For underdetermined systems, instead of solving (JᵀJ + λD)δ = -Jᵀf, we solve:

1. (JJᵀ + λI) z = -f(x) — This is an m×m system (smaller when n > m)
2. δ = Jᵀz — This gives the minimum-norm step

This approach:

• Keeps the linear system small (m×m instead of n×n)
• Finds the minimum-norm step that solves the linearized equations
• Is numerically well-conditioned with appropriate λ

Changes

lib/NonlinearSolveBase/src/descent/damped_newton.jl

• New mode: Added :minimum_norm to DampedNewtonDescentCache
• New kwarg: min_norm_mode in DampedNewtonDescent:
  • :auto (default): Auto-detect and use minimum-norm for underdetermined systems
  • :minimum_norm: Force minimum-norm mode
  • :disabled: Never use minimum-norm mode
• Auto-detection: Checks length(fu) < length(u) to detect underdetermined systems
• New caches: JJᵀ_cache, z_cache, J_raw for minimum_norm computations
• Helper functions: dampen_jacobian_minimum_norm!! and _extract_scalar_damping

lib/NonlinearSolveFirstOrder/test/underdetermined_tests.jl

New test file with comprehensive tests:

• Linear underdetermined systems
• Nonlinear underdetermined systems
• In-place formulation
• Overdetermined regression tests (ensure we didn't break anything)
• Large underdetermined systems (50 unknowns, 10 equations)

Usage

using NonlinearSolve

# Underdetermined: 2 equations, 4 unknowns
function f(u, p)
    return [
        u[1] + u[2] - 1.0,
        u[3] + u[4] - 2.0,
    ]
end

u0 = zeros(4)
prob = NonlinearLeastSquaresProblem(f, u0)
sol = solve(prob, LevenbergMarquardt())
# Automatically uses minimum-norm mode

References

• Issue: Levenberg–Marquardt variant for underdetermined least-square problems? #851
• Discussion on Discourse
• Chen et al. (2026), appendix B - Optica Express 34(4):5729
• Reference implementation: LeastSquaresOptim fork

Checklist

• Implementation in damped_newton.jl
• Tests for underdetermined systems
• Documentation updated in docstring
• CI passes (needs testing)

[stevengj]

Thanks Chris, I didn't realize that you had worked on this.

cc @aristeidis-karalis and @mochen4

[mochen4]

Thanks! This looks great!

[stevengj]

A simple test case would be to put in a linear function $f(x) = Ax - b$ for a "wide" matrix $A$ and check that the result of minimizing $\Vert f(x) \Vert_2^2$ starting at $x=0$ matches the minimum-norm solution A \ b.

[aristeidis-karalis]

@ChrisRackauckas-Claude, it's awesome you implemented this!

Just a comment:
I could not figure out (from the code) what D is in your original (m>n) formulation. Usually, it's best if it's something that scales with J², so that the result is independent from the scale of the problem. Some commonly used options are D = diag(JᵀJ) or D ~ ||J||²∙I (where ||J|| can be the Frobenius norm).

The same holds for the minimum-norm solution, it's often preferable if it is independent from the scale of f, J. Therefore, instead of only the λI option, it would be good to give the user the ability to choose λD also with, e.g., D = diag(JJᵀ) or D ~ ||J||²∙I. We used the latter, in particular D ~ ||J||²/m∙I

@stevengj

[stevengj]

Couldn't all three of these methods be replaced by

_extract_scalar_damping(D) = tr(D) / size(D,1)

assuming D is always a square matrix as seems to be the case?