Least Squares

Least squares is the method of choosing model parameters so that the sum of squared differences between the model’s predictions and the measurements is as small as possible. It is the standard way to solve an overdetermined system, one with more equations than unknowns, in which noise makes the equations mutually inconsistent. In computer vision it fits lines and curves to points, estimates motion in the Lucas–Kanade method, initializes homographies and camera matrices, and, in its nonlinear form, refines entire 3D reconstructions.

Adrien-Marie Legendre published the method in 1805, in an appendix to a book on computing the orbits of comets. Carl Friedrich Gauss published his own treatment in 1809, stating that he had used the principle since 1795 and deriving it as the most probable estimate under what is now called the Gaussian distribution of errors. Their priority dispute is the most famous in the history of statistics. Stigler (1981) weighed the documentary and numerical evidence and concluded, while stressing that it is not decisive, that Gauss probably did use the method before Legendre but failed to make it known to his contemporaries.

Definition

Given $m$ measurements and a model with $n < m$ unknown parameters that predicts each measurement, the least-squares estimate is the parameter vector that minimizes the sum of the squared residuals, the differences between measured and predicted values. When the predictions are linear in the parameters, the problem is linear least squares and has a closed-form solution; otherwise it is nonlinear least squares and is solved iteratively.

Intuition

Fitting a line to twenty noisy points gives twenty equations $y_i = a x_i + c$ in two unknowns. No line passes through all the points, so the question is how to trade errors off against each other. Squaring makes every error count positively and produces a smooth, bowl-shaped objective whose minimum is found by setting the derivative to zero, which yields a linear system. The same squaring has a cost: a single point with a large error contributes so heavily to the sum that it can drag the whole fit toward itself.

Formal Definition

Write the linear model as $A\mathbf{x} \approx \mathbf{b}$, with design matrix $A \in \mathbb{R}^{m \times n}$, unknowns $\mathbf{x} \in \mathbb{R}^n$, measurements $\mathbf{b} \in \mathbb{R}^m$, and $m > n$. The least-squares solution minimizes the squared norm of the residual $\mathbf{r} = A\mathbf{x} - \mathbf{b}$:

$$\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \lVert A\mathbf{x} - \mathbf{b} \rVert^2 = \arg\min_{\mathbf{x}} \sum_{i=1}^{m} \big( \mathbf{a}_i^\top \mathbf{x} - b_i \big)^2,$$

where $\mathbf{a}_i^\top$ is the $i$-th row of $A$. Setting the gradient $2A^\top(A\mathbf{x} - \mathbf{b})$ to zero gives the normal equations

$$A^\top A \, \hat{\mathbf{x}} = A^\top \mathbf{b}.$$

If $A$ has full column rank, $A^\top A$ is invertible and the solution is unique, $\hat{\mathbf{x}} = (A^\top A)^{-1} A^\top \mathbf{b}$. If it does not, for example when the data cannot distinguish two parameters, infinitely many minimizers exist, and the pseudoinverse $A^{+}$ selects the one of smallest norm, $\hat{\mathbf{x}} = A^{+}\mathbf{b}$.

Geometric view

The products $A\mathbf{x}$ for all $\mathbf{x}$ form the column space of $A$, an $n$-dimensional subspace of $\mathbb{R}^m$. The measurement vector $\mathbf{b}$ generally lies outside it, and the closest point of the subspace to $\mathbf{b}$ is its orthogonal projection $A\hat{\mathbf{x}}$. The normal equations say exactly that the residual is orthogonal to every column of $A$, $A^\top(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}$. The projection onto the column space is $P = A(A^\top A)^{-1}A^\top$.

Properties

Solving it in practice

Forming $A^\top A$ explicitly squares the condition number, $\kappa(A^\top A) = \kappa(A)^2$, so on an ill-conditioned problem the normal equations can lose about twice as many significant digits as methods that work with $A$ directly, especially when the residual is small. Numerical libraries therefore factor $A$ itself:

• QR decomposition. Factor $A = QR$ with orthonormal columns in $Q$ and upper-triangular $R$; then $\hat{\mathbf{x}}$ solves the triangular system $R\hat{\mathbf{x}} = Q^\top \mathbf{b}$. This is the usual method for well-conditioned, full-rank problems.
• Singular value decomposition. Factor $A = U S V^\top$ with the singular values on the diagonal of $S$; then $\hat{\mathbf{x}} = V S^{+} U^\top \mathbf{b}$, where $S^{+}$ inverts the nonzero singular values. The SVD is slower but reveals rank deficiency, and small singular values can be truncated to stabilize nearly singular problems.

