The directions that matter

Karl Pearson, in 1901, was after a stubbornly practical thing: given a swarm of measurements — heights against weights, exam scores against each other, whatever a Victorian biometrician could tabulate — what is the best-fitting line through the cloud? Not the regression line, which minimizes vertical error and secretly treats one variable as the boss, but a line that treats every axis as an equal citizen and minimizes the perpendicular distance to the points. He found it, and he noticed something more: the same construction hands you not just one line but a whole set of perpendicular axes, ranked by how spread out the data is along each. Three decades later Harold Hotelling, in 1933, gave the method its modern name and its statistical spine — principal component analysis. And here is the quiet miracle this essay is about: Pearson's best-fit axes, Hotelling's directions of maximum variance, the eigenvectors of the covariance matrix, and the singular vectors of the data are all the same objects, arrived at from four different doors. The whole arc of the last six essays — eigenvectors and the spectral theorem, positive-definite matrices, the singular value decomposition — converges here, and it all cashes out at once.

This is the capstone, so let me say up front what it pays off. PCA is the statement that a cloud of data, however high-dimensional and tangled, has natural axes — perpendicular directions, ordered by importance — and that keeping only the first few of them throws away almost nothing. That claim rests on three results you've already met, and the joy of this essay is watching them lock together.

First, center the cloud

Everything begins with a step so humble it's easy to skip. Stack your data as a matrix X with one row per observation and one column per feature — n points in \mathbb{R}^d. Before doing anything else, subtract the mean of each column, so the cloud is recentered on the origin:

\tilde{x}_i = x_i - \bar{x}, \qquad \bar{x} = \frac{1}{n}\sum_{i=1}^n x_i.

Why insist on this? Because PCA is a theory about variation, about how the points spread away from their own center of mass — and spread is measured from the middle, not from some arbitrary origin the coordinate system happened to inherit. An uncentered cloud sitting far from zero has enormous "variance" that is really just its distance from the origin, a fact about where you parked the axes, not about the shape of the data. Center first, and every number that follows is honestly about shape. From here on, assume X is already centered.

The direction of maximum variance

Now the question, made precise. Take any unit direction w (so \lVert w \rVert = 1) and project every point onto the line it spans: the projection of \tilde{x}_i is the scalar w^\top \tilde{x}_i, its signed distance along w. How spread out are those projections? Their variance, since the data is centered, is

\operatorname{Var}(w) = \frac{1}{n}\sum_{i=1}^n \left(w^\top \tilde{x}_i\right)^2 = \frac{1}{n}\, w^\top X^\top X\, w = w^\top C\, w, \qquad C = \frac{1}{n}X^\top X.

That matrix C is the covariance matrix of the data, and everything about PCA is written in it. The variance along direction w is exactly the quadratic form w^\top C w — and C is no ordinary matrix. It's symmetric by construction (X^\top X always is), and it's positive semidefinite, because w^\top C w = \tfrac1n\lVert Xw\rVert^2 \ge 0 is an average of squares and a variance can't be negative. Those two properties are not decoration; they are the licence for everything that follows.

The first principal component is the direction of most spread — the unit w that maximizes w^\top C w. Written as an optimization:

w_1 = \arg\max_{\lVert w \rVert = 1} w^\top C\, w.

Maximizing a quadratic form over the unit sphere is a problem we have already solved, twice over. Because C is symmetric, the spectral theorem hands it an orthonormal eigenbasis C = Q\Lambda Q^\top with real eigenvalues \lambda_1 \ge \lambda_2 \ge \cdots \ge 0 (nonnegative because C is positive semidefinite). Rewrite w in that eigenbasis as y = Q^\top w, and the quadratic form becomes a plain weighted sum of squares, w^\top C w = \sum_i \lambda_i y_i^2, with \sum_i y_i^2 = 1. You're spending a fixed budget of one unit across the coordinates, and each unit of budget buys \lambda_i of variance — so of course you pour it all onto the largest eigenvalue. The maximum is \lambda_1, achieved by pointing w straight along the top eigenvector.

\boxed{\;w_1 = \text{top eigenvector of } C = \tfrac1n X^\top X, \qquad \operatorname{Var}(w_1) = \lambda_1.\;}

The direction of maximum variance is the top eigenvector of the covariance matrix, and the variance it captures is the eigenvalue. That is the entire idea of PCA in one line.

Principal components as orthogonal axes

One direction is a beginning, not an answer. Where does the second most important direction live? Ask for the unit vector that maximizes variance among all directions perpendicular to w_1 — because the spread already explained by w_1 shouldn't be double-counted. The spectral theorem answers before you finish the question: with the top eigen-direction ruled out, the same budget argument pours everything onto the next eigenvalue, so w_2 is the second eigenvector, capturing variance \lambda_2. Continue, and the principal components are simply the eigenvectors of C in descending order of eigenvalue — w_1, w_2, \ldots, w_d — a full set of orthonormal axes for the data, each one perpendicular to all the rest and each capturing strictly less variance than the last. This is Pearson's set of best-fit axes and Cauchy's principal axes of a quadratic surface, the same perpendicular frame the spectral theorem promises any symmetric matrix. The covariance ellipsoid — the cloud's shape, idealized — has these eigenvectors as its axes and \sqrt{\lambda_i} as its semi-axis lengths. PCA is nothing more exotic than finding the axes of the ellipse the data is shaped like.

The SVD says the same thing, without ever forming C

There's a second door to exactly these axes, and it's the one you'd actually use in practice. Take the singular value decomposition of the centered data itself, X = U\Sigma V^\top. Recall the identity that powers the whole SVD: X^\top X = V\,(\Sigma^\top\Sigma)\,V^\top. Divide by n and compare to the eigendecomposition C = Q\Lambda Q^\top we just used:

