Gauss-Newton and Levenberg-Marquardt for Nonlinear Least Squares

Intuition Beginner

A great many problems boil down to the same request: pick the numbers that make a model match a pile of measurements as closely as possible. You have a curve with a few dials on it, you have data points, and for each data point there is a residual — the gap between what the model predicts and what was actually measured. The goal is to turn the dials so the gaps are as small as possible overall, measured by adding up the squares of all the gaps.

If the dials entered the model in a straight-line way, this would be the ordinary least-squares fit you already know: one clean solve. The catch is that the dials usually enter in a bent, curved way. So the trick is to pretend, for one step, that things are straight. Stand at your current dial settings, replace each curved residual by its best straight-line stand-in, solve the easy fit that results, and move the dials there. Then repeat from the new spot. This is the Gauss-Newton idea: a sequence of ordinary least-squares solves that chase the true bent problem.

Why not just use the general curvature-aware method from the previous units? Because least-squares problems hand you a gift. The very same slopes you compute to build the straight-line stand-in also give you most of the curvature information for free, with no extra second-derivative work. You exploit the structure instead of ignoring it.

There is one place this gift runs thin. When the model can match the data well, the leftover gaps are small and the shortcut is excellent. When the model is a poor fit and the gaps stay large no matter how you turn the dials, the part of the curvature you threw away starts to matter, and the plain method slows down or wanders. The fix is to add a brake — a damping dial — that shortens and steadies each step when the shortcut is misbehaving, and releases when it is working well.

Visual Beginner

Picture fitting a curve through scattered dots. At your current guess the curve sits a bit off; each dot has a short vertical line connecting it to the curve — that line is its residual. Gauss-Newton looks at how each residual would change if you nudged the dials, builds a flat (straight-line) picture of all those changes, and solves for the nudge that shrinks the squared gaps the most.

The damped version, Levenberg-Marquardt, adds a brake. With the brake off, you take the full straight-line-fit step, which is bold and fast when the fit is good. With the brake on hard, you take a short, cautious step pointed almost straight downhill. The method turns the brake automatically: a step that worked releases the brake, a step that overshot tightens it.

situation	what the method does	brake setting
good fit, small gaps	full Gauss-Newton step, fast	loose (small damping)
step improved the fit	accept, loosen the brake	decreases
step made things worse	reject, tighten the brake	increases
brake very tight	short step, nearly straight downhill	large damping

The takeaway: Gauss-Newton turns a bent fitting problem into a sequence of straight-line fits, and the damping brake of Levenberg-Marquardt keeps those steps trustworthy when the shortcut would otherwise mislead.

Worked example Beginner

Fit the one-dial model $\varphi (t;x)=xt$ — a line through the origin with slope $x$ — to two points: at $t=1$ the value is $3$, and at $t=2$ the value is $5$. The residuals are the model-minus-data gaps,

$$r_1(x)=x\cdot 1-3=x-3,r_2(x)=x\cdot 2-5=2x-5.$$

Nudging $x$ by one unit changes $r_1$ by $1$ and $r_2$ by $2$; those two numbers are the slopes (the Jacobian). Start at $x=0$, where the residuals are $r_1=-3$ and $r_2=-5$.

The Gauss-Newton step solves the straight-line fit: find the nudge $p$ that makes $1\cdot p+(-3)$ and $2\cdot p+(-5)$ as small as possible in squared total. The recipe (the normal equation) divides the slope-weighted residual sum by the slope-squared sum:

$$p=-\frac{(1)(-3)+(2)(-5)}{(1{)}^2+(2{)}^2}=-\frac{-3-10}{1+4}=-\frac{-13}{5}=\mathrm{2.6.}$$

So the new dial setting is $x=0+2.6=2.6$.

Because this model is straight-line in $x$, the slopes never change, so one step lands on the exact best fit. Check it: at $x=2.6$ the residuals are $r_1=2.6-3=-0.4$ and $r_2=5.2-5=0.2$. The total squared gap is $(-0.4{)}^2+(0.2{)}^2=0.16+0.04=0.20$. Trying a nearby value like $x=2.5$ gives $(-0.5{)}^2+(0{)}^2=0.25$, which is larger — so $2.6$ really is better.

What this tells us: Gauss-Newton converts each fitting step into a weighting of residuals by their slopes, and when the model is straight in its dials it finishes in a single step, exactly like ordinary least squares.

Check your understanding Beginner

Formal definition Intermediate+

Let $r:{\mathbb{R}}^n\to {\mathbb{R}}^m$ with $m\ge n$ be a twice continuously differentiable residual map, $r(x)=(r_1(x),\dots ,r_m(x))$, and consider the nonlinear least-squares objective

$$f(x)=\frac{1}{2}\Vert r(x){\Vert }^2=\frac{1}{2}\sum _{j=1}^m{r}_j(x{)}^2.$$

Write $J(x)=[\partial r_j/\partial x_i{]}_{j,i}\in {\mathbb{R}}^{m\times n}$ for the Jacobian of $r$. Differentiating $f$ gives the gradient and Hessian in the structured form

