Linear Least Squares and the Normal Equations

Derivation

Linear Least Squares and the Normal Equations

D-121 Home PU-106 Threads chance Depends on Least Squares from Gaussian Likelihood, matrix-inverse-and-transpose

Statement

For a model that is linear in its parameters, y = Ap + noise, the weighted least-squares estimate p̂ that minimises χ² = (y − Ap)^TW(y − Ap) satisfies the normal equations (A^TWA)p̂ = A^TWy, where A is the N×M design matrix, y is the vector of N observations, and W is a symmetric positive-definite weight matrix.

Why it matters

Fitting a straight line, a polynomial, a spectral baseline, or a sum of basis functions to data are all the same computation once written in matrix form. The normal equations turn "find the best curve" into "solve one linear system," making the whole of linear regression a single closed-form result rather than an iterative search.

The weight matrix W is where the physics of measurement enters. Taking W = C⁻¹, the inverse of the data covariance, is exactly what maximising a Gaussian likelihood demands, so the estimator is not an ad-hoc recipe but the statistically optimal (minimum-variance unbiased) linear estimator. It also delivers the parameter uncertainties for free as (A^TWA)⁻¹.

Assumptions

The model is linear in the parameters p.If y depends non-linearly on p (e.g. an exponential decay rate), no fixed matrix A exists; the objective is no longer quadratic and has no closed-form minimiser — you must iterate (Gauss–Newton / Levenberg–Marquardt). The weight matrix W is symmetric and positive-definite.If W is not symmetric, the cross-term identity used to combine the two linear pieces fails and the gradient picks up an extra antisymmetric part. If it is not positive-definite, χ² is not bounded below and the stationary point need not be a minimum. The columns of A are linearly independent (rank A = M, so N ≥ M).If dropped, A^TWA is singular and non-invertible: the fit is degenerate (e.g. two identical basis functions), the solution is not unique, and a pseudo-inverse or regularisation is required to pick one. W is known independently of the data y.If the weights are estimated from the same residuals being fitted, the estimator becomes biased and the tidy covariance (A^TWA)⁻¹ understates the true variance.

Derivation

χ²(p) = (y − Ap)^T W (y − Ap)

Objective: the weighted sum of squared residuals. For Gaussian noise with covariance C and W = C⁻¹ this is −2 ln (likelihood) up to a constant, so minimising it maximises the likelihood. A

χ² = y^TWy − y^TWAp − p^TA^TWy + p^TA^TWAp

Multiply out the quadratic form, using (Ap)^T = p^TA^T. No approximation — pure distribution of the bilinear form. A

y^TWAp = p^TA^TWy ⇒ χ² = y^TWy − 2p^TA^TWy + p^TA^TWAp

Each middle term is a 1×1 scalar, so it equals its own transpose: (y^TWAp)^T = p^TA^TW^Ty, and W^T = W by symmetry. The two cross-terms are therefore identical and combine. B

∂(a^Tp)/∂p = a, ∂(p^TMp)/∂p = (M + M^T)p

Standard vector-calculus gradients. Here a = A^TWy and M = A^TWA. Because W is symmetric, M^T = A^TW^TA = M, so (M + M^T) = 2M. B

∂χ²/∂p = −2A^TWy + 2A^TWAp

Apply the two gradient rules from step 4 term by term; the constant y^TWy vanishes. A

∂χ²/∂p = 0 at p = p̂

Stationarity: at a minimum the gradient vanishes. Since χ² is a quadratic with Hessian 2A^TWA, which is positive-definite when A has full column rank and W > 0, the stationary point is the unique global minimum. B

−2A^TWy + 2A^TWAp̂ = 0 ⇒ (A^TWA)p̂ = A^TWy

Divide by 2 and move the data term to the right-hand side. These are the normal equations. A

p̂ = (A^TWA)⁻¹ A^TWy

Left-multiply by the inverse of A^TWA, which exists precisely under the full-rank assumption. A