C = \frac{1}{n}X^\top X = V\left(\frac{\Sigma^\top\Sigma}{n}\right)V^\top.

Read off the correspondence. The right singular vectors V of X are the eigenvectors Q of C — they are the principal components. And the eigenvalues are the squared singular values, scaled: \lambda_i = \sigma_i^2 / n, so the variance along the i-th principal component is \sigma_i^2/n and the standard deviation is \sigma_i/\sqrt{n}.

\boxed{\;\text{principal components of } X \;=\; \text{right singular vectors } v_i \;=\; \text{eigenvectors of } C.\;}

This is not a coincidence to memorize; it's the same fact seen twice. Forming C = \tfrac1n X^\top X and diagonalizing it, versus taking the SVD of X directly, are two routes to one destination — and the SVD route is numerically the sane one, because squaring the data into X^\top X squares its condition number and drowns the small singular values in roundoff. Either way, the arc closes: the SVD essay called PCA "the SVD wearing a statistician's hat," and now you can see the hat is just the factor of 1/n and the word variance.

Projection is dimensionality reduction

Here is where the ranking earns its keep. You have d orthonormal axes, sorted by how much variance each carries. To compress the data, keep the top k and drop the rest. Projecting a centered point onto the span of the first k principal components gives its coordinates in that reduced space:

z_i = W_k^\top \tilde{x}_i, \qquad W_k = [\,w_1 \mid w_2 \mid \cdots \mid w_k\,],

a k-dimensional summary of a d-dimensional point. When k is far smaller than d — a hundred features crushed to two or three — you've done dimensionality reduction: a smaller representation that keeps the loud directions and discards the quiet ones. To go back, you reconstruct by re-expanding in the same axes (and adding the mean back): \hat{x}_i = W_k z_i + \bar{x}. The reconstruction lands on the k-dimensional flat spanned by the top components — the best-fitting k-plane through the cloud, Pearson's line generalized to a plane and beyond.

Why the dropped directions are safe: Eckart–Young

Truncating begs the obvious worry — how much did we lose? The answer is precise and beautiful, and it's an old friend from the SVD essay. The total variance in the data is the sum of all the eigenvalues, \sum_i \lambda_i = \tfrac1n\sum_i \sigma_i^2 (the trace of C, invariant under any rotation of axes). Keep the top k and the variance you retain is \sum_{i\le k}\lambda_i; the variance you sacrifice, the mean-squared reconstruction error, is exactly the tail you dropped:

\frac{1}{n}\sum_i \lVert \tilde{x}_i - \hat{x}_i \rVert^2 = \sum_{i > k} \lambda_i = \frac{1}{n}\sum_{i>k}\sigma_i^2.

And that this is the smallest possible error for any k-dimensional projection is not a hope — it's the Eckart–Young theorem. The rank-k truncation of the SVD is the best rank-k approximation of X there is, so no other set of k axes, perpendicular or not, can reconstruct the cloud more faithfully. PCA doesn't find a good low-dimensional summary; it finds the optimal one, and the price of admission is precisely the singular values you chose to leave behind. If those trailing \sigma_i are small — and for real, correlated data they almost always decay fast — the leftover directions genuinely carry almost none of the story, and dropping them costs almost nothing. That is why the quiet directions are safe to drop: Eckart–Young guarantees the cost is exactly the variance they held, and by construction they held the least.

Explained variance: how much story each axis carries

The natural way to report all this is a single ratio. The fraction of the total variance living along the i-th principal component is its explained variance:

\text{explained variance}_i = \frac{\lambda_i}{\sum_j \lambda_j} = \frac{\sigma_i^2}{\sum_j \sigma_j^2}.

It's a proportion between 0 and 1, and the components sum to 100% of the story. For an elongated, cigar-shaped cloud the first component might explain 95% on its own — a two-dimensional dataset that is really one-dimensional wearing a disguise — and the standard practice of "keep enough components to explain 90% (or 99%) of the variance" is just reading a running total of these ratios until it crosses your threshold. The cumulative explained variance is the dial that turns "how much can I compress?" into a number.

Where this goes — and where it came from

Principal component analysis is the place the whole foundations arc pays off at once. Center the data, form the positive-semidefinite covariance matrix, and the spectral theorem hands you orthonormal axes ranked by variance; the SVD hands you the very same axes without ever squaring the data, and Eckart–Young certifies that keeping the top few is the best compression that exists. Symmetric-means-orthogonal-axes, positive-semidefinite-means-nonnegative-variances, and best-low-rank-means-safe-to-truncate — three theorems, one method. From here the road runs outward into everything that reaches for structure in high-dimensional data: whitening and preprocessing, latent-factor models, the eigenfaces and the recommender systems, and the long tail of nonlinear cousins that keep PCA's ambition while trading its linearity. But the linear core, the thing all of them are bending, is the humble claim you can now derive from scratch: data has natural axes, and most of them don't matter.

The fastest way to feel it is to reshape a cloud and watch the axes chase it. The PCA playground scatters a 2-D point cloud, computes the 2\times2 covariance and its two eigenvectors analytically, and draws the principal axes straight through the centroid — PC1 along the long direction, PC2 perpendicular, each drawn at two standard deviations for visibility. Regenerate the cloud and watch the axes swing to follow the new shape; the readout shows the explained variance each axis carries. Then hit project and watch the second dimension vanish — every point collapsing onto the PC1 line, the best one-dimensional shadow the cloud can cast, with exactly the reconstruction error Eckart–Young predicted. It's one thing to read that a cloud of data has directions that matter more than others. It's another to squash the cloud flat and watch how little you lose.