Least Squares

The best-fit solution to an overdetermined system, minimizing the sum of squared residuals.

The problem

When a system Ax = b has more equations than unknowns, it usually has no exact solution. The least-squares approach instead finds the x that makes Ax as close to b as possible, minimizing the squared length of the residual |Ax - b|^2. Geometrically, Ax is the orthogonal projection of b onto the column space of A.

The normal equations

Setting the gradient of the squared residual to zero gives the normal equations A^T A x = A^T b. When A has independent columns, A^T A is symmetric positive definite and the solution is unique. These equations express that the residual must be orthogonal to the column space, the defining property of a projection.

Solving it well

Forming and solving the normal equations is simple but squares the condition number of A, which can wreck accuracy. The numerically preferred route factors A = QR and solves R x = Q^T b, or uses the SVD for rank-deficient or ill-conditioned problems. Software such as lstsq chooses a stable method automatically.

Where it appears

• fitting a model curve to noisy measurements
• linear regression in statistics and machine learning
• reconstructing a signal from more observations than parameters
• calibrating simulation parameters against experimental data

import numpy as np
# Fit a line y = m x + c to noisy data
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([0.1, 0.9, 2.1, 2.9])
A = np.vstack([x, np.ones_like(x)]).T
(m, c), *_ = np.linalg.lstsq(A, y, rcond=None)
print(m, c)   # slope near 1, intercept near 0

Least squares is how measured diagnostic signals are reconciled with model predictions, yielding the parameter estimates and uncertainties that ground a simulation in observation.