Levenberg–Marquardt Algorithm

The Levenberg–Marquardt algorithm (LM) minimizes a sum of squared residuals by repeating one move: linearize the residuals, solve a damped version of the Gauss–Newton system for a step, try the step, and keep or reject it depending on how well the linear model predicted the decrease. The damping is raised after a poor prediction and lowered after a good one, so the method behaves like a short gradient-descent step far from the solution and like Gauss–Newton near it. Optimization calls Gauss–Newton and LM the standard solvers of calibration and bundle adjustment, and LM also refines poses, as in OpenCV’s default solver described in Perspective-n-Point. Nonlinear least squares derives the Gauss–Newton step, the damped normal equations and the gain ratio; this article treats LM as a procedure: what an iteration does, how the damping and the scaling matrix are chosen, what the stopping tests measure, and how three implementations (MINPACK, reached in Python through SciPy, and Ceres Solver) differ.

The method is named after Kenneth Levenberg, who proposed damping in 1944, and Donald Marquardt, who proposed it again in 1963, as Nonlinear Least Squares reports. Ceres Solver’s documentation calls it “the most popular algorithm for solving non-linear least squares problems” and says it was also the first trust-region algorithm to be developed, citing both papers. Both papers, and Jorge Moré’s 1978 paper on the algorithm’s implementation and theory, are cited here without having been read: the statements below about MINPACK come from SciPy’s C translation of its source.

Problem

Given $m$ residuals $\mathbf{r}(\boldsymbol{\theta}) \in \mathbb{R}^m$ of $n \le m$ parameters, find a local minimizer of

$$F(\boldsymbol{\theta}) = \tfrac{1}{2} \lVert \mathbf{r}(\boldsymbol{\theta}) \rVert^2 ,$$

using only the residuals and their Jacobian $J = \partial \mathbf{r} / \partial \boldsymbol{\theta}$. Gauss–Newton solves the same problem but can take steps that increase $F$, and fails when $J^\top J$ is singular; LM adds the safeguards that make the iteration descend from a wider range of starting points. It remains a local method: it finds a stationary point near its start, not the global minimum.

Inputs and Outputs

• Inputs: a residual function, its Jacobian (analytic, automatic, or by finite differences), a starting point, an initial damping value or trust-region radius, a scaling rule for the parameters, tolerances, and an evaluation budget.
• Outputs: the final parameters, the final cost, the reason for stopping (which test fired, or the budget ran out), and counts of residual and Jacobian evaluations. SciPy’s least_squares also returns the final Jacobian, from which a covariance $\sigma^2 (J^\top J)^{-1}$ of the estimate can be formed when the noise level $\sigma$ is known or estimated.

Intuition

Each iteration trusts a linear model of the residuals, but only so far. The damping $\lambda$ (or, equivalently, a trust-region radius $\Delta$) expresses how far: a large $\lambda$ gives a short step along the negative gradient, a small one gives the full Gauss–Newton step. After each trial the algorithm compares the decrease it got with the decrease the linear model promised. A ratio near one means the model was right and can be trusted further; a small or negative ratio means it was wrong, the step is shortened and, if the cost went up, thrown away.

The second ingredient is the scaling matrix $D$, which decides what “short” means. Parameters in vision come in very different units: an amplitude in thousands of grey levels, a position in pixels, a rotation in radians, a focal length in hundreds of pixels. Damping with the identity treats a change of one count in amplitude like a change of one pixel in position. Damping with a diagonal $D$ built from the Jacobian’s columns measures every parameter by how strongly it moves the residuals, which, by the derivation below, makes the step independent of the units the parameters happen to be written in.

Algorithm

The damped form, with an explicit damping parameter, as Nonlinear Least Squares presents it and as hand-written implementations often code it:

1. Choose a start $\boldsymbol{\theta}_0$, a damping $\lambda_0 > 0$, a scaling rule $D$, and thresholds $0 < \eta_0 < \eta_1 < \eta_2 < 1$.
2. Evaluate $\mathbf{r}$, $J$, the gradient $\mathbf{g} = J^\top \mathbf{r}$ and $D$ at the current point. Stop if a convergence test holds (see below).
3. Solve $(J^\top J + \lambda D^\top D)\,\boldsymbol{\delta} = -\mathbf{g}$.
4. Evaluate the cost at $\boldsymbol{\theta} + \boldsymbol{\delta}$ and the gain ratio $\rho$ (actual over predicted decrease).
5. Update the damping: decrease $\lambda$ if $\rho > \eta_2$, increase it if $\rho < \eta_1$.
6. If $\rho > \eta_0$, accept the step and go to 2 (a new Jacobian); otherwise keep $\boldsymbol{\theta}$ and go to 3 with the larger $\lambda$ (the same Jacobian).

MINPACK’s LMDER, the Fortran routine of Moré, Garbow and Hillstrom that SciPy’s least_squares(method='lm') calls in a C translation, is written as a trust-region method instead. In the SciPy v1.16.3 source each outer iteration does the following:

1. Evaluate the Jacobian and factor it once by QR with column pivoting.
2. On the first iteration, set the scaling $D$ to the column norms of the Jacobian (or to user-supplied values) and the radius to $\Delta = \text{factor} \cdot \lVert D \boldsymbol{\theta}_0 \rVert$, with factor 100 in SciPy; afterwards, keep each entry of $D$ as the running maximum of that column’s norm.
3. Test the gradient: stop if the largest cosine between the residual vector and a Jacobian column is at most gtol.
4. Inner loop: find the Levenberg–Marquardt parameter $\lambda$ (MINPACK’s par) whose step satisfies $\lVert D\boldsymbol{\delta} \rVert \approx \Delta$ (routine LMPAR, below), try the step, and compute the ratio of actual to predicted reduction, both relative to $\lVert \mathbf{r} \rVert^2$. On the first iteration, shrink $\Delta$ to $\lVert D\boldsymbol{\delta} \rVert$ if that is smaller.
5. If the ratio is at most 0.25, set $\Delta$ to a factor between 0.1 and 0.5 of $\min(\Delta, 10\lVert D\boldsymbol{\delta} \rVert)$, the factor chosen by a quadratic interpolation of the cost along the step, and raise $\lambda$ by its inverse; otherwise, if the ratio is at least 0.75 or $\lambda = 0$, set $\Delta = 2\lVert D\boldsymbol{\delta} \rVert$ and halve $\lambda$.
6. Accept the step if the ratio is at least $10^{-4}$; otherwise repeat the inner loop with the smaller $\Delta$, reusing the factorization.
7. Stop when both the actual and predicted relative reductions are at most ftol (with the ratio at most 2), when $\Delta \le \texttt{xtol} \cdot \lVert D\boldsymbol{\theta} \rVert$, or when the number of residual evaluations reaches the budget; other exit codes report tolerances too small for machine precision.

