Principal component analysis

What this chapter covers

k-means asked where the groups are. This chapter asks a different unsupervised question: which directions in the data actually carry information, and which ones are just noise you can throw away. That's dimensionality reduction, and principal component analysis is the workhorse everyone learns first.

The pitch fits in one line. Find the axes along which the data varies most, keep the top few, and drop the rest. A four-column dataset becomes two columns you can plot, and if you chose right you lose almost nothing worth keeping. We build it by hand in NumPy — center the data, form the covariance matrix, take its eigendecomposition, sort by eigenvalue, project onto the top directions — then hand the same job to scikit-learn and check that the two agree down to the sign. And because the whole thing is a rotation in the plane, we get to watch it: a tilted cloud of points, the two principal axes snapping into place, the cloud spinning until the long axis lies flat, and every point dropping onto that one line.

The concept runs on a toy 2-D cloud so you can see each point move. The real result runs on iris — 150 flowers, four measurements each — squashed from four dimensions down to two and colored by species, to show that the structure survives the squeeze.

A bit of history

PCA has been invented at least twice. Karl Pearson got there first, in 1901, in a paper with the wonderfully literal title "On lines and planes of closest fit to systems of points in space." His framing was geometric: given a cloud of points, find the line, then the plane, that sits closest to all of them in the least-squares sense. That best-fit line is the first principal component. He was doing it with pencil and biometric data, no computer in sight, and the whole argument is about perpendicular distances to a line.

Harold Hotelling came at it from the other side in 1933, in the Journal of Educational Psychology, and this is where the modern name and machinery come from. Hotelling wasn't fitting a line; he was looking at a pile of correlated variables — test scores, mostly — and asking for a smaller set of uncorrelated variables that captured the same information. He derived them as the directions of maximum variance, showed they were the eigenvectors of the covariance matrix, and called them the principal components. Two descriptions, one algorithm: the line closest to the points is the same as the direction the points spread out along most. That equivalence is the quiet heart of the method, and it's why the same tool gets pitched sometimes as compression and sometimes as decorrelation.

The intuition

Look at a cloud of points that's stretched and tilted, like a cigar lying at an angle. Two numbers describe each point, its x and its y, but they're redundant: tell me where a point sits along the long direction of the cigar and I can guess the rest, because the cloud is thin the other way. The two original axes are a bad description of this data. They're not wrong, they're just not aligned with anything the data cares about.

PCA finds a better pair of axes. The first one, PC1, points along the direction the cloud is most spread out — the length of the cigar. The second, PC2, is forced to be perpendicular to it and picks up whatever variance is left, the thin width. These aren't two of your original features; they're new directions, mixtures of the old ones, chosen so that the first captures as much spread as possible and each next one mops up the leftovers. Once you're on those axes, the data tells you how much it cares about each: nearly all the variance sits on PC1, almost none on PC2. So you keep PC1, drop PC2, and you've gone from two numbers per point to one while barely disturbing the picture.

Here's the cloud we'll work on — 150 points, a correlated Gaussian blob tilted off the axes on purpose:

Your eye already found the long direction. PCA is how a computer finds the same line, and it does it with a covariance matrix and an eigendecomposition.

The math

Start with $N$ points $x_1, \dots, x_N$, each a vector in $D$ dimensions. PCA measures variance, and variance is measured around the mean, so the first move is to center. Let $\bar{x} = \frac{1}{N}\sum_{i=1}^{N} x_i$ be the mean and subtract it from every point:

Now the cloud is centered at the origin. Stack the centered points as the rows of a matrix $\tilde{X}$ (shape $N \times D$) and form the covariance matrix:

$C$ is $D \times D$ and symmetric. Its diagonal entry $C_{aa}$ is the variance of feature $a$; the off-diagonal $C_{ab}$ is the covariance between features $a$ and $b$, how much they move together. Everything PCA knows about the data's shape is in this matrix. The principal directions are its eigenvectors — the vectors $v$ that $C$ merely stretches without rotating:

Each eigenvector $v_j$ is a principal direction, and its eigenvalue $\lambda_j$ is exactly the variance of the data along that direction. Because $C$ is symmetric the eigenvectors are orthogonal, so the $v_j$ form a clean new set of perpendicular axes. Sort them so $\lambda_1 \ge \lambda_2 \ge \dots$ and $v_1$ is the direction of greatest variance, $v_2$ the next, and so on.

