Least Squares as Orthogonal Projection onto a Design Matrix

A design matrix $X$ with $n$ rows and $p$ columns defines a $p$-dimensional subspace of $\mathbb{R}^n$: the set of all vectors $X\beta$. Fitting a linear model by least squares means finding the point in that subspace closest to the observed response $y$. The closest point in a subspace is reached by dropping a perpendicular, so the fitted vector is the orthogonal projection of $y$ onto the column space of $X$. Every property that follows (the normal equations, the behaviour under collinearity, the sensitivity to scaling, the meaning of residual degrees of freedom) is a restatement of that one geometric fact.

The projection view also explains why the obvious algebraic route is the wrong one numerically. Writing down $X^\top X \beta = X^\top y$ and solving it is the shortest path from the optimisation problem to a formula, and it is also the path that discards roughly half the available precision. The two facts are connected: the normal equations are the projection condition expressed in a basis that is badly conditioned whenever the columns of $X$ are close to each other.

The residual is perpendicular to every column

Minimising $\lVert y - X\beta \rVert_2^2$ expands to $y^\top y - 2\beta^\top X^\top y + \beta^\top X^\top X \beta$. The quadratic is convex, so a stationary point is the global minimum, and differentiating with respect to $\beta$ gives the stationarity condition

$$X^\top X \beta = X^\top y.$$

Read as geometry rather than algebra, this says $X^\top (y - X\beta) = 0$: the residual $r = y - X\hat\beta$ has zero inner product with each column of $X$, hence with the entire column space. If some direction $Xv$ in the column space had a nonzero inner product with $r$, then moving the coefficient vector a small distance along $v$ would shorten the residual, contradicting optimality. For a quadratic the first-order condition is exact, so no second-order check is needed.

Two consequences follow immediately. First, the decomposition is Pythagorean: $\lVert y \rVert^2 = \lVert X\hat\beta \rVert^2 + \lVert r \rVert^2$, which is the algebraic identity behind the coefficient of determination. Second, the fitted values are a linear function of the response, $\hat y = Hy$ with

$$H = X(X^\top X)^{-1}X^\top .$$

$H$ is symmetric and idempotent (that is, $H^2 = H$), which is exactly the statement that projecting twice changes nothing. Its rank equals $\operatorname{rank}(X)$, and since the trace of an idempotent matrix equals its rank, $\operatorname{tr}(H) = p$ for a full-rank design. The diagonal entries $h_{ii}$, the leverages, therefore sum to $p$ and lie in $[0,1]$. An observation with $h_{ii} = 1$ is fitted exactly no matter what its response value is, which is the algebraic form of “this row is a singleton in the column space”. The residual degrees of freedom, $n - p$, is the dimension of the orthogonal complement that the residual lives in, not a bookkeeping convention.

Without an intercept column, the projection is onto a subspace that excludes the constant vector, so the residual is orthogonal to the columns but not to $\mathbf{1}$. The mean residual is then generally nonzero, and the usual centred $R^2$ is not comparable across models. This is the whole reason a column of ones is prepended: it makes the constant part of $y$ part of the projection rather than part of the residual.

Rank deficiency leaves the fit unique but the coefficients not

If the columns of $X$ are linearly dependent, $X^\top X$ is singular and the normal equations have infinitely many solutions. The fitted vector does not share that ambiguity. The column space is unchanged by adding a vector from the null space to $\beta$, so every minimiser produces the same $\hat y$ and the same residual sum of squares. Only the coefficients are unidentified, and only in the null-space directions.

The standard resolution is the Moore–Penrose pseudo-inverse, computed from the singular value decomposition $X = U\Sigma V^\top$. The solution $\hat\beta = V\Sigma^{+}U^\top y$ sets to zero every coefficient component along a direction whose singular value falls below a threshold, and solves the rest exactly. Among all minimisers it returns the one of smallest Euclidean norm. That is a choice, not a theorem about the data: the minimum-norm solution is the one that a shrinkage prior centred at zero would prefer, and a different prior would select a different member of the same affine set.

Any implementation therefore needs a rank decision, and the decision is a relative threshold on singular values. numpy.linalg.lstsq exposes it as rcond, interpreted relative to the largest singular value; a value that is too generous silently discards directions that carry real signal, and a value that is too strict amplifies directions that carry only rounding error. Neither failure announces itself. The first produces a model that underfits in a way the training loss barely registers, and the second produces coefficients whose magnitudes are unrelated to anything in the data.