LMPAR returns $\lambda = 0$, the Gauss–Newton step, when that step already satisfies $\lVert D\boldsymbol{\delta} \rVert \le 1.1\Delta$; otherwise it searches for $\lambda > 0$ with $\lvert \lVert D\boldsymbol{\delta} \rVert - \Delta \rvert \le 0.1\Delta$, giving up after ten iterations with the best value found, and solves each damped system from the QR factor already computed, so trying a new $\lambda$ does not refactor the Jacobian.

Mathematical Formulation

The damped step and the trust region

The step solves the damped normal equations

$$\left(J^\top J + \lambda D^\top D\right) \boldsymbol{\delta} = -J^\top \mathbf{r},$$

which are the normal equations of the augmented linear least-squares problem $\min_{\boldsymbol{\delta}} \lVert J\boldsymbol{\delta} + \mathbf{r} \rVert^2 + \lambda \lVert D\boldsymbol{\delta} \rVert^2$. MINPACK solves that augmented problem, with $\sqrt{\lambda}\,D$ stacked under $J$ and zeros under $\mathbf{r}$, from its QR factor of $J$ rather than forming $J^\top J$, and Ceres’ documentation writes its subproblem the same way. With $D = I$ this is the system of Nonlinear Least Squares; with $D^\top D = \operatorname{diag}(J^\top J)$ it is the variant that article mentions in one sentence, which damps each parameter in proportion to its own curvature.

The trust-region problem is

$$\min_{\boldsymbol{\delta}} \ \tfrac{1}{2}\lVert J\boldsymbol{\delta} + \mathbf{r} \rVert^2 \quad \text{subject to} \quad \lVert D\boldsymbol{\delta} \rVert \le \Delta .$$

In our derivation, its Lagrangian conditions are the damped equations above with a multiplier $\lambda \ge 0$ and $\lambda(\lVert D\boldsymbol{\delta} \rVert - \Delta) = 0$: either the Gauss–Newton step lies inside the region and $\lambda = 0$, or the step lies on its boundary with some $\lambda > 0$. Every damped step is therefore a trust-region step for the radius $\lVert D\boldsymbol{\delta}_\lambda \rVert$, the equivalence Nonlinear Least Squares attributes to Madsen, Nielsen and Tingleff. LM thus belongs to the trust-region family that optimization contrasts with line searches. The two descriptions differ only in which quantity the algorithm updates: a damping-form implementation changes $\lambda$ directly; MINPACK changes $\Delta$ and solves for $\lambda$.

Gain ratio

With $\mathbf{g} = J^\top\mathbf{r}$, the linear model predicts the decrease

$$\text{pred} = F(\boldsymbol{\theta}) - \tfrac{1}{2}\lVert J\boldsymbol{\delta} + \mathbf{r} \rVert^2 = -\mathbf{g}^\top\boldsymbol{\delta} - \tfrac{1}{2}\boldsymbol{\delta}^\top J^\top J \boldsymbol{\delta}, \qquad \rho = \frac{F(\boldsymbol{\theta}) - F(\boldsymbol{\theta} + \boldsymbol{\delta})}{\text{pred}} .$$

For a damped step $\text{pred} = \tfrac{1}{2}\boldsymbol{\delta}^\top(J^\top J + 2\lambda D^\top D)\boldsymbol{\delta} > 0$ (substituting the damped equations, our derivation), so a negative ratio always means the cost rose. MINPACK writes both reductions relative to $\lVert\mathbf{r}\rVert^2$; Ceres’ documentation writes the ratio with differences of $\lVert F \rVert^2$, which is the same quantity.

Invariance under a rescaling of the parameters

Write the parameters in other units, $\boldsymbol{\theta}' = S\boldsymbol{\theta}$ with $S$ diagonal (kilo-counts instead of counts, millimetres instead of metres). Then $J' = JS^{-1}$, $J'^\top J' = S^{-1}J^\top J S^{-1}$ and, because $S$ is diagonal, $\operatorname{diag}(J'^\top J') = S^{-1}\operatorname{diag}(J^\top J)S^{-1}$. With $D^\top D = \operatorname{diag}(J^\top J)$ the damped system in the new units is

$$S^{-1}\left(J^\top J + \lambda \operatorname{diag}(J^\top J)\right)S^{-1}\boldsymbol{\delta}' = -S^{-1}J^\top \mathbf{r} \quad\Longrightarrow\quad \boldsymbol{\delta}' = S\boldsymbol{\delta},$$

