Prev
Unit 5 · Least Squares and Data Fitting
Linear Least Squares and the Normal Equations
60m100 XP
Next
📐Least Squares and Data Fitting

Linear Least Squares and the Normal Equations

Overdetermined systems, minimising ‖Ax − b‖, projection, normal equations and weighting

Palermo, 1 January 1801
On the first night of the nineteenth century, Giuseppe Piazzi, director of the Palermo Observatory, notices a faint star in Taurus that is not in his catalogue. The next night it has moved. He follows it for 41 days and records 24 positions — then it disappears into the glare of the Sun. By autumn it should be visible again, but nobody knows where. Piazzi's observations cover a tiny arc of a long orbit; each is slightly wrong; and the orbit has six unknown elements. The astronomers of Europe look and do not find it. In Brunswick, Carl Friedrich Gauss, 24 years old, takes up the problem. He develops a new method for determining an orbit from a few observations, fits it to Piazzi's data, and sends his prediction to Franz Xaver von Zach. On 31 December 1801, von Zach and Heinrich Olbers find the object close to where Gauss said it would be. It is Ceres, the largest body in the asteroid belt. When Gauss finally explained his methods in 1809, at their heart was a principle he said he had used since 1795: when the data cannot all be satisfied, choose the answer that makes the sum of the squares of the misfits as small as possible.
A measurement model that is linear in its unknowns gives a system Ax=bA\mathbf x = \mathbf b with AA of size m×nm \times n. With m>nm > n — more measurements than unknowns, which is what every careful experiment does — the system is overdetermined, and with real data it has no exact solution. The residual r=b−Ax\mathbf r = \mathbf b - A\mathbf x cannot be made zero. The least-squares solution makes it as short as possible:
min⁡x ∥b−Ax∥22=min⁡x ∑i=1m(bi−∑jaijxj)2\min_{\mathbf x}\ \|\mathbf b - A\mathbf x\|_2^2 = \min_{\mathbf x}\ \sum_{i=1}^{m}\Big(b_i - \sum_j a_{ij}x_j\Big)^2

Squares rather than absolute values: the problem becomes smooth, has a unique solution when A has full column rank, and is solved by linear algebra alone. (With independent Gaussian errors it is also the maximum-likelihood estimate.)