Result

(A^TWA) p̂ = A^TWy ⇒ p̂ = (A^TWA)⁻¹A^TWy, Cov(p̂) = (A^TWA)⁻¹ (when W = C⁻¹)

Reading. The best-fit parameters are a linear map applied to the data. The matrix A^TWA (the "curvature" or Fisher information matrix) measures how sharply χ² rises away from the minimum; its inverse is the parameter covariance, so directions the data constrain poorly show up as large eigenvalues of (A^TWA)⁻¹. The right-hand side A^TWy is the weighted projection of the data onto each basis function.

Units check. Let observation i have units [y] and column j of A have units [y]/[p_j], and W = C⁻¹ have units [y]⁻². Then a diagonal entry (A^TWA)_jj has units ([y]/[p_j])²·[y]⁻² = [p_j]⁻², so (A^TWA)⁻¹_jj carries [p_j]² — correctly a variance in the parameter. And A^TWy has units ([y]/[p_j])·[y]⁻²·[y] = [p_j]⁻¹, so p̂_j comes out in [p_j]. Consistent.

Limiting cases

Uniform weights, W = I: reduces to ordinary least squares, A^TAp̂ = A^Ty.
Diagonal W = diag(1/σ_i²): independent measurements; each equation is weighted by inverse variance.
Single column of ones, A = (1,…,1)^T: the estimate collapses to the weighted mean p̂ = Σ_i w_iy_i / Σ_i w_i.
Square, invertible A (N = M): the fit interpolates exactly, p̂ = A⁻¹y, and the weights drop out.
One dominant weight w_k → ∞: the fit is forced through point k, an exact linear constraint on p.

Breaks when

Rank-deficient design (collinearity). If two columns of A are equal or nearly proportional, A^TWA is singular or ill-conditioned; the inverse blows up and the fitted parameters swing wildly on tiny data changes. Forming the normal equations squares the condition number, so even a moderately collinear problem loses precision — QR or SVD on A directly is safer.
Non-linear parameter dependence. If the model is y = f(p) with f non-linear (decay rates, frequencies, peak positions), there is no design matrix and χ² is not quadratic; the single-shot normal equations do not apply and one must linearise and iterate.
Outliers / non-Gaussian noise. The L₂ (squared) penalty assumes Gaussian errors; a single gross outlier can dominate A^TWy and drag the whole fit. Heavy-tailed data need robust (e.g. Huber, L₁) losses, which do not reduce to linear normal equations.
Wrong or misspecified weights. If W ≠ C⁻¹ (e.g. correlated errors treated as independent), p̂ stays unbiased but is no longer minimum-variance, and the reported covariance (A^TWA)⁻¹ is simply wrong.

Failure modes

Forgetting W is symmetric and keeping the two cross-terms separate, producing a spurious factor or an asymmetric gradient in step 3–4.
Solving by explicit inverse p̂ = (A^TWA)⁻¹A^TWy numerically instead of Cholesky-solving the linear system — wastes precision and can fail on ill-conditioned problems.
Putting the intercept in the wrong place: omitting the column of ones in A, which silently forces the fit through the origin.
Using w_i = 1/σ_i instead of 1/σ_i² — the weight is the inverse variance, not inverse standard deviation.
Reading (A^TWA)⁻¹ as the covariance when W ≠ C⁻¹; in general the covariance is the sandwich (A^TWA)⁻¹(A^TWCWA)(A^TWA)⁻¹.
Transposing the design matrix the wrong way so that A is M×N — then A^TWA has the wrong shape and the code errors or, worse, runs on a transposed problem.

Discussion

Geometrically the normal equations say the residual r̂ = y − Ap̂ is orthogonal to every column of A in the metric defined by W: A^TWr̂ = 0. The fit is the W-orthogonal projection of the data vector onto the column space of A. Best-fit does not mean the residual is zero; it means the residual has no component left inside the model subspace, so no adjustment of p could reduce it further. This projection picture is the reason least squares is so robust and universal — it is linear algebra, not statistics, once W is fixed.