To reduce to $k$ dimensions, keep the top $k$ eigenvectors as the columns of a matrix $V_k$ and project each centered point onto them:

The vector $z_i$ holds the point's coordinates on the new axes — its scores. That's the reduction: $D$ numbers in, $k$ numbers out. To judge what you kept, the explained-variance ratio of component $j$ is its eigenvalue over the total:

These sum to 1 across all $D$ components. The cumulative sum over the top $k$ is the fraction of the data's total variance your reduced version still carries — the single number you use to decide how many components are enough.

What it's good at, what it isn't

The reason PCA is everywhere is that it does several unrelated-looking jobs with one computation. It compresses: swap $D$ correlated features for $k$ that carry almost all the variance. It visualizes: force high-dimensional data down to two or three axes you can actually plot. It decorrelates: the new axes are orthogonal by construction, so the components are uncorrelated even when the original features were tangled together, which is a gift for any model that chokes on collinearity. And it denoises: the small-eigenvalue directions are usually where the noise lives, so dropping them can leave a cleaner signal than you started with. One eigendecomposition, four uses.

The catch is one word: linear. PCA only knows how to rotate and project. Its components are straight-line combinations of your features, and the structure it finds is whatever a covariance matrix can see — variance and pairwise correlation, nothing else. Hand it data that curls, a spiral or a shape that lives on a bent sheet, and PCA flattens it wrong, because no rotation of the axes untangles a curve. It also equates variance with importance, which isn't always true; a low-variance direction can carry the signal you actually care about while a high-variance one is just the loudest, most useless feature drowning everything else out. And the components are combinations of every original feature, so they're often hard to name — PC1 is some weighted blend of sepal length, petal length, and the rest, which summarizes well and explains poorly.

The data

Two datasets, for two jobs. The concept animation runs on the tilted 2-D cloud above: 150 points from a Gaussian stretched along one axis and rotated about 33 degrees off horizontal, seeded so the snapshot is stable. Two dimensions so the whole rotation fits on the page and you can follow individual points.

The real result runs on iris — Fisher's 1936 measurements of 150 flowers, 50 each from three species, four features apiece: sepal length, sepal width, petal length, petal width. Four dimensions, which you can't draw. PCA's whole promise is to fix that: collapse the four measurements to two principal components and plot every flower on a plane. The species labels come along for the ride but PCA never sees them — it only gets the four numbers per flower. We use the labels at the very end, to color the plot and check whether the reduction kept the three species apart.

Build it, one function at a time

Seven short functions, and they line up exactly with the math. Center, covary, eigendecompose, then the three things you do with the result: measure variance, project, reconstruct. The first step is the mean subtraction, and we hand back the mean because we'll need it later to undo the shift:

def center(X):
    """Subtract the column means so the cloud sits at the origin.

    PCA is about variance — directions of spread — and spread is measured
    around the mean. Returns the centered data and the mean vector (kept so we
    can undo the shift when we reconstruct).
    """
    mean = X.mean(axis=0)
    return X - mean, mean

With the cloud centered on the origin, the covariance matrix is one matrix product. Divide by $N-1$, the unbiased estimator, so we match what numpy and sklearn report:

def covariance(Xc):
    """The (D, D) covariance matrix of already-centered data.

    Entry (a, b) is the covariance between feature a and feature b. The
    diagonal is each feature's variance; off-diagonals say how features move
    together. We divide by N-1 (the unbiased estimator), which is also what
    numpy and sklearn use.
    """
    n = len(Xc)
    return (Xc.T @ Xc) / (n - 1)

That $D \times D$ matrix holds the entire shape of the data. Now the one line that turns it into directions. The covariance matrix is symmetric, so eigh is the right routine — it returns real eigenvalues and orthonormal eigenvectors, and we flip its ascending order to descending so component 0 is the direction of greatest variance:

def eigendecompose(C):
    """Eigenvalues and eigenvectors of the covariance matrix, sorted.

    C is symmetric, so `eigh` is the right tool — real eigenvalues, orthonormal
    eigenvectors. `eigh` returns them in ascending order; we flip to descending
    so component 0 is the direction of greatest variance. Eigenvectors are the
    COLUMNS of the returned matrix, one principal direction each.
    """
    eigvals, eigvecs = np.linalg.eigh(C)
    order = np.argsort(eigvals)[::-1]
    return eigvals[order], eigvecs[:, order]