$$\nabla f(x)=J(x{)}^{\top }r(x),{\nabla }^{2}f(x)=J(x{)}^{\top }J(x)+\sum _{j=1}^m{r}_j(x){\nabla }^2{r}_j(x)=:J^{\top }J+S(x),$$

where $S(x)={\sum }_j{r}_j(x){\nabla }^2{r}_j(x)$ is the residual-curvature term. The first piece $J^{\top }J$ uses only first derivatives — the same data as the gradient — while $S(x)$ requires the second derivatives of each residual and is weighted by the residuals themselves.

The Gauss-Newton step $p^{\mathrm{GN}}$ drops $S(x)$, approximating ${\nabla }^{2}f$ by $J^{\top }J$, and solves

$$J(x{)}^{\top }J(x)p^{\mathrm{GN}}=-J(x{)}^{\top }r(x)\text{(the normal equations [43.04.01])}.$$

Equivalently, $p^{\mathrm{GN}}$ minimises the linearised residual $\frac{1}{2}\Vert J(x)p+r(x){\Vert }^2$: it is the linear least-squares solution of ${\min }_p\Vert J(x)p+r(x){\Vert }^2$, the exact subproblem one would solve if $r$ were replaced by its first-order model $r(x)+J(x)p$. When $J(x)$ has full column rank the normal-equation matrix $J^{\top }J$ is positive definite, $p^{\mathrm{GN}}$ is unique, and $\nabla f^{\top }{p}^{\mathrm{GN}}=-(J^{\top }r{)}^{\top }(J^{\top }J{)}^{-1}(J^{\top }r)\le 0$, so $p^{\mathrm{GN}}$ is a descent direction whenever $\nabla f\ne 0$. In practice the normal equations are not formed; instead a QR factorisation $J=QR$ 43.04.01 reduces ${\min }_p\Vert Jp+r\Vert$ to the triangular solve $Rp=-Q_1^{\top }r$, halving the condition number relative to $J^{\top }J$. The Gauss-Newton iteration is $x_{k+1}=x_k+{\alpha }_k{p}_k^{\mathrm{GN}}$ with ${\alpha }_k$ from a line search 44.03.01.

The Levenberg-Marquardt step damps the normal equations,

$$(J(x{)}^{\top }J(x)+\lambda I)p^{\mathrm{LM}}=-J(x{)}^{\top }r(x),\lambda \ge 0.$$

Since $J^{\top }J⪰0$, the shift $\lambda I$ makes $J^{\top }J+\lambda I\succ 0$ for every $\lambda >0$, so $p^{\mathrm{LM}}$ exists and is a descent direction even when $J$ is rank-deficient. As $\lambda \to 0$, $p^{\mathrm{LM}}\to p^{\mathrm{GN}}$; as $\lambda \to \infty$, $p^{\mathrm{LM}}\to -\frac{1}{\lambda }{J}^{\top }r$, a short step along the steepest-descent direction $-\nabla f$. The damping $\lambda$ is the Lagrange multiplier of a trust region 44.03.02: $p^{\mathrm{LM}}$ solves ${\min }_{\Vert p\Vert \le \Delta }\Vert Jp+r{\Vert }^2$ for a radius $\Delta$ tied to $\lambda$ by the trust-region characterisation $(J^{\top }J+\lambda I)p=-J^{\top }r$, $\lambda (\Delta -\Vert p\Vert )=0$. The parameter $\lambda$ is updated by the gain ratio

$$\rho =\frac{f(x)-f(x+p^{\mathrm{LM}})}{\frac{1}{2}\Vert r(x){\Vert }^2-\frac{1}{2}\Vert J(x)p^{\mathrm{LM}}+r(x){\Vert }^2},$$

the ratio of actual to predicted decrease: decrease $\lambda$ (loosen damping) when $\rho$ is large, increase $\lambda$ (tighten) when $\rho$ is small or negative. The new symbols $r$ (residual), $J$ (Jacobian), $S$ (residual-curvature term), $p^{\mathrm{GN}}$, $p^{\mathrm{LM}}$, and $\lambda$ (LM damping) are recorded in _meta/NOTATION.md.

Counterexamples to common slips Intermediate+

• "Gauss-Newton is Newton's method for $f$." It is Newton applied to the approximate Hessian $J^{\top }J$, omitting $S(x)=\sum r_j{\nabla }^2{r}_j$. The two coincide only when $S$ vanishes — at a zero-residual solution, or when every residual is affine in $x$. With large residuals at the solution the methods differ and Gauss-Newton loses Newton's quadratic rate.

• "$J^{\top }J$ is always invertible, so the Gauss-Newton step always exists." Only when $J$ has full column rank. A rank-deficient or ill-conditioned $J$ makes $J^{\top }J$ singular or nearly so, and the Gauss-Newton step is undefined or wildly large. This is precisely the failure Levenberg-Marquardt repairs through the shift $\lambda I$.

