Levenberg–Marquardt QR

This Levenberg–Marquardt variant solves the damped linear least-squares problem with column-pivoted QR. It avoids forming the Gram matrix JᵀJ, while retaining the Cholesky variant’s damping choices and convergence rules.

API reference

When to use it

Consider pivoted QR when forming normal equations loses too much accuracy, particularly for poorly conditioned dense Jacobians. Compare the Cholesky variant when sparse Jacobians are needed. QR changes the linear solve; it does not make a nonlinear fit globally optimal or remove sensitivity to parameter and residual scaling.

The method follows Madsen, Nielsen, and Tingleff (2004), §3.2, and MINPACK’s qrfac and qrsolv routines.

How it works

The solver factors the Jacobian with column pivoting and solves the stacked system [J; sqrt(μD)] h ≈ [-r; 0]. Diagonal regularization is applied through the QR solve rather than by constructing JᵀJ. The factorization can be reused while damping changes at the same iterate.

Scaled trust-region damping is the default, following the radius updates of Moré (1978) and MINPACK’s lmder with Basin’s bracketed parameter search. Nielsen damping is also available. Both variants use the same gain-ratio acceptance test and stopping settings as Cholesky LM.

Example

The same residual problem can run with either LM factorization. This example selects QR through with_pivoted_qr(); LevenbergMarquardtQr::new() is also available.

use basin::{DenseMatrix, Executor, Jacobian, LevenbergMarquardt, Residual};
use std::convert::Infallible;

struct Fit;

impl Residual for Fit {
    type Param = Vec<f64>;
    type Output = Vec<f64>;
    type Error = Infallible;

    fn residual(&self, x: &Self::Param) -> Result<Self::Output, Self::Error> {
        Ok(vec![x[0] - 1.0, 2.0 * (x[1] - 2.0)])
    }
}

impl Jacobian for Fit {
    type Jacobian = DenseMatrix;

    fn jacobian(&self, _: &Self::Param) -> Result<Self::Jacobian, Self::Error> {
        Ok(DenseMatrix::from_row_slice(2, 2, &[1.0, 0.0, 0.0, 2.0]))
    }
}

fn main() {
    let solver = LevenbergMarquardt::new()
        .with_pivoted_qr()
        .with_absolute_gradient_tolerance(1e-8);
    let result = Executor::from_start(Fit, solver, vec![0.0, 0.0])
        .max_iter(100)
        .run()
        .unwrap();

    assert!((result.param()[0] - 1.0).abs() < 1e-6);
    assert!((result.param()[1] - 2.0).abs() < 1e-6);
    assert!(result.cost() < 1e-12);
}

Practical considerations

Treat numerical rank as a solve safeguard

Rank is checked after regularization and column equilibration. The default relative threshold is epsilon(F) * (m + n) for m residuals and n parameters; with_relative_rank_tolerance changes it. This check concerns the regularized solve and does not establish that the original parameters are identifiable. There is no truncated solution or normal-equation fallback.

Nielsen damping increases after rank loss; trust-region damping uses rank loss to bound its parameter search. An unusable solve or exhausted retries can produce SolverFailed. For bounded problems needing minimum-norm steps under rank deficiency, compare trust-region reflective.

Keep damping and stopping choices explicit

The default absolute gradient threshold is 1e-8. The damping, optional progress tests, and numerical no-progress safeguard behave as described in the LM guide. A small damped step does not establish stationarity. Put iteration and evaluation budgets on the executor, and inspect both termination and fit quality.

Implementation

Levenberg–Marquardt QR publishes matching parameters and cost through PointState. It computes the gradient internally; the state does not expose it through GradientState. Residuals, Jacobians, and model caches belong to the solver. Fresh runs reset them, while exact solver-and-state checkpoints retain them.

For ordinary residuals, the reported cost is 0.5 * sum(r[i]^2). The solver requires Residual and Jacobian, without a CostFunction implementation. Wrapping the problem in RobustLeastSquares selects a robust loss and residual scale. The adapter supplies a safeguarded model, while costs and convergence use the robust objective and its gradient. Custom Jacobians additionally need ScaleRowsInPlace for robust losses.

Supported dense pairs are Vec<F>/DenseMatrix<F>, nalgebra::DVector<F>/DMatrix<F>, ndarray::Array1<F>/Array2<F>, and faer::Col<F>/Mat<F>, with F = f32 or f64. The default Vec path needs no BLAS/LAPACK provider.

Sparse Jacobians are unsupported: the sparse backends do not provide the required numerical column pivoting. The dense paths use pure Rust. Use LevenbergMarquardtQr::default() for f32.