BrevFeed
Sign in

48. Principal Component Analysis (PCA) via SVD/eigen-decomposition

📖 Reading · 10 min
💡 Every code box below is live — edit it and hit Run.

§21.4 called SVD's low-rank approximation "the mathematical foundation of PCA" without elaborating. This lesson makes good on that promise: Principal Component Analysis finds the directions along which a dataset varies most — and it turns out to be nothing more than an eigendecomposition of a covariance matrix, or equivalently, an SVD of centered data.

Setup: centering and covariance

Given n data points in \mathbb{R}^d, stacked as rows of a matrix X (n\times d), center the data by subtracting the mean of each column: \tilde X=X-\vec1\bar{\vec x}^T (every row shifted so the new column means are all 0). The covariance matrix is

C=\frac1{n-1}\tilde X^T\tilde X

C is symmetric (C^T=C, from (\tilde X^T\tilde X)^T=\tilde X^T\tilde X, §16.2) and PSD (§21.2's B^TB fact, with B=\tilde X/\sqrt{n-1}) — so by §21.0's Spectral Theorem, C has an orthonormal eigenbasis with nonnegative eigenvalues.

Principal components are eigenvectors of the covariance matrix

The principal components are C's eigenvectors \vec v_1,\dots, \vec v_d, ordered by decreasing eigenvalue \lambda_1\ge\cdots\ge \lambda_d\ge0. \vec v_1 is the direction of maximum variance: the variance of the data projected onto a unit direction \vec u (§20.4's projection) is exactly \vec u^TC\vec u — a quadratic form (§21.3) — and §21.3's diagonalization argument shows this is maximized exactly at \vec u=\vec v_1, with maximum value \lambda_1. Each subsequent \vec v_i captures the most remaining variance, orthogonal to every earlier component (guaranteed by §21.0's automatic orthogonality of eigenvectors from distinct eigenvalues).