The eigenvectors come back as columns, one principal direction each, already sorted. The eigenvalues are the variances along them, and turning those into the fraction each component explains is just division:

def explained_variance_ratio(eigvals):
    """Fraction of total variance each component accounts for.

    Each eigenvalue IS the variance along its eigenvector, so the ratio is just
    the eigenvalue over their sum. These add to 1 across all components; the
    cumulative sum tells you how much you keep by stopping at k.
    """
    return eigvals / eigvals.sum()

That's the number that decides how many components to keep. Now the reduction itself — project the centered points onto the top-k directions to get the scores:

def project(Xc, components, k):
    """Project centered data onto the top-k principal directions.

    `components` holds eigenvectors as columns (from eigendecompose). Taking the
    first k columns and multiplying gives the scores: each point's coordinates
    in the new k-dimensional basis. This is the dimensionality reduction —
    (N, D) in, (N, k) out.
    """
    return Xc @ components[:, :k]

(N, D) in, (N, k) out. That's the whole point of the method in one matrix multiply. Going the other way is reconstruction: take the low-dimensional scores back up into the original space and add the mean back. With k < D it can't be exact — what returns is the closest rank-k approximation, and the gap is the variance you dropped:

def reconstruct(scores, components, mean, k):
    """Map the k-dim scores back into the original feature space.

    Undo the projection (multiply by the components' transpose) and add the
    mean back. With k < D this is lossy: what comes back is the closest
    rank-k approximation of the original point. The gap between the two is the
    reconstruction error, and it equals the variance in the components we
    dropped.
    """
    return scores @ components[:, :k].T + mean

Last, the wrapper that runs the pipeline end to end and returns everything the chapter needs — mean, components, eigenvalues, explained-variance ratio, and the scores — in one dict:

def pca(X, k):
    """Full PCA: center, covariance, eigendecompose, project onto top-k.

    Returns everything the chapter needs in one dict — the mean, the sorted
    principal directions (columns), the eigenvalues (per-component variance),
    the explained-variance ratio, and the k-dim scores. This is the whole
    algorithm; every field comes from the small functions above.
    """
    Xc, mean = center(X)
    C = covariance(Xc)
    eigvals, components = eigendecompose(C)
    scores = project(Xc, components, k)
    return {
        "mean": mean,
        "components": components,          # (D, D), eigenvectors as columns
        "eigenvalues": eigvals,            # variance along each component
        "explained_variance_ratio": explained_variance_ratio(eigvals),
        "scores": scores,                  # (N, k) reduced coordinates
    }

No loops, no iteration, no convergence to wait on. Unlike k-means, PCA is a closed-form computation: one eigendecomposition and you're done. Which makes it odd to animate — there are no steps to replay. So instead we animate the meaning.

Watch it work

PCA doesn't iterate, so there's nothing to step through the way k-means marches its centroids. What there is to see is the geometry: what "project onto the principal components" actually does to the cloud. The animation below runs the real pca function on the tilted 2-D data and plays back the transform in three acts. Every coordinate is computed, not staged — the axes are the true eigenvectors, the rotation is the real change of basis, the residuals are the actual discarded component.

First the axes appear. The long orange arrow is PC1, laid down the direction the cloud is most spread out; it alone accounts for 92.7% of the variance. The shorter cyan arrow is PC2, perpendicular to it, holding the remaining 7.3%. Then the whole cloud rotates, arrows and all, until PC1 lies flat along the horizontal — that's the change of basis, the moment the data's natural axes become the plot's axes. Finally every point drops straight down onto the PC1 line and the gray stub left behind is its residual, the PC2 coordinate we're choosing to throw away. Points are shaded by their position along PC1 so you can track where each one lands.

Watch the last act. Almost every point barely moves when it drops — the residual stubs are tiny — because there was almost nothing in the PC2 direction to begin with. That's the 92.7% made visible: after the projection each point is a single number, its spot on the line, and the picture is nearly unchanged. Reset and run it again; it's deterministic, the same rotation every time, because there's no seed and no random start here, just linear algebra.

That's PCA in two dimensions, where dropping to one is a party trick. The value shows up when the cloud lives in four dimensions, or four hundred, and the same move buys you a plot you couldn't otherwise draw.

The full implementation

The whole file, no library, top to bottom. This is exactly what the animation ran and what the iris numbers below come from:

