Principal Component Analysis (PCA)

pca is the linear algebra exam question that escaped into industry. given a cloud of points in \(\mathbb{R}^d\), it finds the orthogonal directions along which the cloud spreads the most, and lets you throw away the rest. 𐃏 two apparently different questions — “which directions carry the most variance?” and “which subspace loses the least when i project onto it?” — turn out to have the same answer, and that answer is an eigendecomposition.

setup#

data \(x_1, \dots, x_n \in \mathbb{R}^d\), assumed centred: \(\sum_i x_i = 0\) (subtract the mean first; everything below silently requires this). the sample covariance matrix is

\begin{equation} S = \frac{1}{n}\sum_{i=1}^{n} x_i x_i^\top \in \mathbb{R}^{d \times d}, \end{equation}

symmetric and positive semi-definite, so it has an orthonormal eigenbasis \(S v_j = \lambda_j v_j\) with \(\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_d \ge 0\). 𐃏 the variance of the data projected onto a unit vector \(w\) is \(\operatorname{Var}(w^\top x) = w^\top S w\).

variance maximisation#

the first component#

find the unit direction of maximal projected variance:

\begin{equation} \max_{w \in \mathbb{R}^d} \; w^\top S w \quad \text{subject to} \quad w^\top w = 1. \end{equation}

the constraint matters — without it, scale \(w\) up and the “variance” grows forever. form the lagrangian

\begin{equation} \mathcal{L}(w, \lambda) = w^\top S w - \lambda\,(w^\top w - 1), \end{equation}

and set the gradient to zero:

\begin{equation} \nabla_w \mathcal{L} = 2 S w - 2 \lambda w = 0 \;\;\Longleftrightarrow\;\; S w = \lambda w. \end{equation}

every stationary point is an eigenvector of \(S\), and at such a point the objective evaluates to

\begin{equation} w^\top S w = w^\top (\lambda w) = \lambda. \end{equation}

so the achievable variances are exactly the eigenvalues, and the maximiser is the top eigenvector \(v_1\): the first principal component, carrying variance \(\lambda_1\) (Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon, 2020).

the rest, by deflation#

the second component solves the same problem with the extra constraint \(w \perp v_1\). restricting \(S\) to the subspace orthogonal to \(v_1\) leaves a symmetric matrix whose top eigenvector is \(v_2\), and by induction the \(k\)-th component is the \(k\)-th eigenvector. 𐃏 collecting the top \(k\) as columns of \(W_k = [v_1 \cdots v_k]\), the scores \(z_i = W_k^\top x_i \in \mathbb{R}^k\) are the low-dimensional code, with the tidy properties that the coordinates of \(z\) are uncorrelated and \(\operatorname{Var}(z_{(j)}) = \lambda_j\).

the reconstruction view#

now the other question: among all \(k\)-dimensional subspaces, which one loses the least when we project onto it? for an orthonormal \(W \in \mathbb{R}^{d \times k}\) (so \(W^\top W = I_k\)), the projection of \(x\) is \(W W^\top x\) and the average squared reconstruction error is

\begin{align*} J(W) &= \frac{1}{n}\sum_{i=1}^{n} \big\lVert x_i - W W^\top x_i \big\rVert^2 \\ &= \frac{1}{n}\sum_{i=1}^{n} \left( \lVert x_i \rVert^2 - x_i^\top W W^\top x_i \right) \\ &= \operatorname{tr}(S) - \operatorname{tr}\!\big(W^\top S W\big), \end{align*}

using \(\lVert x - WW^\top x\rVert^2 = x^\top x - x^\top W W^\top x\) (pythagoras: projection and residual are orthogonal). the first term is fixed, so minimising reconstruction error is exactly maximising the projected variance \(\operatorname{tr}(W^\top S W)\) — the two views coincide. the trace is maximised over orthonormal \(W\) by the top-\(k\) eigenvectors,¹ giving optimal error

\begin{equation} J(W_k^\ast) = \sum_{j=k+1}^{d} \lambda_j, \end{equation}

