Principal Component Analysis

Start with two features and eleven people. Shorter people tend to be lighter, taller people heavier, and around that overall relationship there is variation, someone heavy for their height, someone light. Each feature has its own number line, and the scatter is what you get by pairing them.

So does each feature contribute entirely new information? Not entirely. Much of what the weight line says, the height line already said, and a single combined measurement of how big someone is would describe most of the shared pattern. Not all of it, since being heavy for your height is real and is not on that combined axis, which is the part that is unique to each feature.

The idea grows with the number of columns. Measure sixty people five ways, height, arm span, leg length, sitting height and shoe size, and the five are genuinely different measurements that all follow one shared pattern, overall body size, each with something of its own besides.

…

The original features are the axes we happened to measure along, a height axis and a weight axis. A direction is something else, a weighted combination of the two. Height and weight increasing together is one direction; height increasing while weight decreases is another; any weighting at all is a third. Turn the arrow below through the crowd and it combines the two measurements into one new axis at every angle.

…

The first move is to find the middle of the cloud. For two features it is the pair of means, and for the four people the page works by hand it is a clean number.

mean point = (mean height, mean weight) = (170, 68)

Every person is then some horizontal distance and some vertical distance from that point, which the next widget draws as two dashed legs. The mean is the coordinate-wise centre of the observations, and it is the pivot everything that follows turns about.

Subtract the mean point from every person and each one becomes a deviation, the pair of legs with the pivot moved to the origin.

dᵢ = xᵢ − x̄

…

Why bother? Because the method is about to measure spread along directions, and without centring, a direction can look important just because the whole cloud sits far from the origin. Untick the switch below and the scores along the arrow stop being deviations and start being positions, and the mean of their squares balloons from a measure of spread into a measure of where the cloud is.

…

Choose a unit direction u. For each centred person, drop a perpendicular onto the line through the origin in that direction, mark where it lands, and read off the signed distance of the landing point from the origin. That one number is the person’s score along u, and it is a dot product.

scoreᵢ = dᵢ · u

…

Two things about the scores. First, u has to be a unit vector. Let the arrow grow and every score grows with it, so the apparent spread along a direction could be made as large as you like by stretching the arrow rather than by choosing a better direction. Fixing u · u = 1 makes the comparison between directions a comparison between directions.

…

Second, because the cloud is centred, the scores are centred too, and their variance is simply the mean of their squares.

variance along u = mean of (dᵢ · u)²

Now sweep. Some directions squeeze the people together on the number line and some spread them out, and the curve beneath the cloud records the variance at every angle of the half turn. It has one peak and one trough, a quarter turn apart.

…

That widest direction has a name. The first principal component is the unit direction along which the centred observations have the greatest variance, and on the measured four it sits at 45 degrees, where the variance of the scores is 100 and every other angle gives less. On the crowd it sits at 42.2 degrees. The playground at the top of the page draws it in indigo, and it is not the height axis and not the weight axis. It is the cloud’s own.

At 45 degrees the direction is u₁ = (0.71, 0.71), and a person’s score along it is a weighted sum of their centred features.

score = 0.71 × centred height + 0.71 × centred weight

Both features contribute positively, and the equal coefficients mean a unit of movement in either contributes the same, after centring. The coefficients that define the direction are called the loadings, and the widget builds one person’s score from them, the height part in indigo, the weight part in amber, and their sum.

…

height contribution = 10 / √2 ≈ 7.071
weight contribution = 10 / √2 ≈ 7.071
person 1 score = 20 / √2 ≈ 14.142

person 3 score = 5 / √2 − 5 / √2 = 0

The direction (0.71, 0.71) and the direction (−0.71, −0.71) describe the same axis, and the method has no reason to prefer one. Flipping the sign reverses every score numerically and changes nothing else. Not the distances between projected people, not the variance, not the reconstructions.

…

Which leads to a caution about names. The first component on the crowd resembles overall body size, and it would be natural to call it that, but the method never learned the phrase. It found a mathematical direction. Whether that direction deserves a human label depends on the loadings and on knowing something about bodies, and a component can combine features in ways with no clean real-world label at all. When it does, the honest move is to leave it unnamed.

Project everyone onto the first component and each person leaves a stub behind, the perpendicular from where they are to where they landed. The green landing points are the variation the first component keeps; the rose stubs are the variation it does not, and the readouts total both.

…

The second component is chosen the same way with one extra rule. It must be perpendicular to the first, and among the perpendicular directions it must carry as much of the remaining variance as possible. In two dimensions there is only one perpendicular axis left, apart from its sign, so the second component is forced. In higher dimensions the search continues inside the space perpendicular to every component found so far.