• "Levenberg-Marquardt is just Gauss-Newton with a fixed regularisation." The whole method is in the adaptive control of $\lambda$ via the gain ratio. A fixed $\lambda$ either over-damps (slow) or under-damps (unstable); LM converges because it lets $\lambda \to 0$ near a good solution (recovering the Gauss-Newton rate) and raises $\lambda$ when steps fail.

• "Form the normal equations $J^{\top }Jp=-J^{\top }r$ and solve them." Forming $J^{\top }J$ squares the condition number, $\kappa (J^{\top }J)=\kappa (J{)}^2$, and can lose half the available digits. The stable route is the QR factorisation of $J$ (or, for LM, of the stacked matrix $\left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]$), which works with $\kappa (J)$ directly.

Key theorem with proof Intermediate+

The signature result is the local convergence dichotomy of Gauss-Newton: it explains both the method's speed on well-fitting models and its degradation on poorly fitting ones, isolating the omitted residual-curvature term $S(x^{\ast })$ as the single quantity that decides the rate.

Theorem (local convergence of Gauss-Newton, small- versus large-residual). Let $x^{\ast }$ be a solution with $\nabla f(x^{\ast })=J(x^{\ast }{)}^{\top }r(x^{\ast })=0$, $J(x^{\ast })$ of full column rank, $r$ twice continuously differentiable with Lipschitz Jacobian near $x^{\ast }$, and write $H^{\ast }=J(x^{\ast }{)}^{\top }J(x^{\ast })\succ 0$ and $S^{\ast }={\sum }_j{r}_j(x^{\ast }){\nabla }^2{r}_j(x^{\ast })$. Then the undamped Gauss-Newton iteration $x_{k+1}=x_k-(J_k^{\top }{J}_k{)}^{-1}{J}_k^{\top }{r}_k$ from a start near $x^{\ast }$ satisfies

$$\Vert x_{k+1}-x^{\ast }\Vert \le \Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert \Vert x_k-x^{\ast }\Vert +C\Vert x_k-x^{\ast }{\Vert }^2$$

for a constant $C$. Hence if $\Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert <1$ the iteration converges linearly with asymptotic rate $\Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert$; if in addition $r(x^{\ast })=0$ (the zero-residual case) then $S^{\ast }=0$ and the convergence is quadratic.

Proof. Write $g(x)=\nabla f(x)=J(x{)}^{\top }r(x)$, $H(x)={\nabla }^{2}f(x)=J(x{)}^{\top }J(x)+S(x)$, $e_k=x_k-x^{\ast }$, and $G(x)=J(x{)}^{\top }J(x)$ for the Gauss-Newton matrix. Full column rank of $J(x^{\ast })$ and continuity make $G(x)$ invertible with $\Vert G(x{)}^{-1}\Vert \le 2\Vert G(x^{\ast }{)}^{-1}\Vert$ on a ball about $x^{\ast }$, so the step exists there. Since $g(x^{\ast })=0$,

$$e_{k+1}=e_k-G(x_k{)}^{-1}g(x_k)=G(x_k{)}^{-1}(G(x_k)e_k-(g(x_k)-g(x^{\ast }))).$$

The fundamental theorem of calculus applied to $g$ along the segment from $x^{\ast }$ to $x_k$ 44.03.03 gives $g(x_k)-g(x^{\ast })={\int }_0^{1}H(x^{\ast }+t{e}_k)e_{k}dt$, so

$$G(x_k)e_k-(g(x_k)-g(x^{\ast }))={\int }_0^1(G(x_k)-H(x^{\ast }+t{e}_k))e_{k}dt.$$

Split the integrand: $G(x_k)-H(x^{\ast }+t{e}_k)=(G(x_k)-G(x^{\ast }+t{e}_k))-S(x^{\ast }+t{e}_k)$, using $H=G+S$. The first bracket has norm $O(\Vert e_k\Vert )$ because $G$ is Lipschitz (a product of Lipschitz Jacobians), contributing the $C\Vert e_k{\Vert }^2$ term after the outer factor of $\Vert e_k\Vert$. The second bracket converges to $S^{\ast }$ as $e_k\to 0$, contributing $\Vert G(x_k{)}^{-1}\Vert \cdot \Vert S^{\ast }\Vert \cdot \Vert e_k\Vert +o(\Vert e_k\Vert )$. Collecting,

$$\Vert e_{k+1}\Vert \le \Vert G(x_k{)}^{-1}\Vert (\Vert S^{\ast }\Vert +o(1))\Vert e_k\Vert +C\Vert e_k{\Vert }^2\le \Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert \Vert e_k\Vert +C^{\prime }\Vert e_k{\Vert }^2,$$

since $G(x^{\ast })=H^{\ast }$ (the matrices agree at $x^{\ast }$ because $S$ multiplies $H^{-1}$ via $G^{-1}\to (H^{\ast }{)}^{-1}$). When $r(x^{\ast })=0$ every $r_j(x^{\ast })=0$ so $S^{\ast }=0$, the linear term drops, and $\Vert e_{k+1}\Vert \le C^{\prime }\Vert e_k{\Vert }^2$ is quadratic. $\square$

