Introduction to the Gauss-Newton Method

Fitting a nonlinear model to data, an exponential decay, a growth curve, a sinusoid with an unknown frequency, a circle through noisy sensor readings, leads to an optimization problem that looks deceptively similar to ordinary linear regression: minimize a sum of squared errors. The difference is that the model is no longer linear in its parameters, so the tidy closed-form normal equation no longer applies directly. The Gauss-Newton method is what happens when Newton's method is pointed at exactly this kind of problem: it inherits Newton's fast local convergence, but turns out to need only first derivatives, one derivative order cheaper than a generic Newton step.

The discussion below assumes the derivation of Newton's method for optimization, the quadratic Taylor model, the gradient, and the Hessian, is already familiar; the goal here is to build the specific intuition for what changes, and why it is convenient, when the objective happens to be a sum of squares.

Nonlinear Least Squares

Let \(x\in\mathbb{R}^n\) be the parameter vector to be estimated, for example \(x=(a,b)\) in a two-parameter decay model. Suppose \(m\) observations each produce a scalar mismatch \(r_i(x)\), the residual between what the model predicts and what was actually measured. Stacking these gives the residual vector

\[ r(x) = \begin{pmatrix} r_1(x)\\ \vdots\\ r_m(x)\end{pmatrix} \in\mathbb{R}^m, \qquad r:\mathbb{R}^n\to\mathbb{R}^m. \]

In ordinary linear regression, the residual of each data point is \(r_i(x)=\langle a_i,x\rangle-b_i\), affine in \(x\), which is exactly what makes the least-squares loss there a plain quadratic form with a closed-form minimizer. Here, the residual is allowed to depend on \(x\) nonlinearly, \(r_i(x)=g(x;t_i)-y_i\) for a nonlinear model \(g\), for instance \(g(x;t)=a\,e^{bt}\) with \(x=(a,b)\).

With \(m>n\) observations, or noisy data, there is generally no \(x\) that makes every residual vanish at once. As in linear regression, the overall mismatch is measured by the sum of squared residuals,

\[ \phi(x) := \frac12\|r(x)\|^2 = \frac12\sum_{i=1}^m r_i(x)^2, \]

and the fitting problem becomes an unconstrained minimization, \(\min_{x\in\mathbb{R}^n}\phi(x)\), of exactly the kind Newton's method was built for.

Newton's Method, One Level Removed

Applying Newton's method to \(\phi\) directly would produce the update \(x_{k+1}=x_k-H_\phi(x_k)^{-1}\nabla \phi(x_k)\), the minimizer of the local quadratic model built from the gradient and Hessian of \(\phi\). Computing that Hessian is where the sum-of-squares structure starts to pay off. Differentiating \(\phi(x)=\tfrac12\sum_i r_i(x)^2\) once, with \(J(x)\in\mathbb{R}^{m\times n}\) the Jacobian of \(r\), \(J_{ij}=\partial r_i/\partial x_j\), gives

\[ \frac{\partial\phi}{\partial x_j}(x) = \sum_{i=1}^m r_i(x)\frac{\partial r_i}{\partial x_j}(x) \quad\Longrightarrow\quad \boxed{\nabla\phi(x) = J(x)^\top r(x).} \]

Differentiating once more splits the Hessian of \(\phi\) into two structurally different pieces,

\[ \boxed{H_\phi(x) = J(x)^\top J(x) + \sum_{i=1}^m r_i(x)\,\nabla^2 r_i(x).} \]

The first term, \(J^\top J\), needs only first derivatives of the residuals, the same Jacobian already used for the gradient. The second term needs the full Hessian \(\nabla^2 r_i\) of every residual function, weighted by how large that residual currently is. Exact Newton would need both; this is the expensive half.

The Gauss-Newton Approximation

Gauss-Newton makes one structural simplification: drop the second, expensive term and approximate

\[ H_\phi(x) \approx J(x)^\top J(x). \]

Substituting this approximate Hessian into the Newton system \(H_\phi(x_k)\Delta_k=-\nabla\phi(x_k)\) in place of the true one gives the Gauss-Newton system

\[ \boxed{\bigl(J(x_k)^\top J(x_k)\bigr)\Delta_k = -J(x_k)^\top r(x_k), \qquad x_{k+1}=x_k+\Delta_k.} \]

Why is dropping that term reasonable? Its size is controlled by two things: how curved the residual functions are (\(\nabla^2 r_i\)), and how large the residuals themselves currently are (\(r_i\)). Near a good fit, the residuals are small by definition, so the dropped term shrinks together with them, regardless of how curved the underlying model is; the approximation improves as the iteration converges, exactly when accuracy matters most. It is exactly correct whenever every \(r_i\) is affine in \(x\) (then \(\nabla^2 r_i=0\) identically), which is precisely the linear-regression case. This also explains why Gauss-Newton and its damped relative typically show only Q-linear convergence far from a solution, but accelerate toward the near-quadratic rate of exact Newton once the fit is good, the exact behavior tabulated for Gauss-Newton and Levenberg-Marquardt in the convergence-rate overview of the Newton's method chapter.

A Cleaner Route: Linearize the Model, Not the Objective

There is a second, arguably more intuitive way to arrive at the same system, one that never mentions a Hessian at all. Instead of building a quadratic model of \(\phi\), linearize the residual vector itself around \(x_k\),

\[ r(x_k+\Delta) \approx r(x_k) + J(x_k)\Delta =: r_k + J_k\Delta, \]

and minimize the sum of squares of this linear approximation instead of the true, nonlinear one,

\[ \min_\Delta \; \tfrac12\|r_k+J_k\Delta\|^2. \]

