Least-Squares Problems
Many practical problems lead to systems of equations \(A\mathbf{x} = \mathbf{b}\) that have no exact solution, typically because
the system is overdetermined (more equations than unknowns) and the data are noisy. The least-squares approach finds a vector
\(\hat{\mathbf{x}}\) that makes \(A\mathbf{x}\) as close as possible to \(\mathbf{b}\), providing the best approximate solution
in the sense of minimizing the squared error.
Definition: Least-Squares Solution
If \(A \in \mathbb{R}^{m \times n}\) and \(\mathbf{b} \in \mathbb{R}^m\), a least-squares solution of
\(A\mathbf{x} = \mathbf{b}\) is an \(\hat{\mathbf{x}} \in \mathbb{R}^n\) such that
\[
\| \mathbf{b} - A\hat{\mathbf{x}} \| \leq \| \mathbf{b} - A\mathbf{x} \| \quad \forall \mathbf{x} \in \mathbb{R}^n.
\]
The norm \(\| \mathbf{b} - A\hat{\mathbf{x}} \|\) is called the least-squares error of the approximation.
Derivation:
For \(A \neq 0\) the subspace \(\operatorname{Col} A \subseteq \mathbb{R}^m\) is nonzero, so
Basis Extension supplies it with a basis
and the Gram-Schmidt algorithm turns that
basis into an orthogonal one, which is what the two theorems quoted here require. By definition, \(\hat{\mathbf{x}}\) is a
least-squares solution if and only if \(A\hat{\mathbf{x}}\) is a closest point in
\(\operatorname{Col} A\) to \(\mathbf{b}\)
(minimizing \(\|\mathbf{b} - A\mathbf{x}\|\) over \(\mathbf{x} \in \mathbb{R}^n\) is equivalent to minimizing over
\(A\mathbf{x} \in \operatorname{Col} A\)). By the
Best Approximation Theorem, the unique
closest point is the orthogonal projection:
\[
A\hat{\mathbf{x}} = \operatorname{proj}_{\operatorname{Col} A} \mathbf{b} =: \hat{\mathbf{b}}.
\]
By the
Orthogonal Decomposition Theorem,
\(\mathbf{b} - \hat{\mathbf{b}} = \mathbf{b} - A\hat{\mathbf{x}}\) is orthogonal to \(\operatorname{Col} A\). The converse
holds as well, because that decomposition is unique: any \(\mathbf{x}\) for which \(\mathbf{b} - A\mathbf{x}\) is orthogonal
to \(\operatorname{Col} A\) satisfies \(A\mathbf{x} = \hat{\mathbf{b}}\). Let \(\mathbf{a}_j\) denote the \(j\)-th column of
\(A\). Orthogonality to every column means
\[
\mathbf{a}_j \cdot (\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{a}_j^\top(\mathbf{b} - A\hat{\mathbf{x}}) = 0 \quad (j = 1, \ldots, n).
\]
Stacking these \(n\) scalar equations into a single vector equation gives
\[
A^\top(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0},
\]
which rearranges to the normal equations
\[
A^\top A \hat{\mathbf{x}} = A^\top \mathbf{b}.
\]
When \(A^\top A\) is invertible (which holds precisely when the columns of \(A\) are
linearly independent, as the remark below shows),
we can solve explicitly:
\[
\begin{align*}
\hat{\mathbf{x}} &= (A^\top A)^{-1}A^\top\mathbf{b}, \\\\
\hat{\mathbf{b}} &= A\hat{\mathbf{x}} = A(A^\top A)^{-1}A^\top\mathbf{b}.
\end{align*}
\]
Theorem: Normal Equations
The set of least-squares solutions of \(A\mathbf{x} = \mathbf{b}\) consists of all vectors \(\hat{\mathbf{x}}\) that satisfy
the normal equations:
\[
A^\top A\mathbf{x} = A^\top\mathbf{b}.
\]
If the columns of \(A\) are linearly independent, then \(A^\top A\) is invertible and the unique least-squares solution is:
\[
\hat{\mathbf{x}} = (A^\top A)^{-1}A^\top\mathbf{b}.
\]
Remark: \(A^\top A\) is invertible if and only if the columns of \(A\) are linearly independent.
(\(\Leftarrow\)) Suppose the columns of \(A\) are linearly independent, and \(A^\top A\mathbf{x} = \mathbf{0}\). Multiplying
on the left by \(\mathbf{x}^\top\) gives
\[
\mathbf{x}^\top A^\top A \mathbf{x} = (A\mathbf{x})^\top (A\mathbf{x}) = \|A\mathbf{x}\|^2 = 0,
\]
so \(A\mathbf{x} = \mathbf{0}\). By
linear independence of the columns, the only
solution of \(A\mathbf{x} = \mathbf{0}\) is \(\mathbf{x} = \mathbf{0}\). Hence the
null space \(\operatorname{Nul}(A^\top A) = \{\mathbf{0}\}\), so
\(A^\top A\) is invertible by the
Invertible Matrix Theorem.
(\(\Rightarrow\)) Conversely, suppose \(A^\top A\) is invertible. If \(A\mathbf{x} = \mathbf{0}\), then
\(A^\top A\mathbf{x} = A^\top \mathbf{0} = \mathbf{0}\), which forces \(\mathbf{x} = \mathbf{0}\) by invertibility of
\(A^\top A\). Hence \(\operatorname{Nul} A = \{\mathbf{0}\}\), and the columns of \(A\) are therefore linearly independent.
QR Factorization Method:
In practice, forming \(A^\top A\) amplifies numerical errors (see the Computational Insight below). The
QR factorization gives a more stable
route. Assuming the columns of \(A\) are linearly independent, write \(A = QR\) where \(Q\) has orthonormal columns
(\(Q^\top Q = I\)) and \(R\) is upper triangular and invertible with positive diagonal.
Substituting into the normal equations \(A^\top A\hat{\mathbf{x}} = A^\top\mathbf{b}\) and simplifying,
\[
\begin{align*}
(QR)^\top(QR)\hat{\mathbf{x}} &= (QR)^\top\mathbf{b}, \\\\
R^\top \underbrace{Q^\top Q}_{=\,I} R \hat{\mathbf{x}} &= R^\top Q^\top \mathbf{b}, \\\\
R^\top R \hat{\mathbf{x}} &= R^\top Q^\top \mathbf{b}.
\end{align*}
\]
Since \(R\) is invertible,
so is \(R^\top\). Multiplying both
sides by \((R^\top)^{-1}\) yields
\[
R \hat{\mathbf{x}} = Q^\top \mathbf{b},
\]
and hence
\[
\hat{\mathbf{x}} = R^{-1} Q^\top \mathbf{b}.
\]
Note:
In practice, one solves \(R\hat{\mathbf{x}} = Q^\top\mathbf{b}\) by back-substitution rather than computing \(R^{-1}\) explicitly,
which is much faster and more numerically stable.
Computational Insight: Why QR Avoids Squaring the Condition Number
The QR factorization method avoids computing \(A^\top A\), which squares the condition number
of \(A\). This makes QR-based least-squares significantly more numerically stable, especially
for ill-conditioned problems. The back-substitution step is also computationally efficient
since \(R\) is upper triangular.
Moore-Penrose Pseudo-inverse (\(A^\dagger\))
While the normal equations \(A^\top A\hat{\mathbf{x}} = A^\top\mathbf{b}\) solve least-squares problems in principle, they have
serious limitations, one theoretical and one practical. The Moore-Penrose pseudo-inverse provides a unified framework that treats
full-rank, singular, and underdetermined systems alike and always singles out a unique solution. Combined with the singular value
decomposition (defined on a later page), with negligible singular values truncated, the same framework also copes
with ill-conditioned systems.
- Theoretical Failure:
If the columns of \(A\) are linearly dependent (for example, in an underdetermined system where \(n \gt m\), or due to
collinear features), \(A^\top A\) is singular and \((A^\top A)^{-1}\) does not exist.
- Practical Failure:
Even if \(A^\top A\) is invertible, it can be ill-conditioned. As defined in our treatment of the
singular value decomposition, the condition number
\(\kappa(A)\) measures a matrix's numerical sensitivity. The problem is that \(\kappa(A^\top A) = \kappa(A)^2\). This squaring of
the condition number is catastrophic. A moderately ill-conditioned problem (for example, \(\kappa(A) = 10^5\)) becomes a severely
ill-conditioned one (\(\kappa(A^\top A) = 10^{10}\)). Solving for \(\hat{\mathbf{x}}\) through \(A^\top A\) can therefore lose
about twice as many digits as the conditioning of \(A\) alone warrants.
Definition: Moore-Penrose Pseudo-inverse
For any matrix \(A \in \mathbb{R}^{m \times n}\), the Moore-Penrose pseudo-inverse
\(A^\dagger \in \mathbb{R}^{n \times m}\) is the unique matrix satisfying:
- \(A A^\dagger A = A\) (reconstruction property)
- \(A^\dagger A A^\dagger = A^\dagger\) (weak inverse property)
- \((A A^\dagger)^\top = A A^\dagger\) (symmetry of \(AA^\dagger\))
- \((A^\dagger A)^\top = A^\dagger A\) (symmetry of \(A^\dagger A\))
Exactly one matrix satisfies all four conditions for a given \(A\), and we take that existence and uniqueness as given. The power
of \(A^\dagger\) comes from what it does to an arbitrary right-hand side \(\mathbf{b}\).
Theorem: Minimum-Norm Least-Squares Solution
For any \(A \in \mathbb{R}^{m \times n}\) and \(\mathbf{b} \in \mathbb{R}^m\), the vector
\(\hat{\mathbf{x}} = A^\dagger \mathbf{b}\) is the unique minimum-norm least-squares solution
to \(A\mathbf{x} = \mathbf{b}\), meaning:
- Least-Squares Property:
\(\hat{\mathbf{x}}\) minimizes \(\|A\mathbf{x} - \mathbf{b}\|^2\).
- Minimum-Norm Property:
Among all vectors that minimize \(\|A\mathbf{x} - \mathbf{b}\|^2\), \(\hat{\mathbf{x}}\) has the smallest norm
\(\|\hat{\mathbf{x}}\|\).
Both properties follow from the four defining conditions. Conditions 1, 2 and 3 make \(AA^\dagger\) onto \(\operatorname{Col} A\),
idempotent, and symmetric, so it is the orthogonal projection onto \(\operatorname{Col} A\); hence
\(A\hat{\mathbf{x}} = AA^\dagger\mathbf{b}\) is the closest point of \(\operatorname{Col} A\) to \(\mathbf{b}\) and
\(\hat{\mathbf{x}}\) minimizes \(\|A\mathbf{x} - \mathbf{b}\|^2\). Every other minimizer has the form
\(\hat{\mathbf{x}} + \mathbf{z}\) with \(\mathbf{z} \in \operatorname{Nul} A\). Condition 2 gives
\(\hat{\mathbf{x}} = A^\dagger A\hat{\mathbf{x}}\) and condition 4 makes \(A^\dagger A\) symmetric, so
\(\mathbf{z}^\top\hat{\mathbf{x}} = (A^\dagger A\mathbf{z})^\top\hat{\mathbf{x}} = 0\), and therefore
\(\|\hat{\mathbf{x}} + \mathbf{z}\|^2 = \|\hat{\mathbf{x}}\|^2 + \|\mathbf{z}\|^2 \geq \|\hat{\mathbf{x}}\|^2\).
This property is particularly valuable in machine learning. For an underdetermined system (for example, more features than data
points), the minimum-norm solution is the "simplest" one, which is a form of implicit regularization.
The most stable way to compute \(A^\dagger\) is via the Singular Value Decomposition (SVD), which we define in
our treatment of symmetric matrices and the SVD.
Linear Regression
Least-squares theory has its most widespread application in statistical modeling through linear regression. Given data points, we
seek the best linear relationship between predictor variables and an observed response. This is a direct application of the
least-squares framework we have developed.
Definition: Linear Regression Model
Given observed data points \(\{(\mathbf{x}_i, y_i)\}_{i=1}^n\) where \(\mathbf{x}_i \in \mathbb{R}^d\) and
\(y_i \in \mathbb{R}\), the linear regression model is:
\[
\mathbf{Y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon}
\]
where:
- \(X \in \mathbb{R}^{n \times d}\) is the design matrix
- \(\boldsymbol{\beta} \in \mathbb{R}^d\) is the parameter (that is, weight) vector
- \(\mathbf{Y} \in \mathbb{R}^n\) is the observation vector
- \(\boldsymbol{\epsilon} = \mathbf{Y} - X\boldsymbol{\beta}\) is the error (noise) vector; for an estimate
\(\hat{\boldsymbol{\beta}}\), \(\hat{\boldsymbol{\epsilon}} = \mathbf{Y} - X\hat{\boldsymbol{\beta}}\) is the
residual vector
The dimension \(d\) is the number of features, and \(n\) is the number of data points. The residual
\(\hat{\boldsymbol{\epsilon}}\) collects the "differences" between each observed value \(y_i\) and its corresponding predicted
\(y\) value. The least-squares hyperplane represents the set of predicted \(y\) values based on estimated
parameters \(\hat{\boldsymbol{\beta}}\), and \(\hat{\boldsymbol{\beta}}\) must satisfy the normal equations:
\[
X^\top X\hat{\boldsymbol{\beta}} = X^\top\mathbf{Y}.
\]
Note:
The linear model is linear in terms of parameters \(\boldsymbol{\beta}\), not \(X\). We can choose any non-linear
transformation for each \(\mathbf{x}_i\). For example, \(y = \beta_0 + \beta_1 x^2 + \beta_2 \sin(2\pi x)\) is still a linear
model. Thus, \(\mathbf{Y}\) is modeled as a
linear combination of features (that
is, predictors) \(X\) with respect to the coefficients \(\boldsymbol{\beta}\).
Many details remain in this topic. We take them up in
the probability and statistics page on linear regression. Finally,
\(X^\top X\) is a symmetric matrix. We will learn what makes such matrices special.
Interactive Linear Regression Demo
The demo below fits polynomials (and, in 3D mode, surfaces) to data by least squares. It also practices what this page preaches.
The solver runs through the QR factorization (Gram-Schmidt on the columns of the design matrix, then
back-substitution on \(R\hat{\boldsymbol{\beta}} = Q^\top \mathbf{Y}\)) rather than forming \(X^\top X\). On ill-conditioned
problems, however, Gram-Schmidt as implemented here can itself lose accuracy, so avoiding the condition-number squaring of the
Computational Insight above in floating point needs a more careful QR-based solve, for example running the same Gram-Schmidt with
\(\mathbf{Y}\) appended as an extra column of \(X\).
Two readouts turn the theory into checkable numbers. First, \(\|X^\top \hat{\boldsymbol{\epsilon}}\|\) is displayed after every fit. Its
vanishing says that the residual is orthogonal to every column of \(X\), which is exactly the statement of the normal equations,
that is, \(\hat{\mathbf{Y}} = \operatorname{proj}_{\operatorname{Col} X} \mathbf{Y}\). Second, toggle "Show squared residuals" to see
the geometry behind the name. Each residual becomes the side of a square, and the fit is the choice of coefficients minimizing the
total area of those squares.
The demo is also honest about when the theory says no. With fewer points than parameters (\(n \lt p\)), or with fewer
than \(p\) distinct \(x\)-values, the columns of \(X\) are linearly dependent, \(X^\top X\) is singular, and there is no unique
least-squares solution. The demo reports this (pointing back to the pseudo-inverse section) instead of plotting garbage. At
exactly \(n = p\) it flags the interpolation regime, where the residual is zero and "fitting" degenerates into passing through
every point. Watch a quartic do this with five points.
In 3D mode, compare the plane and quadratic surface models on the saddle dataset. The plane underfits (large residuals in a
telltale pattern), yet \(\|X^\top \hat{\boldsymbol{\epsilon}}\|\) stays at zero up to rounding. It is still the best plane,
just not a good model.