Bridge. This dichotomy is the foundational reason Gauss-Newton is the method of choice for least squares yet must be globalised for hard fits: the rate is set entirely by $\Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert$, the size of the curvature the method throws away relative to the curvature it keeps, and this is exactly the term that distinguishes Gauss-Newton from the full Newton method 44.03.03 of the previous unit, whose Hessian retains $S(x)$ and so always converges quadratically. The construction builds toward Levenberg-Marquardt, which appears again in the Advanced results as the trust-region 44.03.02 regularisation that supplies global convergence when $\Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert \ge 1$ makes the undamped iteration diverge. The central insight is that the residual structure ${\nabla }^{2}f=J^{\top }J+S$ generalises the linear least-squares normal equations 43.04.01 — at a zero-residual solution the nonlinear problem is locally a linear least-squares problem, $J^{\top }J$ is the exact Hessian, and the QR solve of ${\min }_p\Vert Jp+r\Vert$ inherits the linear theory verbatim. Putting these together, the damping parameter $\lambda$ is dual to the trust-region radius $\Delta$, and the same gain-ratio test that drove the radius update of 44.03.02 now drives $\lambda$, so the globalisation is the previous unit's trust-region machinery specialised to the Gauss-Newton model.

Exercises Intermediate+

Advanced results Master

The Gauss-Newton model is one approximation — replace ${\nabla }^{2}f=J^{\top }J+S$ by $J^{\top }J$ — and the results below organise its theory around three axes: the descent and convergence consequences of the omitted term $S$, the Levenberg-Marquardt damping that realises the Gauss-Newton model inside a trust region and so converges globally, and the separable structure that lets a subset of parameters be eliminated in closed form.

Theorem 1 (Gauss-Newton step, descent, and the small-residual rate). For full-column-rank $J(x_k)$ the Gauss-Newton step $p_k^{\mathrm{GN}}=-(J_k^{\top }{J}_k{)}^{-1}{J}_k^{\top }{r}_k$ is the unique minimiser of $\Vert J_{k}p+r_k{\Vert }^2$ and satisfies $\nabla f_k^{\top }{p}_k^{\mathrm{GN}}=-\Vert J_k(J_k^{\top }{J}_k{)}^{-1}{J}_k^{\top }{r}_k{\Vert }^2\le 0$, so it is a descent direction whenever $J_k^{\top }{r}_k\ne 0$; the quantity $J_k(J_k^{\top }{J}_k{)}^{-1}{J}_k^{\top }$ is the orthogonal projector onto the column space of $J_k$, so the directional derivative equals $-\Vert {\Pi }_{\mathrm{col}{J}_k}{r}_k{\Vert }^2$, the squared length of the residual's projection onto the reachable directions 43.04.01. At a solution $x^{\ast }$ with $\nabla f(x^{\ast })=0$, $J(x^{\ast })$ full rank, and Lipschitz Jacobian, the iteration converges locally with asymptotic rate $\Vert (J^{\ast }{}^{\top }{J}^{\ast }{)}^{-1}{S}^{\ast }\Vert$, which is $0$ (quadratic convergence) in the zero-residual case and small (fast linear) for nearly-linear or well-fitting models; when $\Vert (J^{\ast }{}^{\top }{J}^{\ast }{)}^{-1}{S}^{\ast }\Vert \ge 1$ the undamped iteration diverges ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 2 (Levenberg-Marquardt as the Gauss-Newton trust-region subproblem). The LM step $p_k^{\mathrm{LM}}$ solving $(J_k^{\top }{J}_k+{\lambda }_{k}I)p=-J_k^{\top }{r}_k$ is exactly the solution of the trust-region subproblem ${\min }_{\Vert p\Vert \le {\Delta }_k}\Vert J_{k}p+r_k{\Vert }^2$ for the radius ${\Delta }_k=\Vert p_k^{\mathrm{LM}}({\lambda }_k)\Vert$, with ${\lambda }_k$ the trust-region multiplier satisfying ${\lambda }_k({\Delta }_k-\Vert p_k^{\mathrm{LM}}\Vert )=0$ and $J_k^{\top }{J}_k+{\lambda }_{k}I⪰0$ 44.03.02. Because $J^{\top }J⪰0$ has its smallest eigenvalue $\ge 0$, the LM "hard case" never arises: $J^{\top }J+\lambda I\succ 0$ for every $\lambda >0$, so the step is always well defined and a descent direction even at rank-deficient $J_k$. As $\lambda$ ranges over $[0,\infty )$ the step interpolates monotonically (Exercise 6) from the full Gauss-Newton step $p^{\mathrm{GN}}$ down to the infinitesimal steepest-descent step $-\frac{1}{\lambda }{J}^{\top }r$, so LM is a blend of Gauss-Newton's speed and steepest descent's robustness selected by the data. Global convergence to a stationary point follows from the trust-region theory of 44.03.02: the predicted reduction $\frac{1}{2}\Vert r{\Vert }^2-\frac{1}{2}\Vert Jp+r{\Vert }^2$ satisfies the Cauchy sufficient-decrease bound, and the gain-ratio update of $\lambda$ is the radius update in disguise ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 3 (damping update and the Marquardt heuristic). The classical Levenberg-Marquardt control of $\lambda$ uses the gain ratio ${\rho }_k={\mathrm{ared}}_k/{\mathrm{pred}}_k$: if ${\rho }_k>\frac{3}{4}$ the model is accurate, accept the step and set ${\lambda }_{k+1}={\lambda }_k/\nu$ (loosen, $\nu >1$); if ${\rho }_k<\frac{1}{4}$ the model overpredicted, set ${\lambda }_{k+1}=\nu {\lambda }_k$ (tighten) and, if ${\rho }_k\le 0$, reject the step; otherwise hold $\lambda$. Marquardt's 1963 refinement scales the shift by the diagonal of $J^{\top }J$ rather than the identity, solving $(J^{\top }J+\lambda D)p=-J^{\top }r$ with $D=\mathrm{diag}(J^{\top }J)$, which makes the damping affine-invariant under coordinate rescaling and matches the local curvature of each parameter — important when parameters have widely different scales. Near a small-residual solution ${\rho }_k\to 1$ so ${\lambda }_k\to 0$, and LM recovers the Gauss-Newton (quadratic, zero-residual) rate of Theorem 1; far away, large ${\lambda }_k$ keeps steps short and safe ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 4 (separable least squares and variable projection). Suppose the residual is linear in a subset of parameters: $r(x,y)=\Phi (y)x-d$ where $x\in {\mathbb{R}}^{n_1}$ enters linearly through the matrix $\Phi (y)\in {\mathbb{R}}^{m\times n_1}$ and $y\in {\mathbb{R}}^{n_2}$ enters nonlinearly. For fixed $y$ the optimal linear part is the linear least-squares solution $x(y)=\Phi (y{)}^{+}d$ (the pseudoinverse), so the problem reduces to the variable-projection functional

