Principal Component Analysis

Principal component analysis reduces the dimension of a dataset by projecting it onto the directions of greatest variance. It centers the data, builds the covariance matrix, and takes its eigenvectors: the first principal component is the eigenvector with the largest eigenvalue, and that eigenvalue is exactly the variance captured along it.

Here is the fact that makes the whole method click, and it is not the one most write-ups lead with. Center your data and pick any direction at all. Every point splits into a piece along that direction and a piece perpendicular to it, and by Pythagoras those two variances add to the same total no matter which direction you picked. There is one budget. Choosing a direction only chooses where to cut it. So "find the direction of greatest spread" and "find the direction of least reconstruction error" are not two goals that happen to agree — they are one goal, said twice.

How Principal Component Analysis Works

Start with why centering is not bookkeeping. Twelve rocks, each measured twice, lengths and widths between 2 and 11 centimetres. Drawn on axes through the origin, that cloud sits up and to the right, nowhere near (0, 0). Square the projections onto a line through the origin and add them up, and the total mixes the cloud's spread with how far its mean sits from the origin, which is an accident of where you happened to put your zero. Variance itself does not care where zero is; the origin-anchored sum of squares does. Subtract the mean from every row and the cloud slides onto the origin. The shape is untouched; only the vantage point moved. From here the sum of squared projections onto a direction through the origin, divided by n − 1, is the variance along it, and the question becomes answerable.

Now fix a unit direction u. Project each centered point onto it: the projection is the foot of a perpendicular dropped from the point to the line, and the coordinate along the line is the dot product x·u. The variance of those twelve coordinates is what u captures. Call it v(u). The leftover — how far each point still sits off the line, squared and summed with the same n − 1 denominator as the sample covariance — is r(u), the reconstruction error you would pay if you kept only the coordinate along u and threw the rest away.

Each point is the hypotenuse of a right triangle whose legs are those two pieces, so v(u) + r(u) = trace(C) for every u, where trace(C) is the total variance of the data and does not depend on u at all. On these twelve rocks that total is 14.55 cm². Every direction you can name splits 14.55 into two numbers that add back to 14.55. Maximizing the first is the same act as minimizing the second.

Worth pausing on the word perpendicular, because the neighbouring lesson trains the opposite reflex. Linear regression measures its residual vertically, y − ŷ, because it is predicting y from x and the miss it cares about is a miss in y. PCA predicts nothing. Both coordinates are just measurements, neither is the answer, so the distance that matters is the honest one: the shortest distance from the point to the line.

That leaves the search. You could sweep the angle and watch v(θ) rise and fall, and the lesson does exactly that so you can see the shape of it. But you never have to. All the information about spread in any direction is already packed into a 2 × 2 matrix, the sample covariance C = XᵀX / (n − 1), whose diagonal holds each feature's own variance and whose off-diagonal holds how the two move together. For any direction, v(u) = uᵀCu. Maximizing that over unit vectors is a textbook eigenvalue problem, and its answer is a formula, not a search: the maximizing direction is the eigenvector of C with the larger eigenvalue, and the eigenvalue is the variance captured along it. Two directions come out, perpendicular to each other, and their eigenvalues sum to the trace — the same fixed budget, now cut at its most lopsided point.

One naming note for when you write the code. Everyone writes the first component as u₁ or PC1, counting from one. np.linalg.eigh returns eigenvalues in ascending order, so without a re-sort that same vector is its last column, vecs[:, -1].

Principal Component Analysis Step by Step

The default run: twelve rocks, integer centimetres, with a trial direction of θ = 0°. That angle is not a neutral placeholder. Keeping only the first coordinate is the direction θ = 0°, and dropping a column is what almost everyone does before they meet PCA. The run prices that move first, then beats it.

