Projections and orthogonality

Two of the most-used procedures in data science — fitting a linear model and reducing dimensions — are the same geometric act pointed in different directions: drop perpendiculars onto a subspace. This lesson builds that geometry once, carefully, and then collects both procedures (and Gram–Schmidt, and the Fourier idea of coordinates-in-a-basis) as corollaries.

Projection onto a line, then a subspace

Given a unit vector $u$ and any point $x$, the multiple of $u$ closest to $x$ is found by demanding the error be perpendicular:

The scalar $u^\top x$ is the coordinate of $x$ along $u$ — the inner product from lesson 1, now wearing its second hat: not just similarity, but how much of this direction does the data contain. For a subspace spanned by the columns of $W$ (full column rank), the same perpendicularity demand $W^\top(x - W\hat c) = 0$ solves to

the projection matrix: symmetric, idempotent ($P^2 = P$ — projecting twice changes nothing), with eigenvalues only 0 and 1 (a direction is either kept or killed — nothing in between, which is a nice sanity-check application of the eigen-lesson). If the columns are orthonormal, $W^\top W = I$ and the formula collapses to $P = WW^\top$: project = read off coordinates ($W^\top x$), rebuild ($W \cdot$). That collapse is why orthonormal bases are the luxury goods of numerical computing, and manufacturing them from arbitrary independent columns is Gram–Schmidt (subtract from each new column its projections onto the previous ones, normalize; done stably, that is the QR factorization your lstsq actually runs).

Pythagoras, the accountant of variance

Because $\hat x \perp (x - \hat x)$, lengths obey

— kept plus lost equals total, exactly. Summed over a centred dataset this becomes the identity that runs half of statistics: variance captured by a direction + mean squared residual = total variance. Maximizing one is minimizing the other; they are a single optimization wearing two costumes. Which resolves a question you can now feel in your hands:

Corollary one: least squares was a projection all along

The Supervised Learning module’s linear regression minimized $\lVert y - X\beta \rVert^2$ — the distance from $y$ to the column space of $X$. The minimizer is the projection $\hat y = P y$ with $P = X(X^\top X)^{-1}X^\top$ (the hat matrix, whose idempotence you proved as an exercise there without being told why it held: now you know — it is a projection, and projections are idempotent by meaning). The normal equations $X^\top(y - X\hat\beta) = 0$ are nothing but the perpendicularity condition; “residuals are orthogonal to every regressor” is the geometric definition of best fit, not a happy accident.

Corollary two: PCA, stated properly

For centred data with covariance $\Sigma$, the top-$k$ principal directions are the orthonormal $u_1, \ldots, u_k$ maximizing captured variance $\sum_j u_j^\top \Sigma u_j$ — equivalently (Pythagoras) minimizing reconstruction error. The solution is the top-$k$ eigenvectors of $\Sigma$, and the projection $W W^\top x$ with $W = [u_1 \cdots u_k]$ is exactly the truncated-SVD reconstruction of the data matrix: PCA = SVD = the widget’s snap button, one object seen from three sides. Two practitioner’s footnotes with teeth: centre first (uncentred “PCA” mostly finds the mean’s direction), and scale matters (variance is unit-dependent — a feature measured in grams will dominate one in kilograms 10⁶-fold; standardize unless the units are genuinely comparable, the same lecture the Regularization lesson gave for penalties).

The residual as signal

One habit separates people who use projections from people who understand them: looking at what the projection throws away. In regression, residual structure is unmodeled signal (the diagnostic argument from the Linear Regression lesson). In PCA, the discarded tail is either noise (good riddance — the denoising story from the SVD lesson) or the anomaly that didn’t fit the main structure — reconstruction error is a workhorse anomaly detector for exactly this reason. And in Gram–Schmidt, the residual is the algorithm: each new basis vector is literally “what’s left after projecting out everything already explained” — which is also, word for word, what a regression coefficient means in the presence of other regressors (the Frisch–Waugh view, for the econometricians in the room).

import numpy as np

rng = np.random.default_rng(0)
X = rng.normal(size=(80, 2)) @ np.array([[1.35, 0.6], [0.0, 0.45]])  # correlated cloud
X -= X.mean(axis=0)                                # centre FIRST

# PCA three ways — identical answers
Sigma = X.T @ X / len(X)
eigval, eigvec = np.linalg.eigh(Sigma)             # 1: eigen of covariance
U, s, Vt = np.linalg.svd(X, full_matrices=False)   # 2: SVD of data
u = Vt[0]                                          # top principal direction
np.testing.assert_allclose(abs(u @ eigvec[:, -1]), 1, atol=1e-10)

P = np.outer(u, u)                                 # 3: projector onto it
kept = np.var(X @ u)
lost = np.mean(np.sum((X - X @ P) ** 2, axis=1))
np.testing.assert_allclose(kept + lost, X.var(axis=0).sum())   # Pythagoras

Exercises