$$\psi (y)=\frac{1}{2}\Vert (I-\Phi (y)\Phi (y{)}^{+})d{\Vert }^2=\frac{1}{2}\Vert {\Pi }_{\Phi (y)}^{\perp }d{\Vert }^2,$$

the squared norm of the residual projected off the column space of $\Phi (y)$. Minimising $\psi$ over $y$ alone — a smaller nonlinear problem — recovers the joint minimiser, and the Golub-Pereyra derivative formula gives the Jacobian of the projected residual in terms of $\Phi$, ${\Phi }^{+}$, and $\partial \Phi /\partial y_k$, so Gauss-Newton or LM applies to $\psi (y)$ directly. Variable projection typically converges in fewer iterations and is less sensitive to the starting guess for the linear parameters than treating $(x,y)$ jointly, because the linear block is solved exactly at every step ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Synthesis. Nonlinear least squares is one structure — the objective $f=\frac{1}{2}\Vert r{\Vert }^2$ with $\nabla f=J^{\top }r$ and ${\nabla }^{2}f=J^{\top }J+S$ — and the foundational reason the whole subject is tractable is that the first Hessian piece $J^{\top }J$ is free from the gradient computation, so the Gauss-Newton step is a linear least-squares solve 43.04.01 disguised as a Newton step. The central insight is the division of labour set by the omitted term $S$: when residuals are small or nearly linear, $S$ is negligible, $J^{\top }J$ is the true Hessian, and Gauss-Newton matches the Newton method 44.03.03 with its quadratic rate; when residuals are large, $S$ dominates and the bare iteration diverges, which is exactly the regime where Levenberg-Marquardt's damping is not optional but structural. The damping parameter $\lambda$ is dual to the trust-region radius $\Delta$ of 44.03.02: the LM step solves the trust-region subproblem of the Gauss-Newton model, $J^{\top }J⪰0$ guarantees the hard case never arises, and the gain-ratio control of $\lambda$ is the radius update specialised to least squares. Putting these together, Gauss-Newton generalises the linear normal equations to the nonlinear setting, Levenberg-Marquardt generalises trust-region globalisation to the Gauss-Newton model, and variable projection generalises the exact-elimination of linear parameters to the separable case — three specialisations of one residual structure, and the bridge from the linear least-squares chapter to the large-scale curve-fitting and bundle-adjustment problems where this method is the workhorse.

Full proof set Master

Proposition 1 (structured gradient and Hessian). For $f(x)=\frac{1}{2}\Vert r(x){\Vert }^2$ with $r\in C^2$ and Jacobian $J(x)$, one has $\nabla f(x)=J(x{)}^{\top }r(x)$ and ${\nabla }^{2}f(x)=J(x{)}^{\top }J(x)+{\sum }_{j=1}^m{r}_j(x){\nabla }^2{r}_j(x)$.

Proof. By the chain rule, $\partial f/\partial x_i={\sum }_j{r}_j\partial r_j/\partial x_i=(J^{\top }r{)}_i$, giving $\nabla f=J^{\top }r$. Differentiating again,