1. The raw data. Twelve points, all in the upper-right quadrant, with the mean marked at (6, 5). The origin sits off in empty space with no data near it.
2. Center. Every row loses (6, 5). The cloud slides onto the origin and the frame closes in. The rock measured (6, 5) lands exactly on the origin and stays there for the rest of the lesson — it is the mean, so it has nothing to say about spread in any direction.
3. Price θ = 0°. Twelve feet drop onto the horizontal axis, twelve dashed segments show the perpendicular gaps. Keeping the first measurement captures 9.27 cm² of spread. What it misses is on screen but not yet counted.
4. The budget. Eleven rocks fade and one right triangle is left standing — the rock whose triangle has the largest area at this angle. The rock measured (11, 9) centers to (5, 4); its coordinate along the line is z = 5, its perpendicular miss is r = 4, and the rock itself is the hypotenuse, ‖x‖² = 25 + 16 = 41. Rotate the line and z and r both change — the hypotenuse cannot. Sum that over all twelve and divide by n − 1: 9.27 + 5.27 ≈ 14.55 cm² (each value is rounded on its own; exactly, 102/11 + 58/11 = 160/11). The budget bar appears with its outer ends pinned, and they never move again.
5. Sweep, 20° at a time. The line rotates through nine directions and each one commits a dot to the spread(θ) panel. The split slides back and forth inside a bar of constant width, which is the identity made visible: 88.25% at 20°, a peak of 94.86% at 40°, then down through 80.47% at 60°, 51.83% at 80°, and a floor of 5.78% at 120° — the direction that captures almost nothing, which is to say the direction perpendicular to where the data actually lives. Back up to 32.81% at 160° and the curve is one full hump over 180°.
6. Covariance. The line and the feet leave. C = [[9.27, 6.27], [6.27, 5.27]] — four numbers, computed in one pass over the data. C₁₁ = 9.27 is the θ = 0° sample already sitting on the curve.
7. The curve. The exact curve v(θ) = trace(C)/2 + d·cos(2(θ − θ₁)) is drawn behind the nine dots, and every dot is sitting on it. This is the beat that makes the matrix worth building: v(u) = uᵀCu returns any direction's captured variance in closed form, so those four numbers already contain every angle the sweep tried and every angle it skipped.
8. Eigen. The maximum of that curve is at θ₁ = 36.16°, capturing λ₁ = 13.86 of 14.55, or 95.26%, with λ₂ = 0.69 left over — and it is found by eigendecomposing C, not by scanning the curve. Look at the sweep again: the nine samples bracketed 36.16° and none of them hit it. The best guess, 40°, came 0.4 percentage points short. That gap is the argument for the closed form.
9. Project. The world rotates by −36.16° so u₁ lies horizontal. Same twelve rocks, same distances between them, new addresses: each row is now (z₁, z₂), how far along the principal direction and how far off it.
10. Reduce. Every z₂ goes to zero and the cloud collapses onto the z₁ axis. Twelve rocks now cost one number each instead of two, and 95.26% of the original variance is still in them. The rock at the origin does not move; it had no z₂ to lose.
11. Reconstruct. Rotate back and rebuild each rock as x̂ = z₁·u₁ + mean. Every rebuilt rock lands on one line, and the gap between it and the real one is precisely the coordinate that was deleted. Sample-normalized residual scale: √λ₂ = 0.83 cm. The per-point RMSE uses n instead of n − 1 and is √((n − 1)λ₂/n).
12. Compare. Two pictures, same storage cost of one number per rock. Drop the second column and the typical miss is 2.30 cm. Drop the second component and it is 0.83 cm — 2.8 times closer, for free, because the second option got to choose its direction instead of inheriting one from whoever named the columns.

What PCA Is Actually For

Four jobs, and the twelve rocks are a scale model of all of them. Compression: store k numbers per row instead of d and accept a known reconstruction error — the trade the last two beats price out. Visualization: project a hundred features onto the top two and plot them, which is one way to inspect high-dimensional data. Decorrelation: the new coordinates are uncorrelated by construction, which is what you want in front of a model that chokes on collinear inputs. Denoising: if the signal is high-variance and the noise is not, dropping the tail components drops mostly noise — a big if, and the second failure case below is what happens when it is false.

The reflex that goes with all four: fit on the training set only. The mean, the per-feature standard deviations, and the components themselves are all estimated from data, so estimating them on the full dataset and then splitting leaks test information into training. Fit StandardScaler and PCA on the training fold, apply both to the test fold, and treat k as a hyperparameter like any other.

Why the Explained Variance Ratio Can Lie

λ₁ / (λ₁ + λ₂) is the number every PCA implementation prints, and on the rocks it says 95.26%. Turn one knob and watch it betray you.

Record width in millimetres instead of centimetres. Nothing about the rocks changed; the second column is now ten times larger. The covariance matrix is dominated by that column, PC1 swings from 36.16° to 83.19°, and the reported explained variance rises to 99.67%. A better score for a worse answer: "the first principal component" has quietly become "width, in millimetres." PCA maximizes variance, variance carries units, and the feature with the bigger numbers wins by default.

The fix is to scale each feature to unit variance before running PCA, which turns the covariance matrix into the correlation matrix. Do that and the answer stops depending on the ruler: 45.00° and 94.85%, at every multiplier you can pick. That is the check worth remembering — not that standardizing produces a nicer number, but that afterwards the scale knob provably does nothing at all.

When Variance Isn't Information

