Prev
Unit 5 · Least Squares and Data Fitting
Nonlinear Least Squares
65m100 XP
Next
📐Least Squares and Data Fitting

Nonlinear Least Squares

Gauss–Newton, the damped step, Levenberg–Marquardt, GNSS point positioning and residual analysis

Philadelphia 1944, Wilmington 1963
Frankford Arsenal, Philadelphia, in the middle of the war. Kenneth Levenberg, a mathematician at the Army's arsenal, faces models whose parameters enter nonlinearly. The standard method — linearise, solve the linear least-squares problem, repeat — often wanders off or diverges from a poor starting guess. In a five-page note in the Quarterly of Applied Mathematics, A Method for the Solution of Certain Non-Linear Problems in Least Squares, he suggests damping each step: penalise long steps, so that far from the answer the method creeps cautiously downhill and near the answer it takes full linearised steps. The note is little read. Nineteen years later, Donald Marquardt, a statistician at the DuPont Company in Wilmington fitting kinetic models to chemical data, arrives at the same idea, adds a scaling that makes it work for parameters of very different sizes, and publishes it in the SIAM Journal in 1963. He distributes a FORTRAN program. It spreads through every laboratory that fits data. Today the Levenberg–Marquardt algorithm is inside scipy.optimize.curve_fit, inside camera calibration, bundle adjustment, orbit determination and the navigation filter in your phone.
A model y=f(t;p)y = f(t;\mathbf p) with parameters p\mathbf p is fitted to data (ti,yi)(t_i, y_i) by minimising
S(p)=12∑iri(p)2=12∥r(p)∥2,ri=yi−f(ti;p).S(\mathbf p) = \tfrac12\sum_i r_i(\mathbf p)^2 = \tfrac12\|\mathbf r(\mathbf p)\|^2, \qquad r_i = y_i - f(t_i;\mathbf p).
When ff is nonlinear in p\mathbf p — an unknown decay rate e−kte^{-kt}, an unknown frequency, a position inside a distance — SS is no longer a quadratic bowl. There is no formula; we iterate, and every iteration is a linear least-squares problem from the earlier lessons.
💡Gauss–Newton
Linearise the residuals about the current pk\mathbf p_k with the Jacobian Jij=∂f(ti)/∂pjJ_{ij} = \partial f(t_i)/\partial p_j: r(pk+δ)≈rk−Jkδ.\mathbf r(\mathbf p_k + \boldsymbol\delta) \approx \mathbf r_k - J_k\boldsymbol\delta. Choose δ\boldsymbol\delta to minimise ∥rk−Jkδ∥\|\mathbf r_k - J_k\boldsymbol\delta\| — a linear least-squares problem, solved by QR — and set pk+1=pk+δ\mathbf p_{k+1} = \mathbf p_k + \boldsymbol\delta. It is Newton's method (lesson 2-4) applied to the gradient of SS, with the second-derivative terms dropped. When the residuals at the solution are small, those terms are small too, and Gauss–Newton converges almost quadratically.
Derivation
Gauss–Newton is Newton with one term dropped
Newton's method for a minimum (lesson 2-4, applied to ∇S=0\nabla S = \mathbf 0) needs the gradient and the second-derivative matrix of S=12∑iri2S = \tfrac12\sum_i r_i^2. Since ri=yi−f(ti;p)r_i = y_i - f(t_i; \mathbf p), its derivatives are minus those of ff, so ∇ri=−(row i of J)\nabla r_i = -(\text{row } i \text{ of } J). By the chain rule:
∇S=∑iri ∇ri=−JTr,∇2S=∑i∇ri ∇riT+∑iri ∇2ri=JTJ−∑iri ∇2fi.\nabla S = \sum_i r_i\,\nabla r_i = -J^{\mathsf T}\mathbf r, \qquad \nabla^2S = \sum_i \nabla r_i\,\nabla r_i^{\mathsf T} + \sum_i r_i\,\nabla^2 r_i = J^{\mathsf T}J - \sum_i r_i\,\nabla^2 f_i.
Newton's step solves ∇2S δ=−∇S\nabla^2S\,\boldsymbol\delta = -\nabla S:
(JTJ−∑iri ∇2fi)δ=JTr.\Big(J^{\mathsf T}J - \sum_i r_i\,\nabla^2 f_i\Big)\boldsymbol\delta = J^{\mathsf T}\mathbf r.
The second-derivative sum is the expensive part, and it is multiplied by the residuals. Drop it, and what remains, JTJδ=JTrJ^{\mathsf T}J\boldsymbol\delta = J^{\mathsf T}\mathbf r, is exactly the normal equations of the linear problem min⁡∥r−Jδ∥\min\|\mathbf r - J\boldsymbol\delta\| — the Gauss–Newton step. When the model fits well the residuals are small, the dropped term is small, and the step is nearly Newton's. Only first derivatives are ever needed.
📐A cooling casting
A casting cools in air: T(t)=Ta+(T0−Ta)e−ktT(t) = T_a + (T_0 - T_a)e^{-kt}, with ambient temperature TaT_a, initial temperature T0T_0 and rate kk all unknown. Twelve noisy readings over 50 minutes. From a rough guess (2020 °C, 300300 °C, 0.10.1/min), Gauss–Newton's step sizes go 77, 5.3, 2.7, 8×10−4, 3×10−6, 10−8, 6×10−1177,\ 5.3,\ 2.7,\ 8\times10^{-4},\ 3\times10^{-6},\ 10^{-8},\ 6\times10^{-11}, landing on Ta=21.14T_a = 21.14 °C, T0=378.41T_0 = 378.41 °C, k=0.0701k = 0.0701/min — a time constant of 14.3 minutes — in agreement with scipy.optimize.curve_fit.
💡Levenberg–Marquardt
Far from the answer the linearisation is poor and a full Gauss–Newton step can overshoot or diverge. Levenberg's fix adds a damping term: (JTJ+μD) δ=JTr,D=I (Levenberg) or diag⁡(JTJ) (Marquardt).(J^{\mathsf T}J + \mu D)\,\boldsymbol\delta = J^{\mathsf T}\mathbf r, \qquad D = I \text{ (Levenberg) or } \operatorname{diag}(J^{\mathsf T}J) \text{ (Marquardt)}. Large μ\mu: a short step along the steepest-descent direction — slow but safe. Small μ\mu: the Gauss–Newton step — fast near the answer. The algorithm adapts μ\mu automatically: if a step reduces SS, accept it and decrease μ\mu; if not, reject it and increase μ\mu. It is a trust-region method in disguise: μ\mu controls how far we trust the linear model.
Derivation
Where the damped equations come from
Penalise long steps: minimise the linearised misfit plus a cost on the step itself,
∥r−Jδ∥2+μ δTD δ.\|\mathbf r - J\boldsymbol\delta\|^2 + \mu\,\boldsymbol\delta^{\mathsf T}D\,\boldsymbol\delta.
Setting the gradient with respect to δ\boldsymbol\delta to zero (as in lesson 5-1) gives −2JT(r−Jδ)+2μDδ=0-2J^{\mathsf T}(\mathbf r - J\boldsymbol\delta) + 2\mu D\boldsymbol\delta = \mathbf 0, which is
(JTJ+μD) δ=JTr.(J^{\mathsf T}J + \mu D)\,\boldsymbol\delta = J^{\mathsf T}\mathbf r.
The two limits explain the behaviour. As μ→0\mu \to 0 it is the Gauss–Newton step. For large μ\mu (with D=ID = I) the μD\mu D term dominates and δ≈JTr/μ=−∇S/μ\boldsymbol\delta \approx J^{\mathsf T}\mathbf r/\mu = -\nabla S/\mu: a short step straight downhill. Marquardt's D=diag⁡(JTJ)D = \operatorname{diag}(J^{\mathsf T}J) scales the penalty to each parameter's own sensitivity, so a decay rate of 0.07 per minute and a temperature of 380 °C are damped equally fairly.
📐A damped vibration
An accelerometer records a decaying oscillation y=Ae−ζωtcos⁡(ωdt+ϕ)y = Ae^{-\zeta\omega t}\cos(\omega_d t + \phi); fit AA, ζ\zeta, ω\omega, ϕ\phi. Start with the frequency guessed 12% low. Plain Gauss–Newton fails at its second iteration: its step drives the damping ratio negative. Levenberg–Marquardt (scipy.optimize.least_squares(method='lm')) converges in 10 function evaluations to A=2.004A = 2.004, ζ=0.0503\zeta = 0.0503, f=1.6996f = 1.6996 Hz, ϕ=0.401\phi = 0.401. The residuals have mean −0.002-0.002 and standard deviation 0.0520.052 — matching the 0.05 noise that was added, which is what a correct model should leave behind.
the theory of probabilities is basically just common sense reduced to calculus
– Pierre-Simon Laplace, Essai philosophique sur les probabilités (1814)
💡GNSS point positioning
A receiver measures pseudoranges to m≥4m \ge 4 satellites at known positions si\mathbf s_i: ρi=∥si−r∥+c δt+εi,\rho_i = \|\mathbf s_i - \mathbf r\| + c\,\delta t + \varepsilon_i, with unknown position r\mathbf r and receiver clock error c δtc\,\delta t (the cheap quartz clock in a receiver is off by milliseconds — hundreds of kilometres of range). Four unknowns, nonlinear in r\mathbf r. The Jacobian rows are [−eiT, 1][-\mathbf e_i^{\mathsf T},\ 1], where ei\mathbf e_i is the unit vector from the receiver to satellite ii: the geometry matrix GG. Started at the centre of the Earth with zero clock error, Gauss–Newton on six satellites moves the estimate 7,674 km, then 1,253 km, 42 km, 48 m, then nothing: five iterations to a position 4.4 m from the truth with 3 m range noise. (GTG)−1(G^{\mathsf T}G)^{-1} gives the dilution of precision: GDOP 2.8, PDOP 2.4 for this sky. Satellites bunched in one part of the sky make GG ill-conditioned and the DOP — the error amplification of lesson 3-2 — large.
Derivation
The geometry matrix and dilution of precision
The Jacobian row. The predicted pseudorange is ρi(r,cδt)=∥si−r∥+c δt\rho_i(\mathbf r, c\delta t) = \|\mathbf s_i - \mathbf r\| + c\,\delta t. Differentiate the distance with respect to the receiver position: ∥si−r∥=(si−r)T(si−r)\|\mathbf s_i - \mathbf r\| = \sqrt{(\mathbf s_i - \mathbf r)^{\mathsf T}(\mathbf s_i - \mathbf r)}, so
∂ρi∂r=−(si−r)T∥si−r∥=−eiT,∂ρi∂(c δt)=1.\frac{\partial\rho_i}{\partial\mathbf r} = \frac{-(\mathbf s_i - \mathbf r)^{\mathsf T}}{\|\mathbf s_i - \mathbf r\|} = -\mathbf e_i^{\mathsf T}, \qquad \frac{\partial\rho_i}{\partial(c\,\delta t)} = 1.
Moving the receiver towards a satellite shortens that range one-for-one along the line of sight; the clock adds to every range equally. Row ii of the Jacobian is [−eiT, 1][-\mathbf e_i^{\mathsf T},\ 1].
Dilution of precision. Near the solution each Gauss–Newton step is a linear least-squares problem with matrix GG. If every range has independent errors of standard deviation σ\sigma, lesson 5-1's covariance formula with W=I/σ2W = I/\sigma^2 gives
Cov⁡(x^,y^,z^,cδt^)=σ2 (GTG)−1.\operatorname{Cov}(\hat x, \hat y, \hat z, c\hat{\delta t}) = \sigma^2\,(G^{\mathsf T}G)^{-1}.
The diagonal holds the variances of the four estimates. Their total, divided by σ2\sigma^2, is a pure number fixed by the geometry alone: GDOP=trace⁡(GTG)−1\text{GDOP} = \sqrt{\operatorname{trace}(G^{\mathsf T}G)^{-1}}, and PDOP is the same over the three position entries. A 3 m range error with PDOP 2.4 means about 3×2.4≈73 \times 2.4 \approx 7 m of position uncertainty.
⚠Local minima and starting values
Nonlinear least-squares problems can have several local minima; LM finds one near where it starts. A frequency guessed 50% wrong may lock onto a harmonic. Use physics for starting values (here: count zero crossings for ω\omega, read AA off the first peak), fit in stages (linear parameters first), or try several starts. For GNSS the problem is benign — almost any start converges — which is why receivers can start from the Earth's centre.
💡Where you will use this
GNSS Navigation: single-point positioning exactly as above; precise point positioning and RTK add carrier phases and more unknowns. Orbital Mechanics: orbit determination from range, Doppler and angle measurements — Gauss's Ceres problem, solved by batch Gauss–Newton. Attitude Determination and Camera calibration: star-tracker and camera models fitted by LM. Materials and Geotechnics: fitting creep laws, consolidation curves and constitutive models to laboratory tests.
📚Further reading
  • Jorge Nocedal and Stephen Wright, Numerical Optimization (2nd ed., Springer, 2006), chapter 10, 'Least-Squares Problems' — Gauss–Newton and Levenberg–Marquardt as trust-region methods.
  • Donald Marquardt, 'An algorithm for least-squares estimation of nonlinear parameters', Journal of the Society for Industrial and Applied Mathematics 11(2), 431–441 (1963).
  • Pratap Misra and Per Enge, Global Positioning System: Signals, Measurements, and Performance (2nd ed., Ganga-Jamuna Press, 2006) — the pseudorange fix and dilution of precision as GNSS engineers use them.
↓But there is a problem
This unit was about solving Ax≈bA\mathbf x \approx \mathbf b for x\mathbf x. Along the way we kept meeting special directions: the singular vectors, the columns of QQ, the axes of the energy ellipse. Some matrices have directions they do not rotate at all — they only stretch them. On 7 November 1940, the main span of the Tacoma Narrows Bridge, less than five months old, began to twist in a 64 km/h wind. Within an hour the deck had torn itself apart. The motion that destroyed it was one particular shape of vibration, at one particular frequency — an eigenvector and an eigenvalue of the bridge.