import numpy as np

def fit_via_svd(X, y, rcond=1e-12):
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    keep = s > rcond * s[0]          # relative cutoff, not absolute
    inv_s = np.zeros_like(s)
    inv_s[keep] = 1.0 / s[keep]
    beta = Vt.T @ (inv_s * (U.T @ y))
    return beta, keep.sum()          # effective rank, useful as a diagnostic

The effective rank returned alongside the coefficients is worth logging. A design that was meant to have $p$ independent columns and reports a smaller number has a duplicated or constant feature somewhere, and the coefficient table will not reveal which one.

Why the Gram matrix halves the precision

The condition number of a matrix measures how much a relative perturbation in the input can be amplified in the solution. In the 2-norm, forming $X^\top X$ squares it: $\kappa_2(X^\top X) = \kappa_2(X)^2$. Since the columns of $X^\top X$ are inner products of the columns of $X$, this is the familiar statement that the cosine between two nearly parallel columns is a number close to one, and the difference from one is a small quantity computed by cancelling large ones.

Double precision carries a 53-bit significand, so the unit roundoff is about $1.1 \times 10^{-16}$. A design with $\kappa_2(X) = 10^4$, unremarkable for unscaled features measured in different units, has $\kappa_2(X^\top X) = 10^8$, and the forward error in $\hat\beta$ from the normal equations scales like $\kappa^2 u \approx 10^{-8}$: eight digits survive out of sixteen. At $\kappa_2(X) = 10^8$ the squared condition number exceeds $1/u$ and the computed Gram matrix can be numerically singular even though $X$ is not. The coefficients are then meaningless while the fitted values may still look reasonable, because the ill-conditioned directions are exactly the ones the projection suppresses.

Orthogonalisation avoids the squaring. Householder QR costs about $2np^2 - \tfrac{2}{3}p^3$ floating-point operations and produces $X = QR$ with $Q$ orthogonal, so $\kappa_2(R) = \kappa_2(X)$. Solving $R\beta = Q^\top y$ by back substitution has forward error of order $\kappa_2(X)u$ rather than its square. For $n \gg p$ the extra cost over the Gram approach, which needs about $np^2$ operations for the cross products alone, is a factor of roughly two in flops. It is paid once, in exchange for doubling the number of correct digits in the worst case.

import numpy as np

def fit_via_qr(X, y):
    Q, R = np.linalg.qr(X)           # reduced QR, Q is n x p
    beta = np.linalg.solve(R, Q.T @ y)   # R is upper triangular
    resid = y - X @ beta
    return beta, np.linalg.norm(resid)

The trade does not always go to QR. The normal equations can be accumulated in a single streaming pass over the data: adding one row updates $X^\top X$ with a rank-one change and $X^\top y$ with one inner product, and the design matrix never has to be held in memory. QR needs the whole matrix, or a carefully maintained orthogonal factor that costs more per row than a rank-one update. For problems where $n$ is large enough that the data cannot be revisited, and $p$ is small enough that the conditioning is comfortable, the Gram matrix remains the practical choice. The decision is between a factor-of-two flop count and the ability to process the data once, and it should be made by looking at the actual condition number rather than by habit.

Where the projection picture stops being useful

Three variants modify the geometry without changing its character. Weighted least squares minimises $(y - X\beta)^\top W (y - X\beta)$ for a positive diagonal $W$, which is an orthogonal projection in the inner product defined by $W$; equivalently it is ordinary least squares on the rows of $X$ and $y$ scaled by $\sqrt{w_i}$. Ridge regression adds $\lambda \lVert \beta \rVert^2$ and is no longer a projection onto the column space at all: the solution $(X^\top X + \lambda I)^{-1}X^\top y$ shrinks each singular direction by the factor $s_i^2/(s_i^2 + \lambda)$, and the estimator is no longer unbiased; the mechanics are in ridge shrinkage in ill-conditioned designs.

The remaining limitation is statistical rather than numerical. Least squares is the minimum-variance linear unbiased estimator when the errors are uncorrelated, have constant variance, and are centred, but the estimate is a sum of squared deviations, so a single observation with a large residual contributes quadratically and can move the fit arbitrarily far. That fragility is a property of the loss, not of the projection, and the alternatives are compared in squared, absolute, and quantile loss. The projection derivation survives only for the squared loss, which is why the closed form does not generalise to classification: replacing the squared error with a likelihood turns the same problem into one with no linear-algebra solution, as logistic regression as conditional likelihood shows.