"""Principal component analysis, built from scratch.

Unsupervised again — no labels, just coordinates. PCA finds the axes along
which the data varies most, so we can rotate onto those axes and keep only the
first few. Pure NumPy: center the data, form the covariance matrix, take its
eigendecomposition, sort by eigenvalue, project onto the top components, and
measure how much variance we kept.

Every function below appears in the chapter one step at a time (the
`# region:` markers are what the book's include directives pull in).
"""

import numpy as np
import pandas as pd


# region: center
def center(X):
    """Subtract the column means so the cloud sits at the origin.

    PCA is about variance — directions of spread — and spread is measured
    around the mean. Returns the centered data and the mean vector (kept so we
    can undo the shift when we reconstruct).
    """
    mean = X.mean(axis=0)
    return X - mean, mean
# endregion


# region: covariance
def covariance(Xc):
    """The (D, D) covariance matrix of already-centered data.

    Entry (a, b) is the covariance between feature a and feature b. The
    diagonal is each feature's variance; off-diagonals say how features move
    together. We divide by N-1 (the unbiased estimator), which is also what
    numpy and sklearn use.
    """
    n = len(Xc)
    return (Xc.T @ Xc) / (n - 1)
# endregion


# region: eig
def eigendecompose(C):
    """Eigenvalues and eigenvectors of the covariance matrix, sorted.

    C is symmetric, so `eigh` is the right tool — real eigenvalues, orthonormal
    eigenvectors. `eigh` returns them in ascending order; we flip to descending
    so component 0 is the direction of greatest variance. Eigenvectors are the
    COLUMNS of the returned matrix, one principal direction each.
    """
    eigvals, eigvecs = np.linalg.eigh(C)
    order = np.argsort(eigvals)[::-1]
    return eigvals[order], eigvecs[:, order]
# endregion


# region: explained_variance
def explained_variance_ratio(eigvals):
    """Fraction of total variance each component accounts for.

    Each eigenvalue IS the variance along its eigenvector, so the ratio is just
    the eigenvalue over their sum. These add to 1 across all components; the
    cumulative sum tells you how much you keep by stopping at k.
    """
    return eigvals / eigvals.sum()
# endregion


# region: project
def project(Xc, components, k):
    """Project centered data onto the top-k principal directions.

    `components` holds eigenvectors as columns (from eigendecompose). Taking the
    first k columns and multiplying gives the scores: each point's coordinates
    in the new k-dimensional basis. This is the dimensionality reduction —
    (N, D) in, (N, k) out.
    """
    return Xc @ components[:, :k]
# endregion


# region: reconstruct
def reconstruct(scores, components, mean, k):
    """Map the k-dim scores back into the original feature space.

    Undo the projection (multiply by the components' transpose) and add the
    mean back. With k < D this is lossy: what comes back is the closest
    rank-k approximation of the original point. The gap between the two is the
    reconstruction error, and it equals the variance in the components we
    dropped.
    """
    return scores @ components[:, :k].T + mean
# endregion


# region: pca
def pca(X, k):
    """Full PCA: center, covariance, eigendecompose, project onto top-k.

    Returns everything the chapter needs in one dict — the mean, the sorted
    principal directions (columns), the eigenvalues (per-component variance),
    the explained-variance ratio, and the k-dim scores. This is the whole
    algorithm; every field comes from the small functions above.
    """
    Xc, mean = center(X)
    C = covariance(Xc)
    eigvals, components = eigendecompose(C)
    scores = project(Xc, components, k)
    return {
        "mean": mean,
        "components": components,          # (D, D), eigenvectors as columns
        "eigenvalues": eigvals,            # variance along each component
        "explained_variance_ratio": explained_variance_ratio(eigvals),
        "scores": scores,                  # (N, k) reduced coordinates
    }
# endregion


def load_cloud(path="../data/cloud.csv"):
    """The 2-D toy cloud: a tilted, correlated Gaussian blob (x, y)."""
    return pd.read_csv(path)


def load_iris_csv(path="../data/iris.csv"):
    """The real dataset: 150 iris flowers, 4 measurements + species label.

    The species column is only used to color the 2-D projection at the end —
    PCA itself never sees it. This is unsupervised.
    """
    return pd.read_csv(path)

The library version

Nobody hand-rolls PCA in production, and once you've built it you don't need to. scikit-learn's PCA computes the same thing a smarter way — it takes the SVD of the centered data directly instead of forming the covariance matrix and eigendecomposing it, which is more numerically stable and lands on the identical answer:

def sklearn_pca(X, k):
    """Fit PCA and return the pieces our scratch version returns.

    `components_` are the principal directions as ROWS (note: sklearn stores
    them transposed relative to our column layout). `explained_variance_ratio_`
    is the per-component variance fraction, and `transform` gives the scores.
    """
    model = PCA(n_components=k)
    scores = model.fit_transform(X)
    return {
        "components": model.components_,                    # (k, D), as rows
        "explained_variance_ratio": model.explained_variance_ratio_,
        "scores": scores,                                   # (N, k)
    }

Two bookkeeping differences to know before comparing. sklearn stores the components as rows, where our version keeps them as columns, so one of them is a transpose of the other. And the important one: an eigenvector is only defined up to sign. Flip $v$ to $-v$ and it's still a unit vector pointing along the same axis with the same eigenvalue — the line doesn't care which way the arrow points. So sklearn's PC1 can come out pointing the opposite direction from ours, which flips the sign of every score on that component, while the axis, the variance, and the explained-variance ratio are all identical. When you compare two PCA implementations and the components look negated, that's not a bug, it's the sign ambiguity, and you align for it before checking.

Scratch versus library

Same iris data, both reduced to two components — our covariance eigendecomposition against sklearn's SVD. The bars are the explained-variance ratio of each component:

They're the same to the digit: PC1 at 0.9246, PC2 at 0.0531, from both. That isn't rounding hiding a gap. After aligning the sign, the largest disagreement between our components and sklearn's is on the order of 1e-14, the explained variance matches to machine precision, and the projected coordinates line up to better than 1e-4 — the two implementations are computing the same thing, one via the covariance eigendecomposition and one via SVD. The one visible difference was the sign: our PC1 came out flipped relative to sklearn's, exactly the ambiguity from the last section, which we correct before the numbers agree.

Now the payoff. Here's iris projected onto those two components, all four measurements collapsed to a plane, each flower colored by its true species:

PCA never saw the colors. It got four numbers per flower and returned two, and yet the three species fall into place — setosa clean off on its own, versicolor and virginica adjacent but mostly separable. The first two components carry 97.77% of the total variance (92.46% on PC1, 5.31% on PC2), so those two axes are almost the whole dataset; the two we dropped hold barely 2%. That's the promise kept: four dimensions down to two, a picture you can actually look at, and the structure that mattered survived the reduction. The average point reconstructs from its two scores with a squared error of about 0.10 across four features — that leftover 2% of variance, and nothing more.

The scree plot is the tool for choosing $k$ without guessing. Plot each component's variance and the running cumulative, and look for where the curve flattens:

The bars crater after PC1 and the cumulative line is already at 0.98 by PC2, then crawls to 1.0 across PC3 and PC4. That flattening is the signal: two components are plenty for iris, and the third and fourth are almost pure diminishing returns. In the wild this curve is your budget — pick the $k$ where the cumulative crosses whatever variance you're willing to keep, 90%, 95%, 99%, and stop.

Takeaways

PCA is the first tool I reach for when a dataset is wide and I need to see it, or when a model is drowning in correlated features. It's cheap, it's a single eigendecomposition with no hyperparameters to babysit and no seed to get unlucky on, and it does four jobs at once: it compresses, it draws high-dimensional data on a plane, it hands a downstream model uncorrelated inputs, and it can shave off noise by dropping the small-variance directions. For visualization and for decorrelation before a linear model, it's close to a free win, and the scree plot tells you how many components to keep instead of leaving it to taste.

The word to keep in mind is linear. PCA rotates and projects, and that's all it does. It reads the world through a covariance matrix, so it sees variance and pairwise correlation and nothing else — no curves, no manifolds, no notion that your data might live on a bent sheet rather than a flat one. When the structure is nonlinear, PCA flattens it and reports high explained variance while quietly destroying the thing you cared about. It also treats variance as importance, which is a useful default and an occasional trap: the loudest direction isn't always the informative one. That's exactly the boundary the next methods push past. t-SNE and UMAP keep the goal — few dimensions you can plot — but bend the map to preserve local neighborhoods instead of global variance, so a spiral stays a spiral. Kernel PCA does the covariance trick in a nonlinear feature space. Autoencoders learn the whole reduce-then-reconstruct pipeline as a neural net, with none of the straight-line constraint. Every one of them is a loosened PCA, which is the best possible reason to have built the linear one by hand first.