the sum of the discarded eigenvalues. small trailing eigenvalues mean the cloud genuinely lives near a \(k\)-flat, and pca finds it (Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon, 2020).

computation via the svd#

stack the centred data as rows of \(X \in \mathbb{R}^{n \times d}\) and take the (thin) singular value decomposition

\begin{equation} X = U \Sigma V^\top, \qquad U \in \mathbb{R}^{n \times r},\; \Sigma = \operatorname{diag}(\sigma_1, \dots, \sigma_r),\; V \in \mathbb{R}^{d \times r}. \end{equation}

then \(S = \frac{1}{n} X^\top X = V \frac{\Sigma^2}{n} V^\top\), so

• the right singular vectors \(V\) are the principal components,
• eigenvalues come from singular values as \(\lambda_j = \sigma_j^2 / n\),
• the score matrix is \(Z = X V = U \Sigma\) — no second pass over the data.

in practice you compute pca via the svd of \(X\), not by forming \(S\) and calling an eigensolver: forming \(X^\top X\) squares the condition number, so half your significant digits evaporate before the eigensolver even starts. 𐃏 when \(d \gg n\) (genomics, images) the trick inverts: eigendecompose the \(n \times n\) gram matrix \(\frac{1}{n} X X^\top\) instead, which shares the nonzero spectrum, then map back — this is also the doorway to kernel pca below.

how many components?#

the eigenvalue spectrum is the budget sheet. define the proportion of variance explained by the top \(k\):

\begin{equation} \mathrm{pve}(k) = \frac{\sum_{j=1}^{k} \lambda_j}{\sum_{j=1}^{d} \lambda_j}. \end{equation}

• scree plot: plot \(\lambda_j\) against \(j\) and look for the elbow where signal gives way to noise floor.
• threshold: keep the smallest \(k\) with \(\mathrm{pve}(k) \ge 0.9\) (or whatever your downstream task forgives).
• scale first: pca is not scale-invariant. a feature measured in millimetres will dominate one measured in kilometres for no honest reason — standardise features (pca on the correlation matrix) unless their units are genuinely commensurate (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).

whitening#

pca decorrelates; whitening goes one step further and equalises. with \(\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_k)\), the pca-whitened code

\begin{equation} z_i = \Lambda^{-1/2}\, W_k^\top x_i \end{equation}

has identity covariance: every direction contributes unit variance. useful before algorithms that assume isotropy (nearest neighbours, some optimisers). zca whitening \(z_i = W \Lambda^{-1/2} W^\top x_i\) (full \(W\), rotated back into the original basis) is the whitening transform closest to the identity — images stay looking like images. caveat: dividing by \(\sqrt{\lambda_j}\) for tiny \(\lambda_j\) amplifies noise \(1/\sqrt{\lambda_j}\)-fold, so regularise with \(\Lambda + \varepsilon I\) or truncate first.

kernel pca, sketched#

pca can only find flat structure. kernel pca runs the gram-matrix formulation in a feature space \(\varphi(x)\) you never compute, touching data only through \(k(x, x’) = \varphi(x)^\top \varphi(x’)\) (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009):

• centre the kernel matrix: \(\tilde{K} = H K H\) with \(H = I - \tfrac{1}{n}\mathbf{1}\mathbf{1}^\top\) (feature-space centring done implicitly).
• eigendecompose \(\tilde{K} \alpha_j = n \lambda_j \alpha_j\); normalise so \(\lambda_j\, n\, \alpha_j^\top \alpha_j = 1\), which makes the implicit feature-space component unit length.
• project a point via kernel evaluations only: \(z_{(j)}(x) = \sum_{i=1}^{n} \alpha_{j,i}\, \tilde{k}(x_i, x)\).

with an rbf kernel this unrolls spirals and separates concentric rings that vanilla pca is congenitally blind to. costs: \(O(n^2)\) memory, \(O(n^3)\) eigendecomposition, and no cheap inverse map from code space back to inputs (the “pre-image problem”). for visualisation-grade nonlinear reduction the modern defaults are t-sne and umap (McInnes, Leland and Healy, John and Melville, James, 2020), which optimise neighbourhood preservation rather than variance — see also kernel methods for the general kernel trick, and the curse of dimensionality for why any of this is necessary.²

