PCA from scratch

The last two lessons kept gesturing at one algorithm without ever naming it in full: eigenvectors sort a covariance matrix into directions of variance, SVD does the same for any matrix, and both, we kept saying, are “basically PCA.” Here is PCA itself, assembled from exactly those parts. You have fifty features and you want to see your data; PCA collapses high-dimensional points down to two or three dimensions while keeping as much of the variation as it can. The machinery is precisely the eigendecomposition of a covariance matrix from the eigenvalues lesson — but the recipe is so clean, and so constantly reached for, that it earns its own walkthrough.

The four-step recipe

1. Center the data:        X_c = X - X.mean(axis=0)
2. Compute covariance:     Σ   = (1/n) X_cᵀ X_c
3. Eigendecompose:         vals, vecs = eig(Σ)      (sort by val descending)
4. Project onto top-k:     Z   = X_c @ vecs[:, :k]

Four lines. Result: Z of shape (n, k) — every row is a low-dim version of the original sample.

From scratch

import numpy as np

rng = np.random.default_rng(0)

# Build 3D data that's mostly 2D -- a flat pancake in 3D space
n = 300
u = rng.normal(0, 3, n)
v = rng.normal(0, 1.5, n)
noise = rng.normal(0, 0.1, n)

# Three features, but the third is almost a function of the first two
X = np.column_stack([u + 0.3 * v, -0.5 * u + v, 0.2 * u - 0.4 * v + noise])

# --- PCA ---
X_c = X - X.mean(axis=0)
cov = (X_c.T @ X_c) / len(X_c)
vals, vecs = np.linalg.eigh(cov)             # eigh: symmetric -> faster, sorted ascending
order = np.argsort(-vals)
vals = vals[order]
vecs = vecs[:, order]

print("Eigenvalues:", vals.round(3))
explained = vals / vals.sum()
print("Explained variance ratio:", explained.round(3))
print("Cumulative:", np.cumsum(explained).round(3))

# Project onto top 2
k = 2
Z = X_c @ vecs[:, :k]
print("Projected shape:", Z.shape)

Eigenvalues: [1.2069e+01 2.5270e+00 7.0000e-03]
Explained variance ratio: [0.826 0.173 0.001]
Cumulative: [0.826 0.999 1.   ]
Projected shape: (300, 2)

The first two eigenvalues capture 99.9% of the variance between them; the third, 0.007, is almost nothing — the data only really lives in a 2-D plane tilted inside 3-D space. The projection Z is the view from straight “above” that plane.

Note: the principal components are new axes — linear combinations of all original features. They are not a subset of the original features.

Variance explained — how to pick k

vals[i] / vals.sum() tells you what fraction of total variance lives along the i-th principal direction. Pick k so the cumulative fraction crosses a threshold like 90% or 95%.

import numpy as np

rng = np.random.default_rng(42)
# Make 200 samples in 10 dimensions, but with only ~3 "real" directions
n, d_high, d_real = 200, 10, 3
A = rng.normal(size=(d_real, d_high))
Z_true = rng.normal(size=(n, d_real))
X = Z_true @ A + 0.1 * rng.normal(size=(n, d_high))

X_c = X - X.mean(axis=0)
cov = (X_c.T @ X_c) / len(X_c)
vals = np.sort(np.linalg.eigvalsh(cov))[::-1]

ratio = vals / vals.sum()
print("Per-component variance ratio:")
for i, r in enumerate(ratio):
    print(f"  PC{i+1:>2}: {r:.4f}   (cum {np.cumsum(ratio)[i]:.4f})")

k90 = (np.cumsum(ratio) < 0.90).sum() + 1
print(f"\nNeed k={k90} components to capture 90% of variance.")

Per-component variance ratio:
  PC 1: 0.5893   (cum 0.5893)
  PC 2: 0.2847   (cum 0.8740)
  PC 3: 0.1217   (cum 0.9957)
  PC 4: 0.0008   (cum 0.9965)
  PC 5: 0.0007   (cum 0.9972)
  PC 6: 0.0006   (cum 0.9978)
  PC 7: 0.0006   (cum 0.9985)
  PC 8: 0.0006   (cum 0.9990)
  PC 9: 0.0005   (cum 0.9995)
  PC10: 0.0005   (cum 1.0000)

Need k=3 components to capture 90% of variance.

The first three components hoover up 99.6% of the variance — exactly as designed, since the data was built from only three real directions. Everything past PC3 is the 0.1-scale noise we added.

Project a 10-feature dataset to 2D and plot it

This is what PCA is most often used for in practice: making a 2D scatter plot of high-dim data, colored by class, to see whether classes are separable.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)

# Three classes living in 10D -- different centers, shared noise
n_per = 100
centers = np.array([
    [3, 3, 0, 0, 0, 0, 0, 0, 0, 0],
    [-3, 3, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, -3, 0, 0, 0, 0, 0, 0, 0, 0],
], dtype=float)