$$\frac{{\partial }^{2}f}{\partial x_i\partial x_l}=\sum _j(\frac{\partial r_j}{\partial x_i}\frac{\partial r_j}{\partial x_l}+r_j\frac{{\partial }^2{r}_j}{\partial x_i\partial x_l})=(J^{\top }J{)}_{il}+\sum _j{r}_j({\nabla }^2{r}_j{)}_{il},$$

where the first sum is the $(i,l)$ entry of $J^{\top }J$ and the second is the residual-weighted curvature. Assembling over $i,l$ gives ${\nabla }^{2}f=J^{\top }J+{\sum }_j{r}_j{\nabla }^2{r}_j$. $\square$

Proposition 2 (Gauss-Newton descent and the projection identity). If $J$ has full column rank and $J^{\top }r\ne 0$, then $p^{\mathrm{GN}}=-(J^{\top }J{)}^{-1}{J}^{\top }r$ satisfies $\nabla f^{\top }{p}^{\mathrm{GN}}=-\Vert \Pi r{\Vert }^2<0$, where $\Pi =J(J^{\top }J{)}^{-1}{J}^{\top }$ is the orthogonal projector onto $\mathrm{col}J$.

Proof. $\nabla f^{\top }{p}^{\mathrm{GN}}=(J^{\top }r{)}^{\top }(-(J^{\top }J{)}^{-1}{J}^{\top }r)=-(J^{\top }r{)}^{\top }(J^{\top }J{)}^{-1}(J^{\top }r)$. Insert $J^{\top }J(J^{\top }J{)}^{-1}=I$: this equals $-r^{\top }J(J^{\top }J{)}^{-1}{J}^{\top }J(J^{\top }J{)}^{-1}{J}^{\top }r=-r^{\top }{\Pi }^{\top }\Pi r=-\Vert \Pi r{\Vert }^2$, using ${\Pi }^{\top }=\Pi$ and ${\Pi }^2=\Pi$ (the matrix $\Pi =J(J^{\top }J{)}^{-1}{J}^{\top }$ is symmetric and idempotent, hence the orthogonal projector onto $\mathrm{col}J$ 43.04.01). It is negative unless $\Pi r=0$, i.e. unless $J^{\top }r=J^{\top }\Pi r=0$, contradicting $J^{\top }r\ne 0$. $\square$

Proposition 3 (local convergence rate). Under $\nabla f(x^{\ast })=0$, $J(x^{\ast })$ full column rank, and $r\in C^2$ with Lipschitz Jacobian near $x^{\ast }$, the Gauss-Newton iterates satisfy $\Vert e_{k+1}\Vert \le \Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert \Vert e_k\Vert +C\Vert e_k{\Vert }^2$ with $H^{\ast }=J^{\ast }{}^{\top }{J}^{\ast }$, $S^{\ast }={\sum }_j{r}_j(x^{\ast }){\nabla }^2{r}_j(x^{\ast })$, $e_k=x_k-x^{\ast }$.

Proof. As in the Key theorem. With $G(x)=J(x{)}^{\top }J(x)$, $g=\nabla f$, full rank gives $\Vert G(x{)}^{-1}\Vert$ bounded near $x^{\ast }$ and $e_{k+1}=G(x_k{)}^{-1}(G(x_k)e_k-(g(x_k)-g(x^{\ast })))$. The identity $g(x_k)-g(x^{\ast })={\int }_0^1{\nabla }^{2}f(x^{\ast }+t{e}_k)e_{k}dt$ and ${\nabla }^{2}f=G+S$ split the integrand into $(G(x_k)-G(x^{\ast }+t{e}_k))$, of norm $O(\Vert e_k\Vert )$ by Lipschitz $G$, and $-S(x^{\ast }+t{e}_k)\to -S^{\ast }$. Hence $\Vert e_{k+1}\Vert \le \Vert G(x^{\ast }{)}^{-1}\Vert \Vert S^{\ast }\Vert \Vert e_k\Vert +C\Vert e_k{\Vert }^2$, and $G(x^{\ast })=H^{\ast }$. The linear coefficient is $\Vert (H^{\ast }{)}^{-1}{S}^{\ast }\Vert$ (up to the working-ball factor on $\Vert G^{-1}\Vert$). If $r(x^{\ast })=0$ then $S^{\ast }=0$ and the bound is purely quadratic. $\square$

Proposition 4 (LM step is the augmented least-squares solution and is always a descent direction for $\lambda >0$). For any $J$ (full rank or not) and $\lambda >0$, the LM step $p^{\mathrm{LM}}=-(J^{\top }J+\lambda I{)}^{-1}{J}^{\top }r$ is the unique minimiser of $\Vert \left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]p+\left[\begin{array}{c}r\\ 0\end{array}\right]{\Vert }^2$ and satisfies $\nabla f^{\top }{p}^{\mathrm{LM}}<0$ whenever $J^{\top }r\ne 0$.