The components are ordered by their variances, largest first, so the first captures the most, the second the most of what remains, and so on down. Ties are possible. When two eigenvalues are equal, any direction inside the plane they span carries the same variance and the component directions inside that plane are not uniquely determined.

λ₁ ≥ λ₂ ≥ … ≥ λₚ

Everything so far was done by rotating and measuring. The covariance matrix is what lets the answer be found without the rotation. For centred height and weight it holds three quantities the statistics primer already defined, the variance of each feature and the covariance of the pair, arranged so that each feature’s spread sits on the diagonal and their paired movement sits off it. The matrix is symmetric, because the covariance of height with weight is the covariance of weight with height.

C = [ Var(height)          Cov(height, weight) ]
    [ Cov(height, weight)  Var(weight)         ]

…

A matrix is also something that acts on a direction. Pick a unit vector u and compute Cu. For most directions Cu points somewhere else, turned away from u. For a few special directions Cu stays on the same axis as u, only longer or shorter, and those directions satisfy the eigenvector equation.

C u = λ u

…

at 0°:  u = (1, 0)
        Cu = (62.5, 37.5)

at 45°: u = (1/√2, 1/√2)
        Cu = 100u ≈ (70.7, 70.7)

Now the two halves meet. The direction the variance sweep found by rotating, 45 degrees on the measured four, is the direction the covariance matrix scales without turning, and the variance it reached, 100, is the scaling factor. That is not a coincidence, and section 29 proves it. The eigenvectors of C are the principal component directions, their eigenvalues are the variances along them, and sorting the eigenvalues sorts the components by explained variance.

The mean is (170, 68) and the deviations are (10, 10), (−10, −10), (5, −5) and (−5, 5), as the table in section 4 had them. Section 13’s builder is the same four people, and clicking a row there finds them in the plot. From the deviations, build the three sums.

Σ (height deviation)²                        = 100 + 100 + 25 + 25 = 250
Σ (weight deviation)²                        = 100 + 100 + 25 + 25 = 250
Σ (height deviation)(weight deviation)       = 100 + 100 − 25 − 25 = 150

Arranged as a matrix, those sums are the scatter matrix, which is the covariance totals before any divisor.

S = [ 250  150 ]
    [ 150  250 ]

Divide by the number of people, 4, for the population covariance, or by 3 for the sample covariance. The choice scales every eigenvalue by the same factor and changes neither the eigenvectors nor the shares.

C = [ 62.5  37.5 ]        C_sample = [ 83.33  50.00 ]
    [ 37.5  62.5 ]                   [ 50.00  83.33 ]

Check the directions on the scatter matrix, where the numbers stay whole, by feeding candidates through it.

S (1, 1)  = (250 + 150,  150 + 250) = (400, 400)  = 400 · (1, 1)
S (1, −1) = (250 − 150,  150 − 250) = (100, −100) = 100 · (1, −1)

Both come back unturned, so the components run along (1, 1) and (1, −1), normalised to u₁ = (0.707, 0.707) and u₂ = (0.707, −0.707), with scatter eigenvalues 400 and 100. On the population covariance the same directions come back with eigenvalues 100 and 25, which are actual variances, and 100 is the number the sweep in section 7 peaked at.

first sample variance = 100 × 4/3 ≈ 133.3
second sample variance = 25 × 4/3 ≈ 33.3

Each component’s share is its eigenvalue over the total, and the divisor cancels in the ratio.

share of component 1 = 400 / 500 = 100 / 125 = 0.8
share of component 2 = 100 / 500 =  25 / 125 = 0.2

Person 1 has deviation d = (10, 10). Their two scores are two dot products.

s₁ = d · u₁ = (10, 10) · (0.707, 0.707)  = 14.14
s₂ = d · u₂ = (10, 10) · (0.707, −0.707) = 0

They lie entirely along the first component, and their second coordinate is zero. Person 3, at (5, −5), scores 0 on the first and 7.07 on the second, the mirror case. The transform does this for everyone at once. With the centred rows stacked into a matrix X and the component directions as the columns of W, the scores are one product, and the widget writes it out for whichever person is selected.

Z = X W

…

Instead of keeping both scores, keep only s₁, and each person is one number. That is dimensionality reduction, and on real data with hundreds of correlated columns it is the difference between unworkable and workable. Tick the flatten toggle on the playground and the crowd moves onto the first component’s line.