Dimensionality reduction: projecting the data onto the top k principal components (§20.4's projection formula, stopped at k terms) gives the best possible k-dimensional summary of the data, in the sense of preserving the most variance — this is precisely §21.4's Eckart-Young "best low-rank approximation" theorem, applied to \tilde X specifically.

The SVD route

Equivalently, and more numerically stable (§22.2's condition-number argument: computing C=\tilde X^T\tilde X squares the condition number, exactly as in §20.5/§20.6's least-squares story), take the SVD of the centered data directly: \tilde X=U\Sigma V^T. Then

C=\frac1{n-1}\tilde X^T\tilde X=\frac1{n-1}V\Sigma U^TU\Sigma V^T=\frac1{n-1}V\Sigma^2V^T

using U^TU=I (§20.4) — so V's columns are exactly the principal components, and \lambda_i=\sigma_i^2/(n-1). This is §21.4's "SVD generalizes the Spectral Theorem" fact, applied concretely: PCA never needs to form C at all, sidestepping the squared condition number by computing the SVD of \tilde X directly.

Doing it in Python

import numpy as np

rng = np.random.default_rng(0)
n = 200
# Data correlated along one dominant direction, with noise
t = rng.normal(0, 3, n)
X = np.column_stack([2*t + rng.normal(0, 1, n), t + rng.normal(0, 1, n)])

X_centered = X - X.mean(axis=0)
C = (X_centered.T @ X_centered) / (n - 1)

eigvals, eigvecs = np.linalg.eigh(C)
order = np.argsort(eigvals)[::-1]   # descending
eigvals, eigvecs = eigvals[order], eigvecs[:, order]

print("eigenvalues (variance along each component):", [round(v, 4) for v in eigvals.tolist()])
print("PC1 direction:", [round(v, 4) for v in eigvecs[:, 0].tolist()])
print("fraction of variance explained by PC1:", round(float(eigvals[0] / eigvals.sum()), 4))

eigenvalues (variance along each component): [42.4627, 0.9604]
PC1 direction: [-0.8883, -0.4592]
fraction of variance explained by PC1: 0.9779

Confirming the SVD route gives the identical components:

import numpy as np

rng = np.random.default_rng(0)
n = 200
t = rng.normal(0, 3, n)
X = np.column_stack([2*t + rng.normal(0, 1, n), t + rng.normal(0, 1, n)])
X_centered = X - X.mean(axis=0)

U, s, Vt = np.linalg.svd(X_centered, full_matrices=False)
variance_from_svd = (s**2) / (n - 1)

print("variance from SVD singular values:", [round(v, 4) for v in variance_from_svd.tolist()])
print("PC1 direction (V's first row):", [round(v, 4) for v in Vt[0, :].tolist()])

variance from SVD singular values: [42.4627, 0.9604]
PC1 direction (V's first row): [-0.8883, -0.4592]

Worked example

For covariance matrix C=\begin{pmatrix}4&2\\2&4\end{pmatrix}, find the first principal component and the fraction of variance it explains.

\operatorname{tr}C=8, \det C=16-4=12. \lambda^2-8\lambda+12=0\Rightarrow(\lambda-2)(\lambda-6)=0\Rightarrow\lambda=6,2.

\lambda_1=6: (C-6I)\vec v=\vec0: \begin{pmatrix}-2&2\\2&-2\end{pmatrix}\vec v=\vec0\Rightarrow v_1=v_2. Normalized: \vec v_1=\tfrac1{\sqrt2}(1,1).

\boxed{\text{PC1}=\tfrac1{\sqrt2}(1,1),\quad\text{variance explained}=\frac6{6+2}=75\%}

Sanity check. \operatorname{tr}C=6+2=8 ✓ matches directly. PC1's direction, (1,1)/\sqrt2, makes sense given C's structure: the two variables have equal individual variance (4 each, the diagonal) and positive covariance (2, off-diagonal), so they tend to move together — the direction of maximum combined variance is exactly the "both increase together" diagonal direction, matching intuition without needing the eigenvector calculation to confirm it.

Your turn

1. If C=\begin{pmatrix}5&0\\0&2\end{pmatrix} (already diagonal), what is the first principal component, without any computation?

2. Why must PCA's eigenvalues always be \ge0?

3. True or false: PCA's principal components are always orthogonal to each other.

Solutions

1. \vec v_1=(1,0)=\vec e_1, with variance 5 (the larger diagonal entry) — for a diagonal covariance matrix, the principal components are exactly the coordinate axes themselves, ordered by variance, since a diagonal matrix's eigenvalues are its diagonal entries (§19.1) and its eigenvectors are the standard basis vectors directly.

2. C is PSD (shown via the B^TB argument above), and §21.2's eigenvalue test says PSD is exactly "every eigenvalue \ge0" — negative variance is nonsensical, and the algebra enforces that directly rather than needing a separate check.

3. True, whenever eigenvalues are distinct (§21.0's automatic orthogonality for symmetric matrices with distinct eigenvalues) — and even when eigenvalues repeat, since §21.0's full Spectral Theorem guarantees an orthonormal eigenbasis can always be chosen (Gram-Schmidt inside a repeated eigenspace, if needed). So yes, unconditionally: PCA components are always orthogonal, by construction.

Check yourself in code

For covariance matrix C=\begin{pmatrix}3&1\\1&3\end{pmatrix}, find the eigenvalues and the fraction of variance explained by the first principal component.

Print exactly this:

eigenvalues: [4.0, 2.0]
variance explained by PC1: 0.6667

import numpy as np

C = np.array([[3., 1.], [1., 3.]])
eigvals = np.linalg.eigvalsh(C)[::-1]   # descending order
print("eigenvalues:", eigvals.tolist())
# print the fraction of total variance explained by the largest eigenvalue

import numpy as np

C = np.array([[3., 1.], [1., 3.]])
eigvals = np.linalg.eigvalsh(C)[::-1]
print("eigenvalues:", eigvals.tolist())

fraction = eigvals[0] / eigvals.sum()
print("variance explained by PC1:", round(float(fraction), 4))

PCA's principal components are the eigenvectors of the (PSD, by construction) covariance matrix, ordered by eigenvalue — each one the direction of maximum remaining variance, always orthogonal to the others by §21.0's Spectral Theorem. Computing them via SVD of the centered data directly, rather than forming the covariance matrix, avoids §22.2's condition-number squaring, exactly as in §20.6's least-squares story.

Next: applying eigenvalues to a discrete structure instead of continuous data — spectral graph theory, where a graph's connectivity is read off the eigenvalues of a matrix built from its edges.

← Previous Next →