Proof. The stacked matrix $\tilde{J}=\left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]$ has ${\tilde{J}}^{\top }\tilde{J}=J^{\top }J+\lambda I\succ 0$ for $\lambda >0$ (the $\lambda I$ block makes it positive definite regardless of $\mathrm{rank}J$), so $\tilde{J}$ has full column rank and the least-squares problem has the unique solution given by its normal equations ${\tilde{J}}^{\top }\tilde{J}p=-{\tilde{J}}^{\top }\tilde{r}$, i.e. $(J^{\top }J+\lambda I)p=-J^{\top }r$. For descent, $\nabla f^{\top }{p}^{\mathrm{LM}}=-(J^{\top }r{)}^{\top }(J^{\top }J+\lambda I{)}^{-1}(J^{\top }r)<0$ since $(J^{\top }J+\lambda I{)}^{-1}\succ 0$ and $J^{\top }r\ne 0$. $\square$

Proposition 5 (step-length monotonicity in $\lambda$). $\Vert p^{\mathrm{LM}}(\lambda )\Vert$ is non-increasing on $[0,\infty )$, strictly decreasing when $J^{\top }r\ne 0$, from $\Vert p^{\mathrm{GN}}\Vert$ (or $+\infty$ if $J$ is rank-deficient) at $\lambda =0^{+}$ to $0$ as $\lambda \to \infty$.

Proof. With the SVD $J=U\Sigma V^{\top }$ and $c=V^{\top }{J}^{\top }r$, Exercise 6 gives $\Vert p^{\mathrm{LM}}(\lambda ){\Vert }^2={\sum }_i{c}_i^2/({\sigma }_i^2+\lambda {)}^2$ and $\frac{d}{d\lambda }\Vert p^{\mathrm{LM}}(\lambda ){\Vert }^2=-2{\sum }_i{c}_i^2/({\sigma }_i^2+\lambda {)}^3\le 0$, strict unless all $c_i=0$. The limits follow from the closed form: as $\lambda \to \infty$ each summand vanishes; as $\lambda \to 0^{+}$ a zero singular value with $c_i\ne 0$ sends a term to infinity (rank deficiency), else the sum is $\Vert p^{\mathrm{GN}}{\Vert }^2$. $\square$

Proposition 6 (variable-projection equivalence). For $r(x,y)=\Phi (y)x-d$ with $\Phi (y)$ full column rank, the joint minimiser $(x^{\ast },y^{\ast })$ of $\frac{1}{2}\Vert r(x,y){\Vert }^2$ satisfies $x^{\ast }=\Phi (y^{\ast }{)}^{+}d$ and $y^{\ast }$ minimises $\psi (y)=\frac{1}{2}\Vert {\Pi }_{\Phi (y)}^{\perp }d{\Vert }^2$, ${\Pi }_{\Phi (y)}^{\perp }=I-\Phi (y)\Phi (y{)}^{+}$.

Proof. For fixed $y$ the inner problem ${\min }_x\frac{1}{2}\Vert \Phi (y)x-d{\Vert }^2$ is linear least squares 43.04.01 with unique solution $x(y)=\Phi (y{)}^{+}d$ and optimal value $\frac{1}{2}\Vert d-\Phi (y)\Phi (y{)}^{+}d{\Vert }^2=\frac{1}{2}\Vert {\Pi }_{\Phi (y)}^{\perp }d{\Vert }^2=\psi (y)$, since $\Phi {\Phi }^{+}$ is the orthogonal projector onto $\mathrm{col}\Phi$. Minimising the joint objective is ${\min }_y{\min }_x\frac{1}{2}\Vert \Phi (y)x-d{\Vert }^2={\min }_y\psi (y)$, attained at $y^{\ast }$, with $x^{\ast }=x(y^{\ast })=\Phi (y^{\ast }{)}^{+}d$. Conversely any joint minimiser has its $x$-block optimal for its $y$-block, hence equal to $\Phi (y{)}^{+}d$, so its $y$-block minimises $\psi$. $\square$

Connections Master

• Newton's method for optimization 44.03.03 is the parent against which Gauss-Newton is defined: Gauss-Newton is Newton applied to the structured Hessian $J^{\top }J$ with the residual-curvature term $S=\sum r_j{\nabla }^2{r}_j$ deleted, so the two coincide at a zero-residual solution where $S^{\ast }=0$ and Gauss-Newton inherits Newton's quadratic rate, and diverge in the large-residual regime where $\Vert (J^{\top }J{)}^{-1}{S}^{\ast }\Vert \ge 1$. The integral-remainder identity $g(x_k)-g(x^{\ast })={\int }_0^1{\nabla }^{2}f(x^{\ast }+t{e}_k)e_{k}dt$ proved there is the engine of the convergence-rate proof here, and the line-search globalisation of that unit's predecessor applies unchanged to the Gauss-Newton direction.