The matrix A^TWA is the Fisher information for a Gaussian linear model. Its inverse saturates the Cramér–Rao bound, so with W = C⁻¹ the least-squares estimator is the best possible unbiased estimator — no cleverer method can do better with the same data. This connects three threads that look unrelated at first: minimising a sum of squares (numerical), maximising a Gaussian likelihood (statistical), and projecting orthogonally in a weighted inner product (geometric). They are one theorem seen from three sides.

The choice W = C⁻¹ is singled out by the Gauss–Markov theorem: among all linear unbiased estimators p̂ = Ly, the one with weight C⁻¹ has the smallest covariance in the positive-definite (Loewner) ordering, independently of whether the noise is actually Gaussian — Gaussianity is only needed to make it the maximum-likelihood and the minimum-variance estimator among all estimators. Any other symmetric positive-definite W still yields an unbiased p̂, but its covariance (A^TWA)⁻¹(A^TWCWA)(A^TWA)⁻¹ is larger, degrading gracefully rather than becoming wrong.

Common misconceptions. "Least squares assumes Gaussian noise" — not to compute p̂; the estimator is defined by the geometry and is unbiased for any zero-mean noise. Gaussianity is what makes it maximum-likelihood and lets you read off exact confidence regions. Also: a large final χ² does not mean the algebra failed; it means the model or the assumed error bars are wrong, and reduced χ² = χ²/(N−M) ≈ 1 is the target.

Worked examples

Example 1 — Unweighted straight-line fit. Fit y = c + mx to the four points (0, 1), (1, 3), (2, 2), (3, 5) with equal errors, so W = I.

A = [[1,0],[1,1],[1,2],[1,3]], p = (c, m)^T, y = (1,3,2,5)^T

First column of ones for the intercept, second column the x-values. A

A^TA = [[N, Σx], [Σx, Σx²]] = [[4, 6], [6, 14]]

Σx = 0+1+2+3 = 6, Σx² = 0+1+4+9 = 14, N = 4. A

A^Ty = [Σy, Σxy]^T = [11, 22]^T

Σy = 1+3+2+5 = 11, Σxy = 0+3+4+15 = 22. A

det = 4·14 − 6·6 = 20, (A^TA)⁻¹ = (1/20)[[14, −6], [−6, 4]]

2×2 inverse: swap diagonal, negate off-diagonal, divide by determinant. A

c = (14·11 − 6·22)/20 = 22/20 = 1.10, m = (−6·11 + 4·22)/20 = 22/20 = 1.10

Apply p̂ = (A^TA)⁻¹A^Ty. A

y = 1.10 + 1.10 x

Reading. Intercept and slope both 1.10 (dimensionless here). Residuals are (−0.1, 0.8, −1.3, 0.6), summing to 0 as the ones-column orthogonality guarantees.

Example 2 — Weighted mean of three measurements. Combine y = 10.2 ± 0.5, 9.8 ± 0.3, 10.5 ± 0.4 (independent). Here A = (1,1,1)^T and W = diag(1/σ_i²).

w_i = 1/σ_i²: w₁ = 1/0.25 = 4.00, w₂ = 1/0.09 = 11.11, w₃ = 1/0.16 = 6.25

Weights are inverse variances, not inverse standard deviations. A

A^TWA = Σw_i = 4.00 + 11.11 + 6.25 = 21.36

With A a single ones-column, A^TWA is the scalar sum of weights. A

A^TWy = Σw_iy_i = 4.00·10.2 + 11.11·9.8 + 6.25·10.5 = 215.3

Weighted projection of the data onto the constant model. A

μ̂ = A^TWy / A^TWA = 215.3 / 21.36 = 10.08

Normal equation in one dimension. A

σ_μ = (A^TWA)^−1/2 = 1/√21.36 = 0.216

Covariance is (A^TWA)⁻¹ = 1/21.36 = 0.0468, and σ is its square root. B