X = np.concatenate([c + rng.normal(0, 1.0, (n_per, 10)) for c in centers])
y = np.repeat([0, 1, 2], n_per)

# PCA to 2D
X_c = X - X.mean(axis=0)
cov = (X_c.T @ X_c) / len(X_c)
vals, vecs = np.linalg.eigh(cov)
order = np.argsort(-vals)
Z = X_c @ vecs[:, order[:2]]
print("Projected to 2D:", Z.shape)

# Plot the 2D projection, coloured by class
fig, ax = plt.subplots(figsize=(5, 4))
for cls in [0, 1, 2]:
    ax.scatter(Z[y == cls, 0], Z[y == cls, 1], label=f"class {cls}", alpha=0.7, s=18)
ax.set_xlabel("PC 1")
ax.set_ylabel("PC 2")
ax.set_title("10D data projected to 2D via PCA")
ax.legend()
plt.tight_layout()
plt.show()

Projected to 2D: (300, 2)

Three blobs in 10-D space become three clearly separated blobs in the 2-D scatter (plt.show() draws it): PCA found the one plane where the classes pull apart the most. The reduction from 10 columns to 2 lost the noise dimensions and kept exactly the structure that matters.

The sklearn version

You’ll almost never write PCA by hand in production. Here’s the same thing in three lines:

# sklearn.decomposition.PCA gives you all of this in three lines. To see
# exactly what it does under the hood, here is the same API in a little NumPy:

import numpy as np

rng = np.random.default_rng(0)
X = rng.normal(size=(100, 8))

class MiniPCA:
    def __init__(self, n_components):
        self.k = n_components
    def fit(self, X):
        Xc = X - X.mean(axis=0)
        self.mean_ = X.mean(axis=0)
        vals, vecs = np.linalg.eigh((Xc.T @ Xc) / len(Xc))
        order = np.argsort(-vals)
        self.components_ = vecs[:, order[:self.k]].T   # shape (k, d)
        self.explained_variance_ratio_ = (vals[order] / vals.sum())[:self.k]
        return self
    def transform(self, X):
        return (X - self.mean_) @ self.components_.T

pca = MiniPCA(n_components=3).fit(X)
print("components_ shape:", pca.components_.shape)
print("variance ratio:", pca.explained_variance_ratio_.round(4))
print("transformed shape:", pca.transform(X).shape)

components_ shape: (3, 8)
variance ratio: [0.1945 0.1676 0.1493]
transformed shape: (100, 3)

(Here X is pure random noise in 8 dimensions, so no direction truly dominates — the top three components each explain a roughly equal ~18%.) sklearn.decomposition.PCA follows exactly this API; under the hood it uses SVD rather than eigendecomposition of XᵀX (more numerically stable, per the last lesson), but conceptually it is the same four-step algorithm.

When PCA falls down

PCA is a linear method. It finds the best linear projection. If your data lives on a curved manifold — a Swiss roll, faces in pose-space, a spiral — the best 2D linear projection might smush points from different parts of the manifold on top of each other.

For non-linear structure, you reach for:

• t-SNE — preserves local neighborhoods, great for visualization.
• UMAP — similar idea, often faster and preserves more global structure.
• Autoencoders — let a neural net learn the nonlinear embedding.

PCA is still the right first move — fast, deterministic, gives you a variance-explained number. Run it first. If the projection looks bad, then reach for t-SNE/UMAP.

In one breath

PCA is four steps: center the data, build the covariance Σ = (1/n)XᵀX, eigendecompose it (use eigh — covariance is symmetric) and sort by eigenvalue, then project onto the top-k eigenvectors. The eigenvalues are variances per direction, so vals / vals.sum() is the variance-explained ratio — choose k where the cumulative sum crosses 90–95%. The principal components are new axes (combinations of all features, not a subset), and the result lets you plot fifty-dimensional data on a flat 2-D scatter. PCA is linear, so it shines on flat structure and fails on curved manifolds (a Swiss roll), where you reach instead for t-SNE, UMAP, or autoencoders.

Practice

A question to carry forward

That closes the linear-algebra half of the toolkit. Step back and notice what every lesson in it shared: vectors, matrices, projections, eigenvectors, SVD, PCA — each was a way of describing structure that already sits still. Given a fixed pile of data, we found its directions, its rank, its best low-dimensional shadow. Nothing in this chapter ever moved.

But a model does not sit still; it learns. It starts with bad parameters and improves them — nudging each weight a little, again and again, until the error stops shrinking. And that raises a question linear algebra simply cannot answer: standing at one setting of a knob with a certain error, which way should you turn the knob, and how far? Here is the thread into the next chapter: what is the derivative — the slope that whispers “tilt this way to go downhill” — and how does following it, one small step at a time, become gradient descent, the single algorithm that trains very nearly everything in machine learning?