the same step expressed in the new units (our derivation). The cost and the predicted decrease do not change, so neither do the gain ratio, the damping updates or any stopping test that compares relative quantities; in exact arithmetic the whole sequence of iterates maps by $S$. With $D = I$ the new system is $(J^\top J + \lambda S^2)(S^{-1}\boldsymbol{\delta}') = -J^\top\mathbf{r}$: LM in the new units is LM in the old units with $D^\top D = S^2$, a different algorithm. MINPACK’s default scaling, the running maximum of column norms, also transforms as $S^{-1}$, and its first radius $100\lVert D\boldsymbol{\theta}_0\rVert$ is unchanged, so by the same argument its iterates are invariant; with $D = I$ they are not, except at steps where $\lambda = 0$, since the undamped Gauss–Newton step is invariant under any linear reparameterization. The example below checks it on 100 starts (H1).

Mapping between conventions

Implementation	Step equation	Damping quantity	Scaling $D^\top D$
Nonlinear Least Squares	$(J^\top J + \lambda I)\boldsymbol{\delta} = -J^\top\mathbf{r}$	$\lambda$	$I$ (or $\operatorname{diag} J^\top J$)
MINPACK (par, radius $\Delta$)	$(J^\top J + \text{par}\, D^2)\boldsymbol{\delta} = -J^\top\mathbf{r}$, $\lVert D\boldsymbol{\delta}\rVert \approx \Delta$	$\lambda = \text{par}$, found from $\Delta$	running max of squared column norms, or $1/x_{\text{scale}}^2$ in SciPy
Ceres Solver (radius $\mu$)	$(J^\top J + \tfrac{1}{\mu}\operatorname{diag} J^\top J)\boldsymbol{\delta} = -J^\top\mathbf{r}$	$\lambda = 1/\mu$	squared column norms clamped to $[10^{-6}, 10^{32}]$

SciPy’s x_scale is the reciprocal of MINPACK’s diag: least_squares.py passes diag = 1 / x_scale, and x_scale='jac' passes no diag, which selects MINPACK’s column-norm scaling. Ceres’ levenberg_marquardt_strategy.cc computes the squared column norms of $J$, clamps them to [min_lm_diagonal, max_lm_diagonal], and passes $\sqrt{\operatorname{diag}(J^\top J)/\mu}$ to the linear solver as the diagonal $D$ of the damped problem; with the documentation’s augmented form, which stacks this $D$ under $J$, the normal equations carry $1/\mu$. The documentation, in our reading, is not consistent with this by a factor of two: it writes the subproblem as “$\frac{1}{2}\|J(x)\Delta x + F(x)\|^2 + \frac{1}{\mu} \|D(x)\Delta x\|^2$”, whose normal equations would carry $2/\mu$, but then appends “$\frac{1}{\sqrt{\mu}} D$” to the Jacobian, which gives $1/\mu$, as the code does.

Parameters

• Initial damping or radius. A hand-written LM needs $\lambda_0$. With $D^\top D = \operatorname{diag}(J^\top J)$, $\lambda$ is dimensionless and a value such as $10^{-3}$ means “start close to Gauss–Newton”; with $D = I$ it carries the units of $J^\top J$ and has to be tied to the problem, for example as a fraction of the largest diagonal entry. MINPACK takes the radius from factor (its documentation recommends values in $(0.1, 100)$, with 100 “generally recommended”), and SciPy fixes factor = 100. Ceres’ initial_trust_region_radius defaults to $10^4$, and its documentation states that the reciprocal, $10^{-4}$, is the initial regularization parameter.
• Update rule. MINPACK: shrink $\Delta$ by 0.1 to 0.5 when the ratio is at most 0.25, double it relative to the step when the ratio is at least 0.75. Ceres’ LM strategy (in levenberg_marquardt_strategy.cc), on an accepted step of quality $q$, divides the radius by $\max(1/3,\, 1 - (2q-1)^3)$, so that it grows by up to a factor of 3 after a good step and shrinks by up to a factor of 2 after a barely acceptable one (capped at max_trust_region_radius, default $10^{16}$) and resets a decrease factor $\nu$ to 2; on a rejected step it divides the radius by $\nu$ and doubles $\nu$, so repeated failures shrink the radius ever faster. The documentation’s generic trust-region template, which doubles or halves the radius, describes the family, not this strategy’s rule.
• Acceptance threshold. MINPACK accepts a step when the ratio is at least $10^{-4}$; Ceres’ min_relative_decrease defaults to $10^{-3}$.
• Scaling. x_scale='jac' (MINPACK’s column norms) or a vector of characteristic parameter sizes. SciPy 1.16.0 changed the default for method='lm' from x_scale=1 to 'jac'; 'trf' and 'dogbox' still default to 1. SciPy’s documentation says that changing the scale is equivalent to solving in variables $\boldsymbol{\theta}/x_{\text{scale}}$.
• Stopping tolerances. In least_squares with method='lm', ftol bounds the relative reduction of the cost, xtol the radius relative to $\lVert D\boldsymbol{\theta}\rVert$, and gtol the largest cosine between the residuals and a Jacobian column; all three default to $10^{-8}$ and must exceed machine epsilon. max_nfev defaults to $100n$ residual evaluations; SciPy 1.16.0 made this the 'lm' default also for a finite-difference Jacobian (previously $100n(n+1)$) and stopped counting the finite-difference evaluations in it.
• Restrictions in SciPy. method='lm' rejects bounds, robust losses other than 'linear', sparse Jacobians and problems with $m < n$; 'trf' and 'dogbox' handle those cases.

Complexity

Each outer iteration costs one Jacobian evaluation and one dense factorization: $O(mn^2)$ for MINPACK’s QR of the $m \times n$ Jacobian. Each further trial of $\lambda$ inside an iteration costs one residual evaluation plus LMPAR’s solves (at most ten) with the stored $n \times n$ triangular factor, each of which eliminates $\sqrt{\lambda}D$ by Givens rotations in $O(n^3)$ operations (our count from qrsolv), so a new $\lambda$ never repeats the $O(mn^2)$ factorization of $J$, and for $m \gg n$ the inner loop is cheap compared with a new Jacobian. An implementation that forms $J^\top J$ once and refactors $J^\top J + \lambda D^\top D$ for every $\lambda$ also pays $O(n^3)$ per trial (our count); Ceres’ documentation notes that its dog-leg strategy needs only one linear solve per successful step, whereas LM may need many. For large sparse problems such as bundle adjustment, the cost of the linear solve dominates and is reduced by the Schur complement, which bundle adjustment describes, or by an inexact step: Ceres’ documentation describes an inexact-step LM that solves the damped system approximately with an iterative solver and stops it with a forcing sequence, citing a convergence result of Wright and Holt.

Implementation

The example fits a 2-D Gaussian spot, the kind of fit used to locate stars, fluorescent markers or calibration dots to a fraction of a pixel, and asks how the scaling matrix behaves when the parameters have very different sizes. A $21 \times 21$ patch holds a spot of amplitude $A = 1000$ counts at $(x_0, y_0) = (10.3, 9.6)$ px with width $s = 1.6$ px on a background $B = 100$ counts, plus Gaussian noise with $\sigma = 10$ counts (one draw, seed 0). The amplitude and the width differ in size by a factor of 625. The same fit is then posed in a second parameterization, P′, with the amplitude in kilo-counts ($A' = A/1000$), and five solvers run from one grid of 100 starts (centre offsets of 0, ±2 and ±4 px in each axis and starting widths of 0.8, 1.6, 3.2 and 6.4 px; amplitude and background from the image’s maximum and median) in both parameterizations:

• LM-I and LM-diag, a hand-written LM with $D = I$ and with $D^\top D = \operatorname{diag}(J^\top J)$, using our own update rule (divide $\lambda$ by 3 when $\rho > 0.75$, double it when $\rho < 0.25$, accept when $\rho > 10^{-4}$) and stopping tests chosen to be invariant to a rescaling;
• MINPACK-jac and MINPACK-1, SciPy’s method='lm' with x_scale='jac' and x_scale=1.0;
• TRF-1, SciPy’s trust-region-reflective method with x_scale=1.0, its default.

Every point at which a solver evaluates the residuals is recorded, converted back to the units of P, and compared between the two parameterizations. Differences are measured in standard errors of the parameters at the reference minimum, $\sigma\sqrt{[(J^\top J)^{-1}]_{jj}}$, a unit that does not depend on how a parameter is written. Conventions: pixel centres at integer coordinates, $x$ the column and $y$ the row (downward); the model is sampled at pixel centres, not integrated over pixels; $s$ and $-s$ give the same spot, so comparisons use $\lvert s \rvert$. The code uses NumPy 2.3.5 and SciPy 1.16.3 and takes about 11 seconds on a shared machine.

import os
for var in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[var] = "1"                       # one BLAS thread, before NumPy loads
import numpy as np
from scipy.optimize import least_squares

np.seterr(all="ignore")                         # far starts can drive s to 0

# A 21 x 21 patch with a Gaussian spot; pixel centres at integers, x = column, y = row (down).
yy, xx = np.mgrid[0:21, 0:21].astype(float)
TRUE = np.array([1000.0, 10.3, 9.6, 1.6, 100.0])   # A (counts), x0, y0 (px), s (px), B (counts)
S = np.array([1e-3, 1, 1, 1, 1])                   # P' = S * P: amplitude in kilo-counts
NAMES = ["A", "x0", "y0", "s", "B"]

def model(t):
    A, x0, y0, s, B = t
    g = np.exp(-((xx - x0) ** 2 + (yy - y0) ** 2) / (2 * s ** 2))
    return B + A * g, g

def jacobian(t):                                   # d model / d (A, x0, y0, s, B), P units
    A, x0, y0, s, B = t
    g, dx, dy = model(t)[1], xx - x0, yy - y0
    cols = [g, A * g * dx / s**2, A * g * dy / s**2, A * g * (dx**2 + dy**2) / s**3, np.ones_like(g)]
    return np.column_stack([c.ravel() for c in cols])

clean = model(TRUE)[0]
noise = np.random.default_rng(0).normal(0.0, 10.0, clean.shape)
noisy = clean + noise

class Problem:
    """Residuals and Jacobian in variables u = scale * theta; records every evaluation in P units."""
    def __init__(self, image, scale):
        self.image, self.scale, self.points = image.ravel(), scale, []
    def res(self, u):
        t = u / self.scale
        self.points.append(t.copy())
        return model(t)[0].ravel() - self.image
    def jac(self, u):
        return jacobian(u / self.scale) / self.scale      # J' = J S^-1

def lm(prob, u0, damping, max_nfev=1000):
    """Hand-written LM: (J^T J + lam D^2) d = -J^T r with D = I or D^2 = diag(J^T J)."""
    u = np.array(u0, float)
    r = prob.res(u); F = 0.5 * r @ r; nfev, it = 1, 0
    J = prob.jac(u); H = J.T @ J
    lam = 1e-3 if damping == "diag" else 1e-3 * H.diagonal().max()
    lam_cap = 1e20 if damping == "diag" else 1e20 * H.diagonal().max()
    while True:
        g, rn = J.T @ r, np.linalg.norm(r)
        cosine = np.abs(g) / np.maximum(np.linalg.norm(J, axis=0) * rn, 1e-300)
        if rn == 0 or cosine.max() <= 1e-12:              # gradient test
            return u, F, nfev, it, True
        D2 = H.diagonal() if damping == "diag" else np.ones(len(u))
        while True:                                       # inner loop: retry with more damping
            if nfev >= max_nfev or lam > lam_cap:
                return u, F, nfev, it, False
            try:
                d = np.linalg.solve(H + lam * np.diag(D2), -g)
            except np.linalg.LinAlgError:
                return u, F, nfev, it, False
            rt = prob.res(u + d); Ft = 0.5 * rt @ rt; nfev += 1
            pred = -(g @ d + 0.5 * d @ H @ d)             # decrease predicted by the linear model
            rho = (F - Ft) / pred if pred > 0 else -1.0   # gain ratio
            if rho > 0.75:
                lam /= 3
            elif rho < 0.25:
                lam *= 2
            if rho > 1e-4:
                break                                     # accept the step
        dec = (F - Ft) / F if F > 0 else 0.0
        u, r, F, it = u + d, rt, Ft, it + 1
        if dec < 1e-12 or np.max(np.abs(d) / np.maximum(np.abs(u), 1e-12)) < 1e-10:
            return u, F, nfev, it, True                   # cost or step test
        J = prob.jac(u); H = J.T @ J

def scipy_fit(prob, u0, method, x_scale):
    try:
        s = least_squares(prob.res, u0, jac=prob.jac, method=method, x_scale=x_scale,
                          ftol=1e-8, xtol=1e-8, gtol=1e-8, max_nfev=1000)
    except (ValueError, np.linalg.LinAlgError):
        return np.full(5, np.nan), np.inf, len(prob.points), 0, False
    return s.x, s.cost, s.nfev, s.njev, s.status > 0

METHODS = {
    "LM-I": lambda p, u0: lm(p, u0, "identity"),
    "LM-diag": lambda p, u0: lm(p, u0, "diag"),
    "MINPACK-jac": lambda p, u0: scipy_fit(p, u0, "lm", "jac"),
    "MINPACK-1": lambda p, u0: scipy_fit(p, u0, "lm", 1.0),
    "TRF-1": lambda p, u0: scipy_fit(p, u0, "trf", 1.0),
}
PARAMS = {"P": np.ones(5), "P'": S}

def fold(t):                                       # s and -s give the same spot
    return np.r_[t[:3], abs(t[3]), t[4]]

def cost(t, image):
    return 0.5 * np.sum((model(t)[0] - image) ** 2)

def moment_start(image, dx, dy, s0):
    return np.array([image.max() - np.median(image), 10.3 + dx, 9.6 + dy, s0, np.median(image)])

# ---- data checks ----
print(f"D1 noise SD = {noise.std():.3f} counts (bar |SD - 10| < 1.0): "
      f"{'held' if abs(noise.std() - 10) < 1.0 else 'FAILED'}")
ratio = np.abs(TRUE).max() / np.abs(TRUE).min()
cn = np.linalg.norm(jacobian(TRUE), axis=0)
print(f"D2 parameter magnitude ratio in P = {ratio:.1f} (bar 300-3000): "
      f"{'held' if 300 <= ratio <= 3000 else 'FAILED'}; Jacobian column-norm ratio at truth: "
      f"P {cn.max() / cn.min():.1f}, P' {(cn / S).max() / (cn / S).min():.1f}")

# ---- exact-data and Jacobian checks ----
worst = 0.0
for name, fit in METHODS.items():
    for key, sc in PARAMS.items():
        prob = Problem(clean, sc)
        x, F = fit(prob, moment_start(clean, 1, -1, 2.4) * sc)[:2]
        err = np.max(np.abs(fold(x / sc) - TRUE) / np.abs(TRUE))
        worst = max(worst, err)
        print(f"C1 noise-free {name:11s} {key:2s}: max rel error {err:.1e}, cost {F:.1e}: "
              f"{'held' if err < 1e-6 and F < 1e-6 else 'FAILED'}")
t1 = moment_start(noisy, 1, -1, 2.4)
Jfd = np.column_stack([(model(t1 + h)[0] - model(t1 - h)[0]).ravel() / (2 * h.sum())
                       for h in np.diag(1e-6 * np.maximum(np.abs(t1), 1))])
c2 = np.abs(Jfd - jacobian(t1)).max() / np.abs(jacobian(t1)).max()
print(f"C2 Jacobian against central differences: {c2:.1e} (bar 1e-6): {'held' if c2 < 1e-6 else 'FAILED'}")
a, b = Problem(noisy, PARAMS["P"]), Problem(noisy, PARAMS["P"])
sa = least_squares(a.res, t1, jac=a.jac, method="lm", ftol=1e-8, xtol=1e-8, gtol=1e-8, max_nfev=1000)
sb = least_squares(b.res, t1, jac=b.jac, method="lm", x_scale="jac", ftol=1e-8, xtol=1e-8,
                   gtol=1e-8, max_nfev=1000)
c4 = np.abs(sa.x - sb.x).max()
print(f"C4 'lm' default x_scale against 'jac': max |dx| = {c4:.1e}, nfev {sa.nfev} and {sb.nfev}: "
      f"{'held' if c4 == 0.0 and sa.nfev == sb.nfev else 'FAILED'}")
for label, kw in (("loss='huber'", {"loss": "huber"}), ("bounds", {"bounds": (0, np.inf)})):
    try:
        least_squares(a.res, t1, method="lm", **kw); raised = False
    except ValueError:
        raised = True
    print(f"C5 method='lm' with {label} raises ValueError: {'held' if raised else 'FAILED'}")

# ---- reference minimum on the noisy image ----
ref_prob = Problem(noisy, PARAMS["P"])
r1 = lm(ref_prob, TRUE, "diag")[0]
r2 = least_squares(ref_prob.res, r1, jac=ref_prob.jac, method="lm", x_scale="jac",
                   ftol=1e-15, xtol=1e-15, gtol=1e-15).x
t_ref = fold(min((r1, r2), key=lambda t: cost(t, noisy)))
F_ref = cost(t_ref, noisy)
Jr = jacobian(t_ref)
SE = 10.0 * np.sqrt(np.diag(np.linalg.inv(Jr.T @ Jr)))
print("\nreference minimum   " + "  ".join(f"{n}={v:.4f}" for n, v in zip(NAMES, t_ref)))
print("standard errors     " + "  ".join(f"{n}={v:.4f}" for n, v in zip(NAMES, SE)))
print("|ref - truth| / SE  " + "  ".join(f"{n}={v:.2f}" for n, v in zip(NAMES, np.abs(t_ref - TRUE) / SE)))
print(f"F_ref = {F_ref:.2f}; 0.5 sigma^2 (m - n) = {0.5 * 100 * (441 - 5):.2f}")

# ---- all methods, both parameterizations, 100 starts ----
starts = [moment_start(noisy, dx, dy, s0) for dx in (-4, -2, 0, 2, 4)
          for dy in (-4, -2, 0, 2, 4) for s0 in (0.8, 1.6, 3.2, 6.4)]
runs = {}
for name, fit in METHODS.items():
    for key, sc in PARAMS.items():
        rows = []
        for u0 in starts:
            prob = Problem(noisy, sc)
            x, F, nfev, it, conv = fit(prob, u0 * sc)
            rows.append(dict(t=fold(x / sc), F=F, nfev=nfev, it=it, conv=conv, path=np.array(prob.points)))
        runs[name, key] = rows

print("\nmethod      par  conv reach else(min F/Fref)  nfev med  max  iters med  capped  "
      "max dF/Fref  max dtheta/SE")
h2 = True
for (name, key), rows in runs.items():
    conv = [r for r in rows if r["conv"]]
    reach = [r for r in rows if r["F"] < 1.01 * F_ref]
    good = [r for r in conv if r["F"] < 1.01 * F_ref]
    other = [r["F"] / F_ref for r in conv if r["F"] >= 1.01 * F_ref]
    dF = max(abs(r["F"] - F_ref) / F_ref for r in good) if good else np.nan
    dt = max(np.max(np.abs(r["t"] - t_ref) / SE) for r in good) if good else np.nan
    h2 &= bool(good) and dF < 1e-9 and dt < 0.01
    nf = [r["nfev"] for r in rows]
    print(f"{name:11s} {key:3s} {len(conv):4d} {len(reach):5d} {len(other):4d}"
          f"({min(other) if other else np.nan:9.3g})  {np.median(nf):8.1f} {max(nf):4d}"
          f"  {np.median([r['it'] for r in rows]):9.1f}  {sum(n >= 1000 for n in nf):6d}"
          f"  {dF:11.1e}  {dt:13.2e}")
print(f"min final cost over all runs / F_ref - 1 = "
      f"{min(r['F'] for rows in runs.values() for r in rows) / F_ref - 1:.1e} "
      f"(C3 bar >= -1e-12): {'held' if min(r['F'] for rows in runs.values() for r in rows) >= F_ref * (1 - 1e-12) else 'FAILED'}")

print("\ninvariance, P against P' (differences in standard errors)")
inv = {}
for name in METHODS:
    eq, maxd, differ, first = 0, 0.0, 0, []
    for a, b in zip(runs[name, "P"], runs[name, "P'"]):
        pa, pb = a["path"], b["path"]
        ia = next((k for k in range(1, len(pa)) if np.any(pa[k] != pa[0])), None)
        ib = next((k for k in range(1, len(pb)) if np.any(pb[k] != pb[0])), None)
        if ia is not None and ib is not None:          # first trial point after the start
            first.append(np.max(np.abs(pa[ia] - pb[ib]) / SE))
        if len(pa) == len(pb):
            eq += 1
            d = np.max(np.abs(pa - pb) / SE)
            maxd = max(maxd, d)
            differ += d > 0.1
        else:
            differ += 1
    inv[name] = (eq, maxd, differ)
    print(f"{name:11s} equal counts {eq:3d}/100  max path diff {maxd:.1e}  "
          f"paths differ {differ:3d}/100  min first-trial diff {min(first):.2e}")

med = {(m, k): np.median([r["nfev"] for r in runs[m, k]]) for m in ("MINPACK-jac", "MINPACK-1") for k in PARAMS}
print()
for h, name in (("H1a", "LM-diag"), ("H1b", "MINPACK-jac")):
    eq, maxd, _ = inv[name]
    print(f"{h} {name} invariant (100/100 equal, max diff < 1e-6): {'held' if eq == 100 and maxd < 1e-6 else 'failed'}")
print(f"H1c LM-I paths differ on 100/100: {'held' if inv['LM-I'][2] == 100 else 'failed'}")
for h, name in (("H1d", "MINPACK-1"), ("H1e", "TRF-1")):
    print(f"{h} {name} paths differ on >= 1 start: {'held' if inv[name][2] >= 1 else 'failed'}")
print(f"H2 every converged run below 1.01 F_ref is at the reference: {'held' if h2 else 'failed'}")
for h, key in (("H3a", "P"), ("H3b", "P'")):
    j, o = med["MINPACK-jac", key], med["MINPACK-1", key]
    print(f"{h} {key}: median nfev jac {j:.1f} <= x_scale=1 {o:.1f}: {'held' if j <= o else 'failed'}")

Output:

D1 noise SD = 10.064 counts (bar |SD - 10| < 1.0): held
D2 parameter magnitude ratio in P = 625.0 (bar 300-3000): held; Jacobian column-norm ratio at truth: P 883.9, P' 135.0
C1 noise-free LM-I        P : max rel error 5.1e-14, cost 5.2e-21: held
C1 noise-free LM-I        P': max rel error 1.4e-15, cost 3.9e-24: held
C1 noise-free LM-diag     P : max rel error 0.0e+00, cost 0.0e+00: held
C1 noise-free LM-diag     P': max rel error 0.0e+00, cost 0.0e+00: held
C1 noise-free MINPACK-jac P : max rel error 1.1e-16, cost 1.1e-25: held
C1 noise-free MINPACK-jac P': max rel error 0.0e+00, cost 0.0e+00: held
C1 noise-free MINPACK-1   P : max rel error 1.1e-16, cost 1.1e-25: held
C1 noise-free MINPACK-1   P': max rel error 0.0e+00, cost 0.0e+00: held
C1 noise-free TRF-1       P : max rel error 1.1e-16, cost 1.1e-25: held
C1 noise-free TRF-1       P': max rel error 0.0e+00, cost 0.0e+00: held
C2 Jacobian against central differences: 9.2e-11 (bar 1e-6): held
C4 'lm' default x_scale against 'jac': max |dx| = 0.0e+00, nfev 7 and 7: held
C5 method='lm' with loss='huber' raises ValueError: held
C5 method='lm' with bounds raises ValueError: held

reference minimum   A=1001.8015  x0=10.2922  y0=9.6048  s=1.5889  B=100.1898
standard errors     A=5.0218  x0=0.0080  y0=0.0080  s=0.0059  B=0.5146
|ref - truth| / SE  A=0.36  x0=0.98  y0=0.60  s=1.90  B=0.37
F_ref = 22008.37; 0.5 sigma^2 (m - n) = 21800.00

method      par  conv reach else(min F/Fref)  nfev med  max  iters med  capped  max dF/Fref  max dtheta/SE
LM-I        P    100    84   16(      169)      14.0  247       13.0       0      8.3e-16       9.81e-08
LM-I        P'    79    64   16(      166)      20.0 1000       15.0      20      9.9e-16       2.87e-07
LM-diag     P     78    64   14(      166)      23.0 1000       12.0      11      9.9e-16       3.76e-07
LM-diag     P'    77    64   13(      166)      23.0 1000       12.0      11      1.2e-15       3.76e-07
MINPACK-jac P     86    63   23(      166)      13.0 1000       10.0      14      2.6e-12       2.89e-05
MINPACK-jac P'    86    63   23(      166)      13.0 1000       10.0      14      2.6e-12       2.89e-05
MINPACK-1   P     83    61   22(      166)     393.0 1000      357.5      17      2.9e-12       2.54e-05
MINPACK-1   P'    81    66   15(      166)      14.0 1000       10.0      19      2.7e-12       2.17e-05
TRF-1       P     43    37    6(      166)    1000.0 1000      967.5      57      2.9e-12       2.63e-05
TRF-1       P'    84    67   17(      166)      13.0 1000       11.0      16      2.6e-12       2.35e-05
min final cost over all runs / F_ref - 1 = -1.1e-15 (C3 bar >= -1e-12): held

invariance, P against P' (differences in standard errors)
LM-I        equal counts   1/100  max path diff 1.4e+02  paths differ 100/100  min first-trial diff 1.05e+01
LM-diag     equal counts  96/100  max path diff 1.9e+157  paths differ  23/100  min first-trial diff 2.26e-14
MINPACK-jac equal counts  98/100  max path diff 6.9e+03  paths differ  16/100  min first-trial diff 3.79e-14
MINPACK-1   equal counts  15/100  max path diff 7.7e+04  paths differ  94/100  min first-trial diff 3.79e-14
TRF-1       equal counts  27/100  max path diff 4.3e+04  paths differ  91/100  min first-trial diff 2.76e-14

H1a LM-diag invariant (100/100 equal, max diff < 1e-6): failed
H1b MINPACK-jac invariant (100/100 equal, max diff < 1e-6): failed
H1c LM-I paths differ on 100/100: held
H1d MINPACK-1 paths differ on >= 1 start: held
H1e TRF-1 paths differ on >= 1 start: held
H2 every converged run below 1.01 F_ref is at the reference: held
H3a P: median nfev jac 13.0 <= x_scale=1 393.0: held
H3b P': median nfev jac 13.0 <= x_scale=1 14.0: held

Checks. The noise has sample standard deviation 10.064 (check D1 held) and the parameter magnitudes differ by a factor of 625 (D2 held); the Jacobian’s column norms at the truth differ by a factor of 883.9 in P and 135.0 in P′, so P is the more badly scaled problem. On the noise-free image every solver recovers the true parameters in both parameterizations (C1 held in all ten rows), the analytic Jacobian matches central differences (C2 held), no run ends below the reference cost (C3 held), SciPy 1.16.3’s method='lm' without x_scale gives exactly the result of x_scale='jac' (C4 held), and method='lm' refuses a Huber loss and bounds (C5 held). The reference minimum lies within 1.90 standard errors of the truth in every parameter, and its cost, 22,008.37, is close to $\tfrac{1}{2}\sigma^2(m - n) = 21{,}800$, approximately the value expected for a correct model at the true noise level (our addition).

The same minimum (H2 held). In every row, every converged run that ends near the reference cost agrees with the reference to within $2.9 \times 10^{-12}$ in relative cost and $2.89 \times 10^{-5}$ standard errors in every parameter, against bars of $10^{-9}$ and 0.01. The hand-written solvers agree more closely (at most $1.2 \times 10^{-15}$ and $3.76 \times 10^{-7}$), which in our reading follows from their tighter stopping thresholds ($10^{-12}$ and $10^{-10}$ against SciPy’s $10^{-8}$). Damping and scaling change the path and the evaluation count, not the estimate. Every solver also stops by its own tests, from some starts, at a point of higher cost (the “else” column, whose smallest cost ratio is 166 in every row except LM-I in P, 169), and the number of the 100 starts that reach the minimum (the “reach” column) depends on the solver and the parameterization; Nonlinear Least Squares’ circle example studies such basins.

Evaluation counts (H3a and H3b held). MINPACK with column-norm scaling needs a median of 13 evaluations in both parameterizations. With x_scale=1.0 it needs 14 in P′ but 393 in P, and TRF-1, whose trust region is also unscaled, needs a median of 1,000 in P, hitting the budget on 57 of 100 starts, against 13 in P′. Both rows of MINPACK-jac print identical summaries; the two LM-diag rows differ by one converged run that ends away from the reference (78 against 77 converged, 14 against 13 in the “else” column) and in the cost-agreement column. LM-I is the exception that the derivation predicts: its two rows differ (reach 84 in P against 64 in P′) because in P′ it is LM with a different damping matrix and a different initial damping, by construction.

Invariance (H1a and H1b failed; H1c, H1d and H1e held). The planned check of the derivation was that, for LM-diag (H1a) and MINPACK-jac (H1b), every start would give the same number of evaluations in both parameterizations and paths that agree to within $10^{-6}$ standard errors. Both checks failed. LM-diag gave equal counts on 96 of 100 starts, a largest path difference of $1.9 \times 10^{157}$ standard errors, and paths that differ by more than 0.1 standard errors, or in length, on 23 of 100 starts. MINPACK-jac gave equal counts on 98 of 100 starts, a largest difference of $6.9 \times 10^3$ standard errors, and differing paths on 16 of 100. The smallest first-trial-point differences for the two are $2.26 \times 10^{-14}$ and $3.79 \times 10^{-14}$ standard errors. The three solvers predicted not to be invariant are not: LM-I’s paths differ on all 100 starts (H1c held), with even its first trial points at least 10.5 standard errors apart, and MINPACK-1’s and TRF-1’s paths differ on 94 and 91 starts (H1d and H1e held, each asking for at least one start). The smallest first-trial differences of MINPACK-1 and TRF-1, $3.79 \times 10^{-14}$ and $2.76 \times 10^{-14}$ standard errors, show that on at least one start their first trial point is the same in both units, which in our reading is the undamped Gauss–Newton step that the derivation predicts to be invariant.

Process note. The set-up, checks, hypotheses and every bar were written before any script ran. A planning script then computed the values the bars rely on, at one start of the grid ($\Delta x_0 = 2$, $\Delta y_0 = -2$, $s_0 = 3.2$): there, the full paths of LM-diag and MINPACK-jac agreed between the parameterizations to $8.33 \times 10^{-13}$ and $1.63 \times 10^{-12}$ standard errors, LM-I’s first step differed by 180 standard errors, and every solver reached the reference to within $9.58 \times 10^{-13}$ in relative cost and $1.72 \times 10^{-5}$ standard errors. D1, D2, C1 and C2 were also computed there and held. H1’s bars were therefore fixed in the frozen plan only after the invariance had been seen to hold on one start, so H1 is an illustration of the derivation rather than a test, and it failed on other starts; H2 had also been seen on that start; H3 was not computed before the pilot and is a test. No bar changed after the planning script. A pilot ran the script unchanged and the first run, which is shown, repeated it with identical output. The example departs from the batch plan’s description in eight ways: (1) the truth, noise level and single noise draw are fixed as above; (2) the start grid and the moment-based amplitude and background are our choice; (3) the hand-written LM’s update, initial damping and stopping rules are ours, and MINPACK’s rules are described, not ported; (4) every SciPy setting is explicit (analytic Jacobian, tolerances $10^{-8}$, a budget of 1,000 evaluations, x_scale=1.0 for 'trf'); (5) the measures are made operational (a reference minimum computed in the script, “reach” as a final cost below 1.01 times the reference, convergence as each method’s own stopping test, path differences in standard errors on every evaluated point, the error against the truth printed once for the reference, $\lvert s\rvert$ in comparisons); (6) the invariance hypothesis is split into H1a–H1e and extended to the three SciPy solvers, H2 is restricted to converged runs that end near the reference, and H3 is scored in both parameterizations; (7) data and exact-data checks D1, D2 and C1–C5 are added; (8) the fraction of starts reaching the minimum is printed without a hypothesis.

Properties and Behavior

• Descent and two limits. Nonlinear Least Squares shows that the damped step is a descent direction and states the Gauss–Newton rates; since a rejected step is never kept, the cost of the accepted iterates decreases. With a scaling matrix, the large-$\lambda$ limit is $-(D^\top D)^{-1}\mathbf{g}/\lambda$, a short scaled-gradient step.
• Same estimator, different paths. Changing the damping matrix, the update rule or the library changes which starts succeed and how many evaluations are spent, not the minimum reached, as H2 in the example shows to within $3 \times 10^{-12}$ in relative cost.
• Scale invariance is a property of the derivation. Diagonal scaling from the Jacobian removes the dependence on parameter units in the derivation; in the example the invariance checks H1a and H1b failed, with the paths of LM-diag and MINPACK-jac differing between the parameterizations on 23 and 16 of 100 starts.
• Ends other than the minimum. Converged runs can stop at points of higher cost; in the example every solver did so from some starts (the “else” column). Comparing the final cost with $\tfrac{1}{2}\sigma^2(m-n)$ is a cheap check when the noise level is known.

Limitations

• Local. LM refines a good start; it does not find one. Vision pipelines start it from a linear method, as camera calibration does with Zhang’s closed form before the refinement that the Zhang paper describes as Minpack’s Levenberg–Marquardt.
• Squared loss. Outliers pull the estimate; robust losses need reweighting or a solver that supports them (SciPy’s 'trf', Ceres’ loss functions), and robust estimation treats the estimators.
• Unbounded parameters. The basic method has no bounds; widths, focal lengths or variances that must stay positive need a reparameterization or a bounded solver.
• Parameters that are not vectors. Rotations need a local parametrization of the update, as pose graph optimization and bundle adjustment describe.
• Dense linear algebra. A dense factorization of $J$ costs $O(mn^2)$ per iteration, which rules out large problems without exploiting sparsity.
• Tuning knobs. The initial damping, the update rule and the tolerances are heuristics; libraries differ in all three, so the same problem can take very different numbers of evaluations, as the x_scale=1.0 rows of the example show.

Variants

• Powell’s dog leg. A trust-region method that combines the Gauss–Newton step and the steepest-descent (Cauchy) step instead of solving for $\lambda$; Ceres offers it with traditional and subspace variants, and SciPy’s 'dogbox' applies it within rectangular trust regions.
• Trust-region reflective. SciPy’s default 'trf' solves trust-region subproblems that respect bounds, with an exact solver similar to MINPACK’s for dense Jacobians, according to its documentation; in common.py its radius is set to a quarter of the step length when the ratio is below 0.25 and doubled when the ratio exceeds 0.75 and the step reached the boundary.
• Inexact-step LM. The damped system is solved approximately by an iterative method such as conjugate gradients, as in Ceres when an iterative linear solver is chosen.
• Sparse and structured LM. Bundle adjustment eliminates the points with the Schur complement before each damped solve; libraries such as SBA, Ceres and g2o implement this, as bundle adjustment describes, and reprojection error gives the cost they minimize.
• Robust LM. A robust loss is handled by rescaling residuals and Jacobian rows so that the gradient and Gauss–Newton Hessian match those of the robust cost, as SciPy’s least_squares documentation describes for its non-'lm' methods; Nonlinear Least Squares explains the reweighting.