Linear Regression as Projection

A model with 40 features fits a training set of 40 points perfectly, residuals exactly zero, and predicts nonsense on anything new. Nothing went wrong in the optimiser. The design matrix spans the entire 40-dimensional space the targets live in, so projecting onto that space leaves nothing behind. The geometry of least squares tells you this before you run anything.

Given \(X \in \mathbb{R}^{n \times p}\) and targets \(y \in \mathbb{R}^n\), least squares seeks \(\hat{\beta}\) minimising \(\|y - X\beta\|_2^2\). The set \(\{X\beta : \beta \in \mathbb{R}^p\}\) is a subspace of \(\mathbb{R}^n\): the column space of \(X\). Minimising Euclidean distance from \(y\) to a subspace is orthogonal projection, and everything else follows.

The normal equations fall out of orthogonality

The residual \(y - X\hat{\beta}\) must be orthogonal to every column of \(X\), otherwise you could move within the subspace and get closer. Writing that condition:

The fitted values are \(\hat{y} = X(X^\top X)^{-1}X^\top y = Hy\), where \(H\) is the hat matrix. \(H\) is symmetric and idempotent (\(H^2 = H\)), which is the algebraic signature of a projection: project twice and nothing changes. Its trace equals \(p\), the number of parameters, and \(\mathrm{tr}(H)\) is the standard definition of effective degrees of freedom, which generalises usefully to ridge and smoothers where "number of parameters" stops being an integer.

The diagonal entries \(h_{ii}\) are leverage: how much observation \(i\) pulls its own fitted value. Since \(\sum_i h_{ii} = p\), average leverage is \(p/n\), and any point several times that is geometrically extreme in feature space regardless of whether its residual looks large.

Why nobody forms \((X^\top X)^{-1}\)

The normal equations are correct mathematics and poor numerics. Forming \(X^\top X\) squares the condition number: if \(X\) has \(\kappa(X) = 10^6\), then \(\kappa(X^\top X) = 10^{12}\), and in double precision you have thrown away most of your significant digits before the solve begins.

Production solvers use a QR decomposition \(X = QR\) with \(Q\) orthonormal and \(R\) upper triangular, reducing the problem to the triangular system \(R\hat{\beta} = Q^\top y\) with condition number \(\kappa(X)\) rather than its square. For rank-deficient or nearly rank-deficient problems, the SVD gives the minimum-norm solution through the pseudoinverse and lets you inspect the singular values to see exactly how deficient the design is. This is the same reason scipy.linalg.lstsq exists and inv(X.T @ X) @ X.T @ y is a code smell.

What the Gauss-Markov theorem does and does not promise

Under the assumptions that errors have zero mean, constant variance, and are uncorrelated, OLS is the best linear unbiased estimator: minimum variance among all linear unbiased estimators. Both qualifiers carry weight. Ridge regression is biased and routinely achieves lower MSE, which contradicts nothing. Normality of errors is not required for Gauss-Markov; it is required only for the exact t and F distributions of the test statistics, and the central limit theorem covers that in large samples.

The assumption that fails most often in practice is uncorrelated errors with constant variance. Heteroscedasticity leaves \(\hat{\beta}\) unbiased but makes the usual standard errors wrong, which is what robust (Huber-White) standard errors repair without touching the coefficients.

When it breaks

Collinearity inflates variance without inflating error. When two features are nearly collinear, \(X^\top X\) is nearly singular and the coefficient estimates have enormous variance and unstable signs, while the predictions remain fine because the projection is onto a well-defined subspace. This is precisely the case where interpreting coefficients is dangerous and predicting is safe, and confusing the two is the most common analysis error in applied regression.

\(p \geq n\) has infinitely many exact solutions. The column space fills \(\mathbb{R}^n\), the residual is zero, and \(X^\top X\) is singular. Least squares stops having a unique answer, and something must break the tie: minimum-norm via the pseudoinverse, or an explicit penalty. Modern overparameterised models live here permanently, which is why the implicit bias of the optimiser is a real research subject rather than a footnote.

Leverage points bend the fit and hide the evidence. A high-leverage observation pulls its fitted value toward itself, so its residual can look small precisely because it is influential. Reading residuals alone will not find it; Cook's distance, which combines leverage and residual, will.

The projection view assumes the model is linear in parameters, not in features. Polynomial and spline bases stay inside this whole framework. What breaks it is a genuinely nonlinear parameterisation, which is where the closed form disappears and iterative optimisation takes over.