One number per person is not the whole of what has to be stored, though. To do anything with those numbers later you also need the original feature means, the feature scales if standardisation was used, and the retained component direction. The scores are meaningless without the axis they are measured along.

Going back is two steps. Scale the retained direction by the score, then add the mean back.

d̂ = s₁ u₁
x̂ = x̄ + d̂

On the playground the green dots are exactly this, each person rebuilt from one number, with a dashed stub to where they really were. If standardisation was used, the scaling is reversed as well, in the opposite order it was applied.

The stub is the error, and its squared length is the natural measure of it.

errorᵢ = ‖xᵢ − x̂ᵢ‖²           total = Σᵢ ‖xᵢ − x̂ᵢ‖²

Now rotate the retained direction instead of accepting the first component, and watch the total. The widget in section 11 is the same instrument, and its discarded readout is this total. It is smallest at the first component, because for centred data the direction that keeps the most variance and the direction that leaves the smallest squared stubs are the same direction.

Pythagoras is why. Each person’s squared distance from the mean splits into a squared score plus a squared stub, and the distances do not depend on the direction, so making the scores as large as possible makes the stubs as small as possible.

‖dᵢ‖² = (dᵢ · u)² + stubᵢ²

With two features the choice is one component or two. The five-measurement dataset from section 1 makes it a real decision. Each component’s explained ratio is its eigenvalue over the total, and the cumulative ratio through m components is their running sum.

rⱼ = λⱼ / Σₖ λₖ            Rₘ = Σⱼ₌₁ᵐ rⱼ

…

The first component alone carries 78.2 percent of the standardised variance, two carry 87.2, three 94.4, four 99.0, and the reconstruction error falls from 65.3 to 38.5 to 16.8 to 2.9 to nothing. Keeping 95 percent is a common rule of thumb and not a law. The right count depends on what the reduced representation is for. A low-variance direction may carry the information a downstream task needs, and reconstruction quality and predictive performance answer different questions.

The crowd’s heights are in centimetres. Write them in metres or in millimetres and nothing about the people changes, but the numbers do, and the covariance matrix is built from the numbers.

…

In millimetres the height variance is 52,178 against weight’s 431, the first component is (−0.996, −0.089), almost pure height, and it claims 99.97 percent of the variance. In metres the height variance is 0.05, the first component is (0.011, 1.000), almost pure weight, and it claims all of it. Same people. Covariance PCA is sensitive to the numerical scale of its features, and a feature can dominate for no better reason than its units produce big numbers.

The usual repair is to standardise each feature first, dividing its deviations by its own standard deviation, so every feature starts with variance one. PCA on standardised features is PCA on the correlation matrix rather than the covariance matrix, and the question it answers changes from how the features vary together in their units to how they vary together relative to their own scales.

zᵢⱼ = (xᵢⱼ − x̄ⱼ) / sⱼ

…

Whether to standardise is a modelling choice, not a rule. Raw covariance PCA is the right instrument when the features share comparable units, when absolute variance means something, and when a feature that varies more should count for more. Standardised PCA is the right one when units differ, when relative variation is what matters, and when no feature should dominate merely because of its measurement scale.

PCA never sees a target. It finds directions of greatest variance, and those are not necessarily the directions that matter for a prediction. Here is a cloud built to make the point. Most of its spread runs along height, and the first component follows it, as it must. The two groups in it differ across the short axis, weight for a given height, which is almost exactly what a one-component reduction throws away.

…

The first component carries 98.5 percent of the variance and the second 1.5 percent. Keep one component and the fourteen people are placed on the first line, where the two groups interleave; the line that separates them cleanly is the one the reduction discarded as almost nothing.

The method runs on means, squared deviations and covariances, and every one of those gives a distant observation more than its share of the say. Add one stranger to the crowd, tall and light, and the mean moves, the covariance changes, the first component turns, and the shares follow.

…

One person in twelve moves the mean from (151.8, 53.4) to (155.4, 51.4), pulls the covariance of height and weight from 464 down to 349, turns the first component by 4.6 degrees, and cuts its share from 0.990 to 0.841. A squared deviation grows with the square of the distance, so a distant person contributes far more than a twelfth.

PCA finds linear structure. Twelve people along a curve have a one-dimensional pattern, their position along the curve, and no straight axis follows it. The first component does what it can, which is to lay a line across the curve and fold the two ends onto the same stretch of it.

…

The first component claims 68 percent of the variance here and the one-component reconstruction error is 3,316, against 108 on the crowd, which has one more person. The numbers say a straight axis is a poor summary; they do not say why, and the picture does.