• Trust-region methods 44.03.02 supply the exact meaning of Levenberg-Marquardt: the LM step $(J^{\top }J+\lambda I)p=-J^{\top }r$ is the trust-region subproblem solution for the Gauss-Newton model $\Vert Jp+r{\Vert }^2$, the damping $\lambda$ is the trust-region Lagrange multiplier, and the gain-ratio update of $\lambda$ is the radius update of that unit specialised to least squares. The secular-equation monotonicity of $\Vert p(\lambda )\Vert$ established there is the same monotonicity proved here for $\Vert p^{\mathrm{LM}}(\lambda )\Vert$, and because $J^{\top }J⪰0$ the trust-region hard case — central to the general theory there — cannot occur for LM.

• Linear least squares: normal equations, QR, SVD, conditioning 43.04.01 is the inner solver of every Gauss-Newton and LM iteration: each step is a linear least-squares problem ${\min }_p\Vert Jp+r\Vert$ (or its augmented form for LM) solved by QR or SVD rather than by forming $J^{\top }J$, the projector $\Pi =J(J^{\top }J{)}^{-1}{J}^{\top }$ that gives the descent identity is that unit's orthogonal projection onto the column space, and the conditioning warning $\kappa (J^{\top }J)=\kappa (J{)}^2$ is what dictates the QR route. Variable projection reduces a separable problem to the pseudoinverse machinery of that unit applied at every outer step.

• Correlation and regression analysis 26.06.01 is the principal statistical consumer: minimising $\frac{1}{2}\Vert r(x){\Vert }^2$ is maximum-likelihood estimation under independent Gaussian errors, so the linear-regression normal equations of that unit are the single-step special case of Gauss-Newton, the Gauss-Newton matrix $J^{\top }J$ is the Fisher-information / observed-information approximation whose inverse gives the asymptotic parameter covariance, and Levenberg-Marquardt is the standard nonlinear-regression fitting algorithm. The same structure powers bundle adjustment in computer vision and the training of small parametric models, where the residual is the prediction error and $J$ is the sensitivity of predictions to parameters.

Historical & philosophical context Master

The least-squares principle itself originates with Carl Friedrich Gauss, who used it to fit orbital parameters to astronomical observations and published the method in his 1809 Theoria Motus Corporum Coelestium ^{[Gauss 1809]}, with Adrien-Marie Legendre having printed the procedure independently in 1805. The iterative linearisation that bears Gauss's name — solving a sequence of linear least-squares problems for a nonlinear fit — is a direct descendant of his Newton-style root-finding for the normal equations, transferred to the over-determined residual setting.

The damping that stabilises Gauss-Newton for ill-conditioned and large-residual problems was introduced by Kenneth Levenberg in 1944, who proposed adding a multiple of the identity to the normal-equation matrix ^{[Levenberg 1944]}, and refined by Donald Marquardt in 1963, who scaled the damping by the diagonal of $J^{\top }J$ to respect parameter scaling and gave the adaptive update of $\lambda$ that made the method practical ^{[Marquardt 1963]}. The recognition that the Levenberg-Marquardt shift $(J^{\top }J+\lambda I)$ is exactly the trust-region subproblem solution for the Gauss-Newton model came later, in the trust-region literature of the 1970s and 1980s, unifying LM with the radius-controlled methods of Powell, Moré, and Sorensen. The separable-least-squares reduction by elimination of the linear parameters was given its definitive algorithmic form by Gene Golub and Victor Pereyra in 1973, whose variable-projection differentiation formula remains the standard tool ^{[Golub-Pereyra 1973]}. The textbook synthesis organising Gauss-Newton, Levenberg-Marquardt, and the small- versus large-residual analysis around the residual Hessian structure is due to Nocedal and Wright ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Bibliography Master

@book{nocedalwright2006gn,
  author    = {Nocedal, Jorge and Wright, Stephen J.},
  title     = {Numerical Optimization},
  edition   = {2},
  series    = {Springer Series in Operations Research and Financial Engineering},
  publisher = {Springer},
  year      = {2006}
}

@article{levenberg1944,
  author  = {Levenberg, Kenneth},
  title   = {A Method for the Solution of Certain Non-Linear Problems in Least Squares},
  journal = {Quarterly of Applied Mathematics},
  volume  = {2},
  number  = {2},
  year    = {1944},
  pages   = {164--168}
}

@article{marquardt1963gn,
  author  = {Marquardt, Donald W.},
  title   = {An Algorithm for Least-Squares Estimation of Nonlinear Parameters},
  journal = {Journal of the Society for Industrial and Applied Mathematics},
  volume  = {11},
  number  = {2},
  year    = {1963},
  pages   = {431--441}
}

@article{golubpereyra1973,
  author  = {Golub, Gene H. and Pereyra, Victor},
  title   = {The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate},
  journal = {SIAM Journal on Numerical Analysis},
  volume  = {10},
  number  = {2},
  year    = {1973},
  pages   = {413--432}
}

@book{bjorck1996,
  author    = {Bj\"{o}rck, \r{A}ke},
  title     = {Numerical Methods for Least Squares Problems},
  publisher = {SIAM},
  year      = {1996}
}

@book{gauss1809,
  author    = {Gauss, Carl Friedrich},
  title     = {Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium},
  publisher = {Friedrich Perthes and I. H. Besser},
  address   = {Hamburg},
  year      = {1809}
}