The normal equations remain acceptable when $A$ is well conditioned and $n$ is tiny, as in the $2 \times 2$ system of Lucas–Kanade, and are often used deliberately when $A^\top A$ is sparse and structured. Centering and scaling the data first, so that columns have similar magnitudes, is cheap and often improves conditioning decisively.

Weighted least squares

When measurements have different reliabilities, each squared residual is weighted:

$$\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} (A\mathbf{x} - \mathbf{b})^\top W (A\mathbf{x} - \mathbf{b}), \qquad A^\top W A \, \hat{\mathbf{x}} = A^\top W \mathbf{b}.$$

With a diagonal $W$, each equation is simply scaled by $\sqrt{w_i}$. The usual choice is the inverse of the noise covariance, $W = \Sigma^{-1}$, so that precise measurements count more. The Gaussian window $g$ in Lucas–Kanade is a weighting of this kind, by distance from the tracked point rather than by noise.

Regularized least squares

When $A$ is ill conditioned or the problem is underdetermined, adding a penalty on the solution makes it well posed:

$$\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \lVert A\mathbf{x} - \mathbf{b} \rVert^2 + \lambda \lVert \Gamma \mathbf{x} \rVert^2, \qquad (A^\top A + \lambda \Gamma^\top \Gamma)\, \hat{\mathbf{x}} = A^\top \mathbf{b}.$$

With $\Gamma = I$ this is ridge regression; with a general $\Gamma$, for example a derivative operator that penalizes rough solutions, it is Tikhonov regularization. The Horn–Schunck method has this structure: a quadratic data term from brightness constancy plus a quadratic smoothness term on the flow field.

Statistical view

Suppose the measurements are $\mathbf{b} = A\mathbf{x} + \boldsymbol{\epsilon}$ with independent Gaussian noise of equal variance $\sigma^2$. The likelihood of $\mathbf{x}$ is proportional to $\exp\!\big(-\lVert A\mathbf{x} - \mathbf{b} \rVert^2 / 2\sigma^2\big)$, so maximizing it is the same as minimizing the sum of squares: least squares is the maximum-likelihood estimate under Gaussian noise. With correlated Gaussian noise of covariance $\Sigma$, the maximum-likelihood estimate is weighted least squares with $W = \Sigma^{-1}$. The estimate’s covariance is $\sigma^2 (A^\top A)^{-1}$, which is how the aperture problem article obtains the elongated uncertainty of flow estimated near an edge.

The Gauss–Markov theorem gives a weaker but distribution-free guarantee: if the noise has zero mean and is uncorrelated with equal variance, ordinary least squares has the smallest variance among all unbiased estimators that are linear in $\mathbf{b}$. No Gaussian assumption is needed, but a nonlinear or biased estimator can still do better.

Nonlinear least squares

Most geometric problems in vision, such as projecting 3D points through a camera, are nonlinear in their parameters. Minimizing $\sum_i r_i(\mathbf{x})^2$ then proceeds iteratively: the Gauss–Newton method linearizes the residuals around the current estimate and solves a linear least-squares problem for the update, and the Levenberg–Marquardt method adds a damping term that blends Gauss–Newton with gradient descent when the linearization is poor. Each iteration is a linear least-squares solve, so everything above applies inside the loop. Iterated Lucas–Kanade is a Gauss–Newton method.

Sensitivity to outliers

Because the penalty grows quadratically, least squares is fragile when some measurements are simply wrong, such as mismatched features or pixels that violate brightness constancy. Robust alternatives replace the square with a penalty that grows more slowly, solved by iteratively reweighted least squares, or separate inliers from outliers by random sampling (RANSAC) and fit only the inliers.

Examples

Line fitting

The code fits a line to twenty noisy points, once with NumPy’s SVD-based lstsq and once through the normal equations, then corrupts a single point.

import numpy as np

rng = np.random.default_rng(0)

# 20 noisy samples of the line y = 2x + 1.
x = np.linspace(0, 10, 20)
y = 2 * x + 1 + rng.normal(0, 0.5, size=x.size)

A = np.column_stack([x, np.ones_like(x)])   # one row [x_i, 1] per point