The kernel PCA page is one answer to a curved pattern. And one more caution belongs here. A component can mix positive loadings, negative loadings and measurements with different meanings, and the mathematics does not care whether the mixture has a name. The loadings show the combination; whether it has a defensible interpretation is domain knowledge’s job, and a component can be useful for compression without ever earning a simple real-world label.

The first component through a cloud looks like a regression line and is not one. Regression chooses an outcome, predicts it from the other feature, and minimises the vertical misses in the picture. PCA treats the features symmetrically after preprocessing, finds a direction through the cloud, and minimises the perpendicular misses. Different objectives, different lines.

…

The third line makes the point twice. Predicting weight from height gives a slope of 0.890; predicting height from weight, drawn on the same axes, gives 0.928; the first component sits between them at 0.907. Reversing which feature is the outcome changes the regression line, because it changes which misses are being minimised, and PCA is not regression without a target. It is a different geometric problem with its own answer.

The claim to earn is that hunting variance leads to the covariance matrix’s eigenvectors, and it walks in four steps. First, the variance along a unit direction u. Each deviation d lands on u by the dot product, and expanding the square of that makes something familiar appear.

(d · u)² = (d₁u₁ + d₂u₂)² = d₁²u₁² + 2 d₁d₂ u₁u₂ + d₂²u₂²

Average over the people and the d sums become the statistics primer’s quantities, the two variances on the squared terms and the covariance on the cross term, which is the covariance matrix evaluated on u.

variance along u = uᵀ C u

Second, maximise that subject to u being a unit vector. Put the constraint in with a Lagrange multiplier and differentiate.

Third, which eigenvector. Multiply the equation by uᵀ on the left and the answer falls out.

uᵀ C u = λ uᵀ u = λ

The variance achieved along an eigenvector is its eigenvalue, so the maximum is the eigenvector with the largest eigenvalue, and that is the first principal component. Fourth, the second component solves the same problem with one more constraint, u₂ᵀu₁ = 0, which prevents it from rediscovering the same direction, and the argument continues for each later component inside the space perpendicular to all the earlier ones. A symmetric matrix’s eigenvectors for distinct eigenvalues are perpendicular to one another already, which is why the constraint costs nothing.

The route the derivation suggests is the one this page has walked. Centre the data matrix, optionally standardise its columns, form the covariance matrix, find its eigenvectors and eigenvalues, sort them by decreasing eigenvalue, keep the ones you want, and transform. It matches the concept step for step.

Practical implementations often skip the covariance matrix and decompose the centred data matrix directly.

X = U Σ Vᵀ

The columns of V are the principal directions, the squared singular values give the explained variances up to the divisor, and the scores are XV or, equivalently, UΣ. The singular value decomposition avoids forming XᵀX explicitly, which is better behaved numerically when the features are nearly collinear, and the result is the same PCA. The full derivation of that route belongs to another page.

A fitted PCA is a transformation, and a new person goes through the same one the training data went through. Use the training mean, the training standard deviations if scaling was fitted, and the retained training directions. Do not refit for each new observation, because then the coordinates would mean something different every time.

z_new = (x_new − x̄_training) W

The fit here keeps the means and, when it standardised, the standardiser, and its transform matches features by name, so a new observation with its columns in a different order is fine and one with a column missing is refused.

A complete implementation states the number of components, whether it centres, whether it scales, whether the covariance uses the population or sample divisor, whether it decomposes the covariance or the data matrix, how it orders components, what sign convention it applies if deterministic output is required, how it treats tied or nearly tied eigenvalues, its tolerances, and which of the scores, components, reconstructions and explained variances it returns.

Most of those choices change a number without changing the fit, which is what makes them worth stating. The divisor scales every variance by one factor and leaves the directions and the shares alone, as Part 7 measured, so two implementations can disagree about every variance they report and agree about every component. The sign is the same kind of difference. The fit on this page keeps whichever sign the solver returned, which is why the millimetre component of Part 11 came back as (−0.996, −0.089) with both loadings negative while the metre one came back as (0.011, 1.000), and Part 4 is the reason neither sign means anything. A reader comparing two implementations has to know which of their disagreements are about the data and which are about these conventions, and only a stated contract tells them.

Constant features are the case that deserves particular care, because the two preprocessing choices treat them oppositely. Raw, a constant column is harmless and simply carries nothing. Standardised, it is a division by zero, and the refusal is the right answer rather than a quiet substitution.