12.3 The Gauss–Newton method
Gauss–Newton replaces the model by its linear approximation at the current estimate \(\boldsymbol \theta ^{(t)}\) and solves the resulting linear least-squares problem. By Taylor’s theorem, \[ \mathbf {f}(\boldsymbol \theta ) \approx \mathbf {f}(\boldsymbol \theta ^{(t)}) + \mathbf {F}^{(t)}(\boldsymbol \theta - \boldsymbol \theta ^{(t)}) , \] with \(\mathbf {F}^{(t)} = \mathbf {F}(\boldsymbol \theta ^{(t)})\). Writing \(\mathbf {e}^{(t)} = \mathbf {y} - \mathbf {f}(\boldsymbol \theta ^{(t)})\) and \(\boldsymbol \delta = \boldsymbol \theta - \boldsymbol \theta ^{(t)}\), the model becomes approximately \(\mathbf {e}^{(t)} = \mathbf {F}^{(t)}\boldsymbol \delta + \boldsymbol \varepsilon \), a linear regression of the current residuals on the columns of \(\mathbf {F}^{(t)}\). Its least-squares solution (Theorem 3.3.1) gives the step \[ \boldsymbol \delta ^{(t)} = \bigl (\mathbf {F}^{(t)\top }\mathbf {F}^{(t)}\bigr )^{-1}\mathbf {F}^{(t)\top }\mathbf {e}^{(t)} , \qquad \boldsymbol \theta ^{(t+1)} = \boldsymbol \theta ^{(t)} + \boldsymbol \delta ^{(t)} . \] The iteration stops when the step, or the relative change in \(S\), is negligible. At a fixed point \(\boldsymbol \delta = \mathbf {0}\), so \(\mathbf {F}^\top \mathbf {e} = \mathbf {0}\): the normal equations hold.
Theorem 12.3.1. Suppose that \(\mathbf {F}^{(t)}\) has full column rank and that \(\mathbf {F}^{(t)\top }\mathbf {e}^{(t)} \ne \mathbf {0}\). Then the Gauss–Newton step \(\boldsymbol \delta ^{(t)}\) is a descent direction: \(S(\boldsymbol \theta ^{(t)} + \lambda \boldsymbol \delta ^{(t)}) < S(\boldsymbol \theta ^{(t)})\) for all small enough \(\lambda > 0\).
Proof. Write \(\mathbf {F}\), \(\mathbf {e}\), \(\boldsymbol \delta \) for the quantities at step \(t\). By Theorem 12.2.1, the gradient of \(S\) is \(\mathbf {g} = -2\mathbf {F}^\top \mathbf {e}\), so the derivative of \(S(\boldsymbol \theta ^{(t)} + \lambda \boldsymbol \delta )\) at \(\lambda = 0\) is \[ \mathbf {g}^\top \boldsymbol \delta = -2\,\mathbf {e}^\top \mathbf {F}(\mathbf {F}^\top \mathbf {F})^{-1}\mathbf {F}^\top \mathbf {e} = -2\,\mathbf {u}^\top (\mathbf {F}^\top \mathbf {F})^{-1}\mathbf {u} , \qquad \mathbf {u} = \mathbf {F}^\top \mathbf {e} \ne \mathbf {0} . \] \((\mathbf {F}^\top \mathbf {F})^{-1}\) is positive definite, so this is negative, and \(S\) decreases for small \(\lambda > 0\). □
A full step can still overshoot when the linear approximation is poor, far from the solution. Two remedies are standard.
- Step halving: if \(S(\boldsymbol \theta ^{(t)} + \boldsymbol \delta ) \ge S(\boldsymbol \theta ^{(t)})\), try \(\lambda = \frac 12, \frac 14, \dots \) until \(S\) decreases. By the theorem, this succeeds after finitely many halvings.
- The Levenberg–Marquardt method takes the step \(\boldsymbol \delta = (\mathbf {F}^\top \mathbf {F} + k\mathbf {I})^{-1}\mathbf {F}^\top \mathbf {e}\) with \(k \ge 0\). This is the ridge estimator of Theorem 8.4.2 applied to the linearised problem. With \(k = 0\) it is the Gauss–Newton step; as \(k\) grows, \(k\boldsymbol \delta \to \mathbf {F}^\top \mathbf {e} = -\frac 12\mathbf {g}\), a short step straight downhill (steepest descent). The value of \(k\) is reduced after a successful step and increased after an unsuccessful one, so the method behaves like steepest descent far from the minimum, where that is safe, and like Gauss–Newton near it, where that is fast.
Starting values
Iterative methods need a starting value, and a poor one can lead to slow convergence, to a local minimum, or to divergence. Good sources are:
- the meaning of the parameters: an asymptote, a value at \(x = 0\), a half-life, read from a plot of the data;
- the fit of a linearised version of the model, when there is one;
- a grid of trial values, keeping the one with the smallest \(S\).
It is good practice to rerun the fit from several starting values and check that it reaches the same estimate.
Example 12.3.2. A drug was injected, and its concentration in the blood, in \(\mu \)g/ml, was then measured every hour:
| Hours, \(x\) | 1 | 2 | 3 | 4 | 5 | 6 |
| \(y\) | 15.0 | 10.8 | 8.3 | 5.8 | 4.6 | 3.2 |
Fit \(y = \theta _1e^{-\theta _2x} + \varepsilon \) by Gauss–Newton, first from the deliberately poor starting value \(\boldsymbol \theta ^{(0)} = (10, 0.1)\), and then from the log-linear fit.
Solution. The first step. The columns of \(\mathbf {F}\) are \(\partial f/\partial \theta _1 = e^{-\theta _2x}\) and \(\partial f/\partial \theta _2 = -\theta _1xe^{-\theta _2x}\). At \((10, 0.1)\), together with the residuals \(e = y - 10e^{-0.1x}\), they are
| \(x\) | \(e^{-0.1x}\) | \(-10xe^{-0.1x}\) | \(e\) |
| 1 | 0.9048 | \(-9.048\) | 5.952 |
| 2 | 0.8187 | \(-16.375\) | 2.613 |
| 3 | 0.7408 | \(-22.225\) | 0.892 |
| 4 | 0.6703 | \(-26.813\) | \(-0.903\) |
| 5 | 0.6065 | \(-30.327\) | \(-1.465\) |
| 6 | 0.5488 | \(-32.929\) | \(-2.288\) |
so that \[ \mathbf {F}^\top \mathbf {F} = \begin {pmatrix} 3.156 & -92.50 \\ -92.50 & 3566.9 \end {pmatrix} , \qquad \mathbf {F}^\top \mathbf {e} = \begin {pmatrix} 5.435 \\ 27.545 \end {pmatrix} . \] The determinant is \(3.156(3566.9) - 92.50^2 = 2702\), and \begin {align*} \delta _1 &= \frac {3566.9(5.435) + 92.50(27.545)}{2702} = 8.117 , \\ \delta _2 &= \frac {92.50(5.435) + 3.156(27.545)}{2702} = 0.218 , \end {align*}
so \(\boldsymbol \theta ^{(1)} = (18.117, 0.318)\), and \(S\) falls from \(51.24\) to \(8.17\).
The iterations. Continuing,
| Step | \(\theta _1\) | \(\theta _2\) | \(S(\boldsymbol \theta )\) |
| 0 | 10.0000 | 0.1000 | 51.242 |
| 1 | 18.1169 | 0.3182 | 8.165 |
| 2 | 20.2301 | 0.3026 | 0.1645 |
| 3 | 20.2380 | 0.3047 | 0.1539 |
| 4 | 20.2386 | 0.3047 | 0.1539 |
so \(\hat {\boldsymbol \theta } = (20.239, 0.3047)\) and \(S(\hat {\boldsymbol \theta }) = 0.1539\). The fitted values are \(14.92\), \(11.00\), \(8.11\), \(5.98\), \(4.41\), \(3.25\), and the residuals \(0.08\), \(-0.20\), \(0.19\), \(-0.18\), \(0.19\), \(-0.05\) show no pattern.
From the linearised fit. Regressing \(\log y\) on \(x\) gives \(\log y = 3.006 - 0.3041x\), so \[ \boldsymbol \theta ^{(0)} = (e^{3.006}, 0.3041) = (20.210, 0.3041) , \] with \(S = 0.1542\) already. One Gauss–Newton step gives \((20.2385, 0.3047)\), and the second agrees with the estimate above to four figures. The log-linear fit is close here because the data are nearly noise-free; it is a starting value, not the answer, since it minimises a different criterion.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.