3 · PCA from scratch
Implement PCA via eigendecomposition and SVD in NumPy — demystify what sklearn does under the hood.
Implement PCA from scratch in numpy — eigendecomposition of the covariance matrix — to demystify what sklearn does in one line.
Without this:
Without seeing PCA from scratch, the 'top-k components' magic feels opaque and you can't debug shape issues or extend the algorithm.
In the previous lesson we used sklearn's PCA as a black box. Here we open it up. The recipe has five steps and maps directly to the math:
Step-by-step PCA from scratch:
- Standardise X → subtract column means, divide by column standard deviations.
- Covariance matrix →
C = np.cov(X_std.T)(d × d symmetric matrix). - Eigendecompose →
eigenvalues, eigenvectors = np.linalg.eigh(C)Useeigh(symmetric) noteig(generic) — eigh guarantees real eigenvalues. - Sort eigenvalues descending; reorder eigenvectors to match.
- Project →
Z = X_std @ eigenvectors[:, :k].
Reconstruction (going back to the original space):
X_reconstructed = Z @ eigenvectors[:, :k].T + X.mean(0)
Reconstruction error measures how much information k components retain:
relative_error = ||X - X_reconstructed||² / ||X||²
As k → d, the error drops to zero — perfect reconstruction. The interesting region is the "elbow" where increasing k stops buying much.
Python (in browser)
Manual PCA via eigendecomposition matches sklearn exactly (up to a sign flip per axis — eigenvectors are unique only up to sign).
Python runs entirely in your browser via Pyodide (~6 MB on first Run, cached after).
Python (in browser)
SVD-based PCA also matches sklearn — the rows of Vt are exactly the principal components (sklearn's components_ attribute).
Python runs entirely in your browser via Pyodide (~6 MB on first Run, cached after).
Python (in browser)
Reconstruction error vs k on the digits dataset — the curve elbows sharply near k=10, showing most pixel variance is concentrated in ~10 PCs.
Python runs entirely in your browser via Pyodide (~6 MB on first Run, cached after).
SVD vs eigendecomp tradeoffs — and how sklearn's svd_solver='auto' picks the right algorithm.
SVD — the full derivation of how singular vectors relate to eigenvectors of XᵀX and XXᵀ. Understanding this makes the 'rows of Vt = principal components' fact obvious.
After standardising X and computing the covariance matrix, the LARGEST eigenvalue represents...
- The five-step PCA recipe: standardise → covariance → eigh → sort descending → project. Every sklearn call reduces to this.
- Always use np.linalg.eigh (not eig) for symmetric matrices — it guarantees real eigenvalues and is faster.
- SVD of X gives the same principal components as eigendecomp of XᵀX/n — with better numerical stability.
- Reconstruction X_rec = Z @ Vt[:k] + mean lets you measure how much information k components retain via the relative reconstruction error.
- The reconstruction error curve shows an 'elbow' — the natural k for most datasets. Beyond the elbow, each additional component buys little.
Knowing PCA's internals lets you debug shapes (PCs are rows of components_?), understand why explained variance can suddenly look weird, and extend the algorithm — kernel PCA, sparse PCA, robust PCA are all variations on the same eigendecomposition theme.
If you remove it: PCA stays a black box you call but can't reason about — you can't debug shape mismatches, explain reconstruction quality, or adapt the algorithm for sparse or robust variants.