Least Squares | LA 101

The Overdetermined Problem

In science and engineering, measurements are abundant and unknowns are few. A GPS receiver takes dozens of satellite range readings to pin down three coordinates. A data scientist fits a line to hundreds of noisy data points. A communications engineer estimates a channel from more pilot symbols than channel taps. In all these cases, we have an overdetermined system Ax = b — more rows than columns — where b lies outside the column space of A and no exact solution exists.

Instead of demanding Ax = b exactly, we ask for the best approximation: find x̂ that minimizes the residual r = b − Ax̂. "Best" means minimizing the squared length of r — hence least squares.

Least Squares Problem

\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \|\mathbf{b} - A\mathbf{x}\|^2

We seek x̂ that minimizes the sum of squared residuals ‖b − Ax‖². A is m×n with m > n (overdetermined). The minimum is achieved when the residual r = b − Ax̂ is orthogonal to every column of A.

Geometric View: Projection onto the Column Space

Here is the key geometric insight. The column space C(A) is a subspace of ℝᵐ. The vector b may not lie in C(A). The closest point to b that does lie in C(A) is the orthogonal projection of b onto C(A), which we call p = Ax̂. The residual r = b − p is then perpendicular to C(A).

Perpendicularity to C(A) means r is orthogonal to every column of A, i.e., Aᵀr = 0. Substituting r = b − Ax̂:

Normal Equations

A^T A\,\hat{\mathbf{x}} = A^T \mathbf{b} \quad \Longrightarrow \quad \hat{\mathbf{x}} = (A^T A)^{-1} A^T \mathbf{b}

Aᵀ(b − Ax̂) = 0 expands to AᵀAx̂ = Aᵀb. These are the normal equations. When A has full column rank, AᵀA is invertible and the unique solution is x̂ = (AᵀA)⁻¹Aᵀb.

Why "Normal"?

The word normal here means perpendicular — the residual is normal (perpendicular) to the column space. It has nothing to do with Gaussian distributions, though the least squares solution coincides with the maximum likelihood estimate when errors are normal (Gaussian).

The Pseudo-Inverse and the Projection Matrix

The matrix A⁺ = (AᵀA)⁻¹Aᵀ is called the pseudo-inverse (or Moore-Penrose inverse) of A. It satisfies x̂ = A⁺b and generalizes the concept of matrix inversion to non-square matrices.

The vector p = Ax̂ = A(AᵀA)⁻¹Aᵀb is the projection of b onto C(A). The matrix that performs this projection is:

Projection Matrix

P = A(A^T A)^{-1} A^T, \qquad P^2 = P, \quad P^T = P

P = A(AᵀA)⁻¹Aᵀ is the orthogonal projector onto C(A). Key properties: P² = P (idempotent — projecting twice changes nothing), Pᵀ = P (symmetric), and (I − P) is the projector onto the orthogonal complement of C(A). ‖b‖² = ‖Pb‖² + ‖(I−P)b‖² (Pythagorean theorem).

Full-Rank Assumption

The formula x̂ = (AᵀA)⁻¹Aᵀb requires AᵀA to be invertible, which happens exactly when A has full column rank (rank n). If the columns of A are linearly dependent, the least squares solution is not unique — any x̂ in a particular affine subspace achieves the same minimum residual. In that case we typically seek the minimum-norm solution via the full pseudo-inverse or regularization.

Linear Regression as Least Squares

Linear regression — fitting a line, plane, or polynomial to data — is exactly least squares. Given m data points (tᵢ, bᵢ), fitting the model b = x₁ + x₂t leads to the system:

Regression System

A = \begin{pmatrix} 1 & t_1 \\ 1 & t_2 \\ \vdots & \vdots \\ 1 & t_m \end{pmatrix}, \quad \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}, \quad \mathbf{b} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{pmatrix}

A is the Vandermonde-style design matrix with a column of ones (intercept) and a column of t values (slope). The normal equations AᵀAx̂ = Aᵀb solve for the intercept x̂₁ and slope x̂₂ that minimize the total squared vertical deviation from the line.

This extends naturally to polynomial regression (columns 1, t, t², …, tᵏ), multiple regression (multiple predictors), and weighted regression (minimize ‖W(b − Ax)‖² for a diagonal weight matrix W).

The Residual Sum of Squares

After finding x̂, the minimum value of the objective is ‖b − Ax̂‖² = ‖(I − P)b‖², the squared distance from b to its projection. Dividing by the degrees of freedom m − n gives the mean squared error — the standard estimate of noise variance in linear regression.

QR Factorization: The Numerically Stable Approach

Forming AᵀA and solving the normal equations directly is conceptually clean but numerically problematic. The condition number of AᵀA is the square of the condition number of A — so if A is even mildly ill-conditioned, forming AᵀA doubles the loss of significant digits.

The QR factorization of A avoids this. Write A = QR where Q is m×n with orthonormal columns (Q in the thin/reduced sense) and R is n×n upper triangular. Then:

QR Solution to Least Squares

A = QR \quad \Rightarrow \quad R\hat{\mathbf{x}} = Q^T \mathbf{b}

A = QR implies AᵀA = RᵀQᵀQR = RᵀR and Aᵀb = RᵀQᵀb. The normal equations simplify to Rx̂ = Qᵀb — a triangular system solved by back-substitution. This avoids forming AᵀA entirely and has condition number κ(A) instead of κ(A)².

This is why all serious numerical software (LAPACK, MATLAB's backslash, NumPy's lstsq) solves overdetermined systems via QR or SVD — never by directly forming the normal equations.

SVD: The Ultimate Least Squares Tool

The Singular Value Decomposition A = UΣVᵀ gives the minimum-norm least squares solution x̂ = VΣ⁺Uᵀb for any A — even rank-deficient ones. It is the most general and robust approach, used when A may not have full column rank or when numerical conditioning is paramount.

Applications in Signal Processing and Communications

Channel Estimation

A wireless channel is estimated by sending a known pilot sequence. If the channel impulse response has n taps and we send m > n pilots, the received signal satisfies y ≈ Φh where Φ is a convolution matrix and h is the unknown channel vector. The least squares estimate ĥ = (ΦᵀΦ)⁻¹Φᵀy is the standard pilot-based channel estimator.

Beamforming Weight Design

In antenna array processing, we want beamforming weights w that approximate a desired spatial response d(θ) at sampled directions. The steering vector matrix A assembles the responses; the least squares weights minimize ‖d − Aw‖².

System Identification

Given input-output measurements from an unknown linear system, fitting an AR or FIR model is a least squares problem. The design matrix is built from delayed input samples; the coefficient vector is estimated via normal equations or QR.

Key Takeaways

An overdetermined system Ax = b (m > n) has no exact solution when b ∉ C(A). The least squares solution x̂ minimizes ‖b − Ax‖² and satisfies the normal equations AᵀAx̂ = Aᵀb — geometrically, it projects b onto C(A). The projection matrix P = A(AᵀA)⁻¹Aᵀ satisfies P² = P and Pᵀ = P. Linear regression is least squares with a design matrix. For numerical stability, solve via QR (Rx̂ = Qᵀb) rather than forming AᵀA. The SVD handles rank-deficient cases and gives the minimum-norm solution.