The second failure mode cannot be fixed by the same variance-only objective, and it is the more important one. Take twelve samples that genuinely come from two groups, arranged so both groups run along the same broad trend and sit slightly offset from each other. PC1 takes 86.91% of the variance, and along it the two groups interleave completely — their projections overlap end to end. PC2, the 13.09% that PCA tells you to discard, separates them with a clean gap and no overlap at all.

Standardizing does not rescue this (83.23%, still fused). Nothing does, because nothing is broken. PCA is unsupervised: it never sees a label, so it cannot possibly optimize for one. It answers "which direction do these numbers vary along most", and you were hoping it would answer "which direction tells my classes apart". Those coincide often enough to be dangerous and diverge often enough to matter. When you have labels and want separation, ask for separation — linear discriminant analysis maximizes between-class spread relative to within-class spread, which is a different question with a different answer.

Time and Space Complexity of PCA

Forming the covariance matrix from n points with d features costs O(n·d²) time; eigendecomposing a d × d symmetric matrix costs O(d³). Projecting all n points onto the top k components costs another O(n·d·k). Space is O(d²) for the matrix plus O(n·d) for the data, which is why the time complexity is usually dominated by whichever of n and d is enormous. With d = 2 the eigendecomposition is a closed-form expression you could evaluate by hand; with d in the thousands the d³ term is what pushes real implementations toward the SVD of the centered matrix instead, which gets the same components with better numerics.

Edge Cases

• A point exactly at the mean. It centers to the origin, both principal coordinates are zero, and it survives any reduction untouched. The default data contains one on purpose.
• A perfectly round cloud (λ₁ = λ₂). The spread curve is flat and every direction captures the same amount, so there is no principal direction. The eigenvectors are not unique and any orthogonal pair is a valid answer. An implementation that reports "the" first component here is reporting an arbitrary choice made by its numerical routine.
• A one-dimensional cloud (λ₂ = 0). All points sit exactly on a line. The reduction to one component is lossless, the explained variance is 100%, and the reconstruction is exact.
• The sign of a component is arbitrary. u₁ and −u₁ are both eigenvectors and capture identical variance. Different libraries — different versions of the same library — will hand you opposite signs on the same data. Never read meaning into which way a component points; this lesson pins the angle to [0°, 180°) so the picture is stable.
• Fewer points than features. With n samples in d dimensions, at most n − 1 eigenvalues can be nonzero. Ask for more components than that and you get numerical noise dressed as structure.

Common Mistakes

• Reporting a variance as a distance. λ₂ = 0.69 is in squared units. The sample-normalized residual scale is √λ₂ = 0.83 cm; the per-point root-mean-square distance is √((n − 1)λ₂/n). Skipping the square root reports cm² as if it were cm: it overstates any error above 1 and understates any below it.
• Skipping the centering. Without it, the first component tends to point at the data's mean rather than along its spread, because the mean's own distance from the origin swamps everything else. Common PCA implementations center by default, but related decompositions may not, which is exactly why it is easy to forget it is happening.
• Running PCA on raw features with different units. Millimetres beat centimetres, dollars beat percentages, and the component you get is a statement about your unit conventions. Scale to unit variance unless the features are genuinely commensurable.
• Reading a high explained-variance ratio as a good representation. The ratio measures how concentrated the variance is, not how useful the result is. The millimetres run scores 99.67%, yet mapped back to centimetres its direction keeps less of the rocks' spread than the 95.26% run does.
• Treating the components as interpretable features. A principal component is a weighted mix of the original ones, chosen for a geometric property, and there is no reason it should correspond to anything you have a word for. Rotating it toward interpretability is a separate technique with separate assumptions.
• Using PCA where the goal is class separation. See above. The direction PCA discards can be the only one that matters.

A Note on Simplification

What is animated is real PCA, not a reduced version of it — but it runs on two features, which is the smallest case where "reduce the dimension" means anything at all. Real datasets carry tens to thousands of features and keep the top k, where choosing k is its own question, usually answered with a scree plot or a cumulative-variance threshold, and only interesting above two dimensions.

Production implementations often take a different route to the same answer. Instead of forming the covariance matrix and eigendecomposing it, many take the singular value decomposition of the centered data matrix directly; some, such as scikit-learn's covariance_eigh solver, still form C when there are few features and many rows. The components are identical; the SVD is better conditioned, because squaring the data to build C squares its condition number along with it. The lesson solves the 2×2 covariance matrix in closed form because it makes the two-variances-and-a-covariance story visible, and this small well-conditioned example makes every method agree to the displayed precision.

Not covered here: whitening (dividing each component by its standard deviation after projection), kernel PCA (running the same argument in a nonlinear feature space), and the incremental and randomized variants that make PCA tractable when the data does not fit in memory.