from scratch#

both routes — eigendecomposition of the covariance and svd of the centred data — on a synthetic 5-d dataset with a planted 2-d signal subspace plus isotropic noise:

import numpy as np

rng = np.random.default_rng(42)

# 200 points in 5-d with a planted 2-d signal subspace
n, d = 200, 5
latent = rng.normal(size=(n, 2)) @ np.diag([3.0, 1.5])   # strong directions
basis, _ = np.linalg.qr(rng.normal(size=(d, 2)))          # random orthonormal 5x2
X = latent @ basis.T + 0.3 * rng.normal(size=(n, d))      # + isotropic noise

def pca_eig(X, k):
    Xc = X - X.mean(axis=0)
    S = Xc.T @ Xc / (len(X) - 1)            # sample covariance
    lam, W = np.linalg.eigh(S)              # ascending order
    lam, W = lam[::-1], W[:, ::-1]          # descending
    return lam[:k], W[:, :k], Xc

def pca_svd(X, k):
    Xc = X - X.mean(axis=0)
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    lam = s**2 / (len(X) - 1)               # singular values -> eigenvalues
    return lam[:k], Vt[:k].T, Xc

lam_e, W_e, Xc = pca_eig(X, 5)
lam_s, W_s, _ = pca_svd(X, 5)

print("eigenvalues (eig):", np.round(lam_e, 4))
print("eigenvalues (svd):", np.round(lam_s, 4))
# columns may differ by sign; align before comparing
signs = np.sign(np.sum(W_e * W_s, axis=0))
print("max |W_eig - W_svd| after sign fix:",
      f"{np.abs(W_e - W_s * signs).max():.2e}")

pve = lam_e / lam_e.sum()
print("proportion of variance explained:", np.round(pve, 4))
print("cumulative:                      ", np.round(np.cumsum(pve), 4))

# reconstruction error for k = 1..5 equals the sum of discarded eigenvalues
for k in range(1, 6):
    Z = Xc @ W_e[:, :k]                     # project
    err = np.mean(np.sum((Xc - Z @ W_e[:, :k].T)**2, axis=1))
    print(f"k={k}: mean sq reconstruction error {err:.4f}   "
          f"sum of discarded eigenvalues {(lam_e[k:].sum() * (n-1)/n):.4f}")

eigenvalues (eig): [8.9681 1.958  0.1099 0.0868 0.072 ]
eigenvalues (svd): [8.9681 1.958  0.1099 0.0868 0.072 ]
max |W_eig - W_svd| after sign fix: 9.55e-15
proportion of variance explained: [0.8011 0.1749 0.0098 0.0078 0.0064]
cumulative:                       [0.8011 0.976  0.9858 0.9936 1.    ]
k=1: mean sq reconstruction error 2.2156   sum of discarded eigenvalues 2.2156
k=2: mean sq reconstruction error 0.2674   sum of discarded eigenvalues 0.2674
k=3: mean sq reconstruction error 0.1580   sum of discarded eigenvalues 0.1580
k=4: mean sq reconstruction error 0.0716   sum of discarded eigenvalues 0.0716
k=5: mean sq reconstruction error 0.0000   sum of discarded eigenvalues 0.0000

everything the theory promised, numerically: the two factorisations agree to machine precision (up to sign — eigenvectors are only defined up to \(\pm 1\)), the planted 2-d subspace announces itself as two dominant eigenvalues carrying \(97.6\%\) of the variance, and the reconstruction error at each \(k\) equals the discarded eigenvalue mass exactly (the \((n-1)/n\) factor in the check just reconciles the two variance conventions).

see also#

• k-means clustering — the other half of ch. 14: compress rows instead of columns
• gaussian mixture models — probabilistic latent structure rather than geometric
• autoencoders — pca with the linearity constraint deleted
• kernel methods — the machinery behind kernel pca
• curse of dimensionality — why high-dimensional clouds need flattening at all