💡Geometry: least squares is projection
As x\mathbf x varies, AxA\mathbf x sweeps out the column space of AA — an nn-dimensional plane inside Rm\mathbb R^m. The data vector b\mathbf b is generally not on it. The closest point of the plane to b\mathbf b is its orthogonal projection Ax^A\hat{\mathbf x}, and the residual b−Ax^\mathbf b - A\hat{\mathbf x} is perpendicular to the plane — perpendicular to every column of AA: AT(b−Ax^)=0⟺ATA x^=ATb.A^{\mathsf T}(\mathbf b - A\hat{\mathbf x}) = \mathbf 0 \quad\Longleftrightarrow\quad A^{\mathsf T}A\,\hat{\mathbf x} = A^{\mathsf T}\mathbf b. These are the normal equations — 'normal' in the sense of perpendicular. ATAA^{\mathsf T}A is n×nn \times n, symmetric, and positive definite when the columns of AA are independent, so Cholesky (lesson 3-3) solves them.
Derivation
The same equations by calculus — and the straight-line fit
Multiply out the squared length of the residual:
S(x)=∥b−Ax∥2=bTb−2 xTATb+xTATA x.S(\mathbf x) = \|\mathbf b - A\mathbf x\|^2 = \mathbf b^{\mathsf T}\mathbf b - 2\,\mathbf x^{\mathsf T}A^{\mathsf T}\mathbf b + \mathbf x^{\mathsf T}A^{\mathsf T}A\,\mathbf x.
At a minimum every partial derivative is zero. Differentiating the linear term gives −2ATb-2A^{\mathsf T}\mathbf b and the quadratic term 2ATAx2A^{\mathsf T}A\mathbf x (the same rule as d(ax2)/dx=2axd(ax^2)/dx = 2ax, since ATAA^{\mathsf T}A is symmetric):
∇S=−2ATb+2ATA x=0⟺ATA x=ATb.\nabla S = -2A^{\mathsf T}\mathbf b + 2A^{\mathsf T}A\,\mathbf x = \mathbf 0 \qquad\Longleftrightarrow\qquad A^{\mathsf T}A\,\mathbf x = A^{\mathsf T}\mathbf b.
The geometry and the calculus agree. The second derivative, 2ATA2A^{\mathsf T}A, is positive definite when the columns are independent, so this stationary point is the minimum.
A straight line. For y≈c0+c1ty \approx c_0 + c_1t through mm points, the rows of AA are [1, ti][1,\ t_i], and the normal equations are only 2×22 \times 2:
[m∑ti∑ti∑ti2][c0c1]=[∑yi∑tiyi]⟹c1=m∑tiyi−∑ti∑yim∑ti2−(∑ti)2,\begin{bmatrix} m & \sum t_i \\ \sum t_i & \sum t_i^2 \end{bmatrix}\begin{bmatrix} c_0 \\ c_1 \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum t_iy_i \end{bmatrix} \quad\Longrightarrow\quad c_1 = \frac{m\sum t_iy_i - \sum t_i\sum y_i}{m\sum t_i^2 - \big(\sum t_i\big)^2},
and c0=(∑yi−c1∑ti)/mc_0 = (\sum y_i - c_1\sum t_i)/m — the formulas behind every calculator's 'linear regression' key, and the tensile-test fit below.
📐Young's modulus from a tensile test
A steel bar, 500 mm long with a 100 mm² cross-section, is loaded from 0 to 35 kN in 5 kN steps; the extensometer reads 0.020, 0.150, 0.271, 0.383, 0.506, 0.627, 0.755 and 0.873 mm. The model is extension =c0+c1⋅= c_0 + c_1 \cdot load, with c0c_0 absorbing the seating of the grips. AA has rows [1, loadi][1,\ \text{load}_i]: eight equations, two unknowns. The normal equations give c0=0.0237c_0 = 0.0237 mm and c1=0.024255c_1 = 0.024255 mm/kN, with a residual norm of 0.0106 mm, and ATr=0A^{\mathsf T}\mathbf r = \mathbf 0 to rounding. Then E=(F/A)/(ΔL/L)=206E = (F/A)/(\Delta L/L) = 206 GPa — a typical structural steel.
When the number of unknown quantities is equal to the number of the observed quantities depending on them, the former may be so determined as exactly to satisfy the latter. But when the number of the former is less than that of the latter, an absolutely exact agreement cannot be obtained, unless the observations possess absolute accuracy.
– Carl Friedrich Gauss, Theoria Motus Corporum Coelestium (1809), tr. C. H. Davis
Weighted least squares. Measurements are rarely equally good. If measurement ii has standard deviation σi\sigma_i, divide its equation by σi\sigma_i before solving — equivalently, minimise ∑(ri/σi)2\sum (r_i/\sigma_i)^2. Precise readings then count more:
ATWA x^=ATWb,W=diag⁡(1/σi2),Cov⁡(x^)=(ATWA)−1.A^{\mathsf T}WA\,\hat{\mathbf x} = A^{\mathsf T}W\mathbf b, \qquad W = \operatorname{diag}(1/\sigma_i^2), \qquad \operatorname{Cov}(\hat{\mathbf x}) = (A^{\mathsf T}WA)^{-1}.
A baseline measured six times by tape (σ=5\sigma = 5 mm), three times by total station (2 mm) and twice by GNSS (10 mm) has a weighted estimate of 125.43335 m ±\pm 1.00 mm — better than the total station alone (±\pm 1.15 mm), because the other instruments still contribute. The covariance matrix is as important as the estimate: it is what a surveyor or a navigation filter carries forward.
Derivation
Weights, and the covariance of the estimate
The language of measurement error. Write each reading as bi=(true value)+εib_i = (\text{true value}) + \varepsilon_i. The error εi\varepsilon_i is random: on average zero, with variance σi2=E[εi2]\sigma_i^2 = E[\varepsilon_i^2], the mean of its square (σi\sigma_i is its standard deviation). Errors of different instruments are independent: E[εiεj]=0E[\varepsilon_i\varepsilon_j] = 0 for i≠ji \ne j. For a vector of errors these averages are collected in the covariance matrix E[εεT]E[\boldsymbol\varepsilon\boldsymbol\varepsilon^{\mathsf T}]: variances on the diagonal, the averages of products of different errors off it. Here it is diag⁡(σi2)=W−1\operatorname{diag}(\sigma_i^2) = W^{-1}. (Unit 12 develops all of this further.)
The weighted fit. Minimising ∑(ri/σi)2=(b−Ax)TW(b−Ax)\sum(r_i/\sigma_i)^2 = (\mathbf b - A\mathbf x)^{\mathsf T}W(\mathbf b - A\mathbf x) is ordinary least squares on rows divided by σi\sigma_i, so the calculus above gives ATWA x^=ATWbA^{\mathsf T}WA\,\hat{\mathbf x} = A^{\mathsf T}W\mathbf b, that is x^=Kb\hat{\mathbf x} = K\mathbf b with K=(ATWA)−1ATWK = (A^{\mathsf T}WA)^{-1}A^{\mathsf T}W.
How errors pass through a linear formula. If y=Kb\mathbf y = K\mathbf b, the error of y\mathbf y is KεK\boldsymbol\varepsilon, and its covariance is
E[(Kε)(Kε)T]=K E[εεT] KT.E\big[(K\boldsymbol\varepsilon)(K\boldsymbol\varepsilon)^{\mathsf T}\big] = K\,E[\boldsymbol\varepsilon\boldsymbol\varepsilon^{\mathsf T}]\,K^{\mathsf T}.
For the weighted estimate, with E[εεT]=W−1E[\boldsymbol\varepsilon\boldsymbol\varepsilon^{\mathsf T}] = W^{-1} and N=ATWAN = A^{\mathsf T}WA:
Cov⁡(x^)=N−1ATW W−1 WAN−1=N−1(ATWA)N−1=(ATWA)−1.\operatorname{Cov}(\hat{\mathbf x}) = N^{-1}A^{\mathsf T}W\,W^{-1}\,WAN^{-1} = N^{-1}(A^{\mathsf T}WA)N^{-1} = (A^{\mathsf T}WA)^{-1}.
The baseline. Every reading measures the same length, so AA is a column of ones and ATWA=∑1/σi2=6/25+3/4+2/100=1.01A^{\mathsf T}WA = \sum 1/\sigma_i^2 = 6/25 + 3/4 + 2/100 = 1.01 per mm². The standard deviation of the estimate is 1/1.01=1.001/\sqrt{1.01} = 1.00 mm; the three total-station readings alone give 1/3/4=1.151/\sqrt{3/4} = 1.15 mm. Precisions add as 1/σ21/\sigma^2, which is why even the poor instruments help a little.
⚠The normal equations square the condition number
κ(ATA)=κ(A)2\kappa(A^{\mathsf T}A) = \kappa(A)^2 (lesson 5-3 proves it with the singular value decomposition). Every digit κ(A)\kappa(A) costs, the normal equations cost twice. Fitting a degree-9 polynomial in the monomial basis on 50 points of [0,1][0, 1]: κ(A)=3.6×106\kappa(A) = 3.6 \times 10^6 but κ(ATA)=1.3×1013\kappa(A^{\mathsf T}A) = 1.3 \times 10^{13}, and the coefficients come out with relative error 2×10−42 \times 10^{-4} from the normal equations against 5×10−115 \times 10^{-11} from an orthogonal-factorisation solver. For small, well-conditioned problems the normal equations are fine and fast. For anything else, use QR or the SVD.
🕰Gauss, Legendre and priority
Adrien-Marie Legendre published the method of least squares first, in an 1805 appendix on comet orbits, and named it. Gauss's 1809 claim to have used it since 1795 started one of the famous priority disputes of mathematics. Gauss's lasting contribution was the justification: in 1809 he connected least squares to the normal distribution of errors, and in 1821–23 he proved that it gives the minimum-variance linear unbiased estimate whatever the error distribution — what is now the Gauss–Markov theorem.
💡Where you will use this
GNSS Navigation: every position fix is a (weighted, then nonlinear) least-squares problem, and orbit determination of the satellites themselves is Gauss's Ceres problem with better data. Attitude Determination: Wahba's problem — the rotation that best maps measured star and Sun directions onto catalogue directions — is a weighted least-squares problem. Materials and Structures: calibration, modulus and creep-law fits. Surveying and Geotechnics: network adjustment, settlement trend fitting.
📚Further reading
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 3, 'Linear Least Squares'.
  • Åke Björck, Numerical Methods for Least Squares Problems (SIAM, 1996) — the standard reference, from normal equations to weighted and constrained problems.
  • Stephen Stigler, 'Gauss and the invention of least squares', Annals of Statistics 9(3), 465–474 (1981) — the Gauss–Legendre priority question, examined from the sources.
↓But there is a problem
The normal equations are elegant: multiply by ATA^{\mathsf T}, and an m×nm \times n problem becomes an n×nn \times n SPD system. But that multiplication squares the condition number, and an ill-conditioned fitting problem — a high-degree polynomial, nearly collinear regressors, a GNSS fix with satellites bunched together — loses twice the digits it should. In 1965 Gene Golub, then a young researcher at Stanford, published a short paper showing how to solve least-squares problems without ever forming ATAA^{\mathsf T}A — using reflections.