theta_lstsq, *_ = np.linalg.lstsq(A, y, rcond=None)   # SVD-based solver
theta_normal = np.linalg.solve(A.T @ A, A.T @ y)       # normal equations
print("lstsq            slope, intercept:", np.round(theta_lstsq, 3))
print("normal equations slope, intercept:", np.round(theta_normal, 3))
print(f"cond(A) = {np.linalg.cond(A):.1f}, cond(A^T A) = {np.linalg.cond(A.T @ A):.1f}")

# Corrupt a single point with a gross error and refit.
y_bad = y.copy()
y_bad[-1] += 30
theta_bad, *_ = np.linalg.lstsq(A, y_bad, rcond=None)
print("with one outlier slope, intercept:", np.round(theta_bad, 3))

Output:

lstsq            slope, intercept: [1.974 1.04 ]
normal equations slope, intercept: [1.974 1.04 ]
cond(A) = 11.5, cond(A^T A) = 132.6
with one outlier slope, intercept: [ 2.788 -1.531]

Both solvers agree on this well-conditioned problem, and the printed condition numbers confirm $\kappa(A^\top A) = \kappa(A)^2$. One corrupted point out of twenty moves the slope by about 40% and the intercept by more than 2.5 units, which is why vision pipelines pair least squares with outlier rejection.

Motion of an image window

Lucas–Kanade collects one linearized brightness constancy equation $I_x u + I_y v + I_t = 0$ per pixel of a window, giving an overdetermined system in the two unknowns $(u, v)$. Its normal equations are a $2 \times 2$ system whose matrix is the structure tensor of the window; when that matrix is singular or nearly so, as on an edge or in a flat region, the flow is undetermined, which is the aperture problem.

Homographies and camera matrices

Each point correspondence gives two equations that are linear in the nine entries of a $3 \times 3$ homography. Stacking them yields a homogeneous system $A\mathbf{h} = \mathbf{0}$; the direct linear transform (DLT) takes the unit vector $\mathbf{h}$ that minimizes $\lVert A\mathbf{h} \rVert$, the right singular vector of $A$ with the smallest singular value. The same construction estimates a $3 \times 4$ camera projection matrix from 3D–2D correspondences. Hartley and Zisserman stress that normalizing the coordinates first is essential, and that the minimized quantity is an algebraic error without direct geometric meaning. The linear solution therefore usually initializes a nonlinear least-squares refinement of a geometric error. Camera calibration follows the same pattern: a closed-form estimate, then nonlinear refinement of the reprojection error.

Bundle adjustment

Bundle adjustment refines all camera parameters and 3D points of a reconstruction together by minimizing reprojection error, the distance between each observed image point and the projection of its 3D point. It is a large nonlinear least-squares problem, solved by Gauss–Newton or Levenberg–Marquardt iterations that exploit the sparsity of the normal equations, since each residual involves only one camera and one point. Triggs and colleagues’ survey covers these methods and the robust cost functions often used in place of the plain square.

Common Misconceptions

• “Least squares assumes Gaussian noise.” The method can be applied to any data; Gaussian noise is what makes it the maximum-likelihood estimate, and the Gauss–Markov theorem gives a guarantee without it. What it does assume implicitly is that large errors are rare, which outliers violate.
• “Solve it with the inverse of $A^\top A$.” Mathematically correct, but numerically the most fragile choice; prefer a QR- or SVD-based solver.
• “The residual to minimize is given by the data.” It is a modeling choice. Algebraic residuals of a DLT, vertical distances in line fitting, and perpendicular distances (total least squares) give different answers when the noise is not where the model assumes it is.
• “A small residual means a good estimate.” A nearly rank-deficient $A$, as on an edge in Lucas–Kanade, gives small residuals and poorly determined parameters; the covariance $\sigma^2 (A^\top A)^{-1}$ measures the latter.

Where It Is Used

• Motion estimation: Lucas–Kanade solves a small least-squares problem per window; variational methods such as Horn–Schunck minimize regularized quadratic energies over the whole image.
• State estimation: the Kalman filter update solves a weighted least-squares problem balancing prediction and measurement by their inverse covariances.
• Geometry: homography, fundamental-matrix, and camera-matrix estimation, calibration, triangulation, and bundle adjustment.
• Model fitting: lines, circles, conics, distortion models, and image alignment.