Writing \(\tilde y := -r_k\), this is precisely an ordinary linear least-squares problem, \(\min_\Delta\tfrac12\|J_k\Delta-\tilde y\|^2\), with the Jacobian \(J_k\) playing the role of the design matrix and \(\tilde y\) the role of the target vector in linear regression. Its normal equation,

\[ J_k^\top J_k\,\Delta = J_k^\top\tilde y = -J_k^\top r_k, \]

is exactly the Gauss-Newton system derived above. Every Gauss-Newton step is therefore nothing more exotic than one ordinary linear-regression fit: linearize the model at the current point, solve the resulting linear least-squares problem for the correction \(\Delta_k\), move there, and repeat with a freshly linearized model.

The Single-Residual Case: Gauss-Newton Is Newton's Method

A useful sanity check is what happens with a single parameter and a single residual, \(n=m=1\). The Jacobian is just the scalar derivative \(J(x)=r'(x)\), and the Gauss-Newton system collapses to

\[ \Delta_k = -\frac{J_k r_k}{J_k^2} = -\frac{r(x_k)}{r'(x_k)}, \qquad x_{k+1}=x_k-\frac{r(x_k)}{r'(x_k)}, \]

which is exactly the ordinary Newton-Raphson root-finding update. Taking \(r(x)=x^2-c\), the same residual used to compute square roots by hand, reproduces the familiar iteration

\[ x_{k+1} = x_k - \frac{x_k^2-c}{2x_k}, \]

the identical update derived from scratch in the square-root example of the Newton's method chapter. Gauss-Newton is best understood as the natural generalization of that single-equation Newton step to several residuals and several parameters at once, replacing ordinary division by \(r'(x_k)\) with solving a small linear system.

Worked Example: Fitting an Exponential Decay

Suppose three noiseless measurements come from an exponential decay \(y(t)=2e^{-t}\), sampled at \(t=0,1,2\):

The model is \(g(t;a,b)=a\,e^{bt}\), so the residuals and their Jacobian entries are

\[ r_i(a,b) = a\,e^{bt_i}-y_i, \qquad \frac{\partial r_i}{\partial a}=e^{bt_i}, \qquad \frac{\partial r_i}{\partial b}=a\,t_i\,e^{bt_i}. \]

Starting from a deliberately poor guess \(x_0=(a_0,b_0)=(1,0)\), the model predicts \(g(t;1,0)=1\) at every \(t\), and a single Gauss-Newton step already makes substantial progress:

Four iterations from a flat, entirely wrong starting curve land essentially exactly on the true \((a,b)=(2,-1)\). The residual norm collapses far faster than linearly between the later rows, the near-quadratic behavior expected once the fit is already close to perfect.

The plot shows the three data points together with the model curve at the initial guess (flat, at height \(1\)) and at the converged fit \(a=2,\,b=-1\).

Algorithm

Given an initial guess \(x_0\), repeat for \(k=0,1,2,\dots\) until convergence:

1. Evaluate residuals \(r_k:=r(x_k)\) and Jacobian \(J_k:=J(x_k)\).
2. Solve the linear system \((J_k^\top J_k)\Delta_k = -J_k^\top r_k\) for \(\Delta_k\).
3. Update \(x_{k+1}:=x_k+\Delta_k\).

In practice, \(J_k^\top J_k\) is rarely formed explicitly; the underlying linear least-squares problem \(\min_\Delta\|J_k\Delta+r_k\|^2\) is instead solved via a QR factorization of \(J_k\), or an SVD when \(J_k\) is ill-conditioned, for the same numerical-stability reasons a design matrix in linear regression is factored rather than inverted directly.

Levenberg-Marquardt: Damping the Step

Far from a solution, or whenever \(J_k^\top J_k\) is singular or nearly so, for example with redundant or correlated parameters, the raw Gauss-Newton step can be huge or entirely undefined. The Levenberg-Marquardt method stabilizes it the same way ridge regression stabilizes an ill-conditioned normal equation: add a multiple of the identity before solving,

\[ \bigl(J_k^\top J_k + \lambda I\bigr)\Delta_k = -J_k^\top r_k, \qquad \lambda\ge 0. \]

This is the least-squares counterpart of a damped Newton step: rather than shrinking the step length with an explicit step size \(\sigma_k\in(0,1]\), Levenberg-Marquardt shrinks it implicitly by inflating the approximate curvature. The two limits are informative. As \(\lambda\to0\), the update reduces to plain Gauss-Newton. As \(\lambda\to\infty\), \(\Delta_k\approx-\tfrac1\lambda J_k^\top r_k=-\tfrac1\lambda\nabla\phi(x_k)\), a heavily damped step in the negative-gradient direction, exactly the direction gradient descent takes. Levenberg-Marquardt therefore interpolates continuously between a fast, curvature-aware Gauss-Newton step and a slow but always-safe gradient step, increasing \(\lambda\) when a trial step fails to reduce \(\phi\) and decreasing it again once progress resumes.

When Gauss-Newton Works Well, and When It Doesn't

• It inherits fast, near-quadratic convergence only close to a good fit; far from a solution, or with residuals that stay large, expect the slower Q-linear rate, and prefer damping via Levenberg-Marquardt.
• \(J^\top J\) can be singular or nearly singular whenever parameters are redundant or only weakly identifiable from the data, the same near-singular-Gramian failure mode as in linear regression.
• Strongly curved residual functions combined with genuinely large residuals violate the approximation \(H_\phi\approx J^\top J\) most severely; the dropped term is then not negligible, and the Gauss-Newton direction can even fail to point downhill.
• As with any Newton-type iteration, a poor initial guess can converge to the wrong local minimum, or diverge outright; damping or an explicit line search is the usual remedy.