μ̂ = 10.08 ± 0.22

Reading. The combined estimate is pulled toward the most precise measurement (9.8 ± 0.3, weight 11.1) and its uncertainty 0.22 is smaller than any single error bar — the payoff of optimal weighting.

Problems

(A) Solve a given system. With A^TWA = [[10, 4], [4, 6]] and A^TWy = [20, 12]^T, find p̂.
Solution
det = 10·6 − 4·4 = 60 − 16 = 44. Inverse = (1/44)[[6, −4], [−4, 10]]. p̂₁ = (6·20 − 4·12)/44 = (120 − 48)/44 = 72/44 = 1.636. p̂₂ = (−4·20 + 10·12)/44 = (−80 + 120)/44 = 40/44 = 0.909. So p̂ = (1.64, 0.91)^T.
(B) Reduce to the weighted mean. Show directly from (A^TWA)p̂ = A^TWy that a one-parameter constant model gives p̂ = Σw_iy_i/Σw_i, and give its variance.
Solution
Take A = (1,…,1)^T (N×1) and W = diag(w_i). Then A^TWA = Σ_i w_i (scalar) and A^TWy = Σ_i w_iy_i. The normal equation (Σw_i)p̂ = Σw_iy_i gives p̂ = Σw_iy_i/Σw_i. Variance = (A^TWA)⁻¹ = 1/Σw_i = 1/Σ(1/σ_i²).
(B) Straight-line fit with data. Fit y = c + mx (equal weights) to (1, 2.1), (2, 3.9), (3, 6.2), (4, 7.8). Report m and c.
Solution
N = 4, Σx = 10, Σx² = 1+4+9+16 = 30, Σy = 2.1+3.9+6.2+7.8 = 20.0, Σxy = 2.1+7.8+18.6+31.2 = 59.7. A^TA = [[4,10],[10,30]], det = 120 − 100 = 20. m = (NΣxy − ΣxΣy)/det = (4·59.7 − 10·20.0)/20 = (238.8 − 200)/20 = 38.8/20 = 1.94. c = (Σx²Σy − ΣxΣxy)/det = (30·20.0 − 10·59.7)/20 = (600 − 597)/20 = 3/20 = 0.15. So y = 0.15 + 1.94x.
(C) Covariance identity. Assuming p̂ = (A^TWA)⁻¹A^TWy with y = Ap_true + n, Cov(n) = C, prove that when W = C⁻¹ the covariance of p̂ equals (A^TC⁻¹A)⁻¹.
Solution
Write L = (A^TWA)⁻¹A^TW. Then p̂ = Ly = LAp_true + Ln. Since LA = (A^TWA)⁻¹(A^TWA) = I, the estimator is unbiased: p̂ = p_true + Ln. Cov(p̂) = LCL^T = (A^TWA)⁻¹A^TW C WA(A^TWA)⁻¹. Set W = C⁻¹: the middle WCW = C⁻¹CC⁻¹ = C⁻¹ = W, so the expression becomes (A^TWA)⁻¹(A^TWA)(A^TWA)⁻¹ = (A^TWA)⁻¹ = (A^TC⁻¹A)⁻¹. ∎
(C) When is the solution unique? State and justify the condition on A and W for A^TWA to be invertible, and give a concrete rank-deficient example.
Solution
A^TWA is invertible iff it is positive-definite, i.e. iff v^TA^TWAv > 0 for all v ≠ 0. Writing u = Av, this is u^TWu > 0, which holds when W is positive-definite and u = Av ≠ 0 for every v ≠ 0 — that is, A must have full column rank M (linearly independent columns, requiring N ≥ M). Example: fitting y = p₁ + p₂ + p₃x gives columns (1,…), (1,…), (x,…); the first two are identical, so v = (1, −1, 0)^T satisfies Av = 0, A^TWA is singular, and p₁, p₂ are individually undetermined (only their sum is fixed).