Part IV — Unsupervised Learning · Lecture 8

Principal Component Analysis

Three completely different questions — which direction has most variance, which subspace loses least information, and what does the SVD of the data matrix say — all have the same answer. That coincidence is PCA, and it is the most useful thing the spectral theorem does.

By the end you can

Derive PCA as variance maximisation using the Rayleigh quotient
Derive it again as minimum reconstruction error, and prove the two agree
Connect both to the SVD of the centred data matrix
Choose \(k\) from the spectrum, and know why 'explained variance' can mislead
State precisely what PCA cannot do

Assumed

Phase A, Lecture 11: spectral theorem, Rayleigh quotient, quadratic forms
Phase A, Lecture 8: projection onto a subspace
Phase A, Lecture 13: variance and covariance

1The setup, and why centring is not optional

Data \(\vx_1,\dots,\vx_n\in\RR^d\), stacked as rows of \(\mX\in\RR^{n\times d}\). First, centre: subtract the mean of each column, giving \(\mX_c\) with \(\frac1n\sum_i\vx_i=\vzero\). The sample covariance is then

\[\mS = \frac{1}{n}\mX_c^\top\mX_c \in \RR^{d\times d},\]

symmetric and positive semi-definite, so the spectral theorem applies: \(\mS=\mQ\Lambda\mQ^\top\) with orthonormal eigenvectors and non-negative eigenvalues (Phase A Lecture 11).

Common trap

Centring is not cosmetic. Without it, the largest "principal component" of a dataset whose mean is far from the origin will simply point at the mean, because \(\frac1n\|\mX\vu\|^2\) then measures second moment, not variance. You would spend your first component describing where the data is rather than how it varies. Every PCA implementation centres; if you write your own, centre first.

2Derivation 1 — maximum variance

Question: which unit direction \(\vu\) makes the projected data most spread out?

The top eigenvector maximises projected variance0 / 6

\[z_i = \vu^\top\vx_i \quad\text{(the projection of point } i \text{ onto direction } \vu)\]

Project. For a unit \(\vu\), \(\vu^\top\vx_i\) is the signed distance along \(\vu\) — Phase A Lecture 8's projection formula, with \(\|\vu\|=1\) so the denominator disappears. Since the data is centred, \(\bar z = \vu^\top\bar{\vx}=0\).

\[\Var(z) = \frac1n\sum_{i=1}^n z_i^2 = \frac1n\sum_i \vu^\top\vx_i\vx_i^\top\vu = \vu^\top\!\left(\frac1n\mX_c^\top\mX_c\right)\!\vu = \vu^\top\mS\vu\]

The projected variance is a quadratic form. Because the mean is zero, variance is just the mean square. Pulling \(\vu\) outside the sum turns the whole thing into \(\vu^\top\mS\vu\) — precisely Phase A Lecture 11's object.

\[\max_{\|\vu\|=1} \vu^\top\mS\vu \quad\Longleftrightarrow\quad \max_{\vu\neq\vzero}\ \frac{\vu^\top\mS\vu}{\vu^\top\vu}\]

The Rayleigh quotient. The constrained and unconstrained-ratio forms are equivalent because the ratio is invariant to rescaling \(\vu\), so we may as well normalise. Note the constraint is necessary: without it we could scale \(\vu\) up and make the variance arbitrarily large.

\[\vu^\top\mS\vu = \sum_{i=1}^d \lambda_i (\vq_i^\top\vu)^2 \le \lambda_1\sum_i(\vq_i^\top\vu)^2 = \lambda_1\|\vu\|^2 = \lambda_1\]

Expand in the eigenbasis. Write \(\vu=\sum_i c_i\vq_i\) with \(c_i=\vq_i^\top\vu\); then \(\mS\vu=\sum_i\lambda_ic_i\vq_i\) and the quadratic form is \(\sum_i\lambda_ic_i^2\). Bounding every \(\lambda_i\) by the largest gives the inequality; \(\sum_i c_i^2=\|\vu\|^2=1\) because \(\mQ\) is orthonormal (Parseval).

\[\boxed{\ \vu^\star = \vq_1, \qquad \max \Var = \lambda_1\ }\]

Equality is attained. Setting \(\vu=\vq_1\) gives \(c_1=1\) and all other \(c_i=0\), so the sum equals \(\lambda_1\) exactly. The top eigenvector of the covariance matrix is the direction of maximum variance, and its eigenvalue is that variance. ∎

\[\vu_2^\star = \vq_2 \quad\text{subject to } \vu_2\perp\vq_1, \qquad \dots\]

Subsequent components. Repeat the argument on the subspace orthogonal to \(\vq_1\): the same bound applies with \(\lambda_1\) removed from the sum, so the maximum is \(\lambda_2\), attained at \(\vq_2\). By induction, the \(k\)-th principal direction is \(\vq_k\). The orthogonality is free — the spectral theorem already guarantees the eigenvectors are perpendicular.

Sweep a direction and watch the projected varianceInteractive

datayour direction ufirst principal componentsecond

direction of u (degrees)10.0

correlation in the data0.85

stretch along one axis1.60

Grey ticks along the orange line are the projections of every point onto your chosen direction. Sweep the angle and read the projected variance below. It is maximised exactly when the orange line lines up with the blue PC1, and minimised at the purple PC2 — 90 degrees away, as the spectral theorem guarantees. Set the correlation to 0 and the stretch to 1 and the variance becomes the same in every direction: PCA has nothing to say about isotropic data.

3Derivation 2 — minimum reconstruction error

A completely different question: which \(k\)-dimensional subspace, when we project the data onto it and back, loses the least? Remarkably, the same answer.

Best rank-k approximation, and why it is the same problem0 / 6

\[\text{Let } \mU_k=[\vu_1|\cdots|\vu_k] \text{ with } \mU_k^\top\mU_k=\mI_k. \quad \hat{\vx}_i = \mU_k\mU_k^\top\vx_i\]

Project and reconstruct. \(\mU_k^\top\vx_i\in\RR^k\) is the coordinate vector in the subspace; \(\mU_k\) maps it back to \(\RR^d\). The composite \(\mP=\mU_k\mU_k^\top\) is the projection matrix, and it is idempotent and symmetric — Phase A Lecture 10, exercise 2.

\[J(\mU_k) = \frac1n\sum_{i=1}^{n}\big\|\vx_i - \mU_k\mU_k^\top\vx_i\big\|^2\]

The reconstruction error. The average squared distance from each point to the subspace. We want to minimise this — a different objective from Derivation 1, which maximised something.

\[\|\vx_i\|^2 = \|\mP\vx_i\|^2 + \|\vx_i-\mP\vx_i\|^2\]

Pythagoras. The projection and the residual are orthogonal (Phase A Lecture 8), so their squared lengths add. This is the pivot of the whole argument.

\[J(\mU_k) = \underbrace{\frac1n\sum_i\|\vx_i\|^2}_{\text{constant}} - \underbrace{\frac1n\sum_i\|\mU_k^\top\vx_i\|^2}_{\text{captured variance}}\]

Rearrange. Using \(\|\mP\vx_i\|=\|\mU_k^\top\vx_i\|\) since \(\mU_k\) has orthonormal columns. The first term does not depend on our choice of subspace. So minimising reconstruction error is exactly maximising captured variance — the two questions were the same question all along, related by a fixed total.

\[\frac1n\sum_i\|\mU_k^\top\vx_i\|^2 = \tr\big(\mU_k^\top\mS\,\mU_k\big) = \sum_{j=1}^{k}\vu_j^\top\mS\vu_j\]

Write the captured variance as a trace. Each column contributes its own Rayleigh quotient, and by the cyclic property of the trace the sum is \(\tr(\mU_k^\top\mS\mU_k)\).

\[\boxed{\ \max_{\mU_k^\top\mU_k=\mI}\ \tr(\mU_k^\top\mS\mU_k) = \sum_{j=1}^{k}\lambda_j, \quad\text{attained at } \mU_k=[\vq_1|\cdots|\vq_k]\ }\]

Ky Fan's theorem — the \(k\)-dimensional generalisation of the Rayleigh bound, proved by the same eigen-expansion argument. So the optimal subspace is spanned by the top \(k\) eigenvectors, the captured variance is \(\sum_{j\le k}\lambda_j\), and the residual error is \(\sum_{j>k}\lambda_j\). The discarded eigenvalues are literally the error you incur. ∎

4Derivation 3 — the SVD, and why you should use it

Every real matrix factors as \(\mX_c = \mU\mathbf{\Sigma}\mV^\top\) with \(\mU,\mV\) orthonormal and \(\mathbf{\Sigma}\) diagonal with non-negative entries \(\sigma_1\ge\sigma_2\ge\cdots\). Then

\[\mS = \frac1n\mX_c^\top\mX_c = \frac1n\mV\mathbf{\Sigma}^\top\mU^\top\mU\mathbf{\Sigma}\mV^\top = \mV\left(\frac{\mathbf{\Sigma}^2}{n}\right)\mV^\top,\]

using \(\mU^\top\mU=\mI\). Comparing with \(\mS=\mQ\Lambda\mQ^\top\) and matching the (unique, up to sign and ties) eigendecomposition:

the right singular vectors \(\vv_j\) are the principal directions \(\vq_j\);
the eigenvalues are \(\lambda_j = \sigma_j^2/n\);
the scores (projected coordinates) are \(\mX_c\mV = \mU\mathbf{\Sigma}\), so column \(j\) of \(\mU\mathbf{\Sigma}\) gives every point's coordinate on PC\(j\).

Practical engineering noteAlways compute PCA from the SVD of \(\mX_c\), never from the eigendecomposition of \(\mX_c^\top\mX_c\). Forming the Gram matrix squares the condition number (Phase A Lecture 10) — a design with \(\kappa=10^{8}\) gives \(\kappa=10^{16}\) for the covariance, which is all of float64. The SVD works on \(\mX_c\) directly. It is also cheaper when \(n\ll d\): the SVD costs \(O(nd\min(n,d))\) whereas forming \(\mS\) alone costs \(O(nd^2)\). Every serious implementation, including scikit-learn's, uses the SVD.

5Choosing k

The spectrum, the cumulative variance, and the reconstruction errorInteractive

eigenvalue λⱼ (scree)cumulative variance explainedyour k

components kept k3

true latent rank of the data4

noise level0.30

Thirty features generated from a low-dimensional latent signal plus noise. The scree plot shows a sharp elbow at the true rank, then a flat noise floor — that gap is the signal you look for. Raise the noise and the elbow blurs until it disappears entirely, which is exactly when 'choose k at the elbow' stops being usable advice and you should cross-validate on the downstream task instead.

Scree plot / elbow. Plot \(\lambda_j\) against \(j\) and look for the knee. Works when there is a genuine gap; degenerates into taste when there is not.
Cumulative variance. Choose the smallest \(k\) with \(\sum_{j\le k}\lambda_j/\sum_j \lambda_j \ge 0.9\) or \(0.95\). Fast, conventional, and arbitrary.
Cross-validate on the downstream task. If PCA feeds a classifier, tune \(k\) by the classifier's validation score. This is the only criterion that is actually about what you want, and it should be your default when a downstream task exists.

Check yourself

You run PCA on a dataset where one feature is measured in metres and another in millimetres. What happens?

Correct. Correct. Variance has units of the feature squared, so rescaling a feature by 1000 multiplies its variance by \(10^6\). PCA maximises variance and has no way to know that metres and millimetres are the same physical quantity — it will simply align PC1 with whichever column happens to carry the biggest numbers.

The fix is to standardise each feature to unit variance first, which makes \(\mS\) the correlation matrix. Then every feature contributes equally to the total variance \(\tr(\mS)=d\), and the components describe correlation structure rather than measurement units.

When NOT to standardise: if all features share a genuine common unit and their differing variances are real information. Pixel intensities, spectra, and repeated measurements of one quantity are the standard cases — there, a high-variance channel really is more informative, and standardising would throw that away by amplifying pure noise channels to the same scale as signal ones.

The choice is a modelling decision. Make it deliberately, and note it in your write-up.

Not quite. Normalising the eigenvectors does not normalise the data. PCA is invariant to rotations and translations of the data (both preserve the covariance structure up to a rotation of the eigenvectors) but emphatically not to per-feature rescaling, which changes \(\mS\) itself. Try it: scale one column of a 2-D dataset by 100 and watch PC1 swing to align with it.

Not quite. It will run perfectly happily and give you a mathematically correct answer to a question you did not mean to ask. That is the dangerous case — a loud failure would be easier to catch. The only numerical symptom is a very large condition number for \(\mS\), which is worth checking but is easy to miss.

6What PCA cannot do

It is linear. Data on a curved manifold — a spiral, a Swiss roll — has no low-dimensional linear subspace that captures it. PCA will report that you need many components when the intrinsic dimension is two. Kernel PCA, t-SNE, UMAP and autoencoders exist for this.
Variance is not information. The discriminative direction in a classification problem can easily be a low-variance one. PCA is blind to labels; if you have labels and want separation, use LDA or a supervised method, and do not assume the top components are the useful ones.
Components are not interpretable by default. PC1 is a dense linear combination of every feature. Reading it as "a factor" is a habit from factor analysis and requires care — the sign is arbitrary (Phase A Lecture 11), and a rotation of the top-\(k\) subspace explains exactly as much variance while looking completely different. Sparse PCA exists precisely to buy interpretability back.
It is sensitive to outliers, because variance is a squared quantity. One extreme point can define PC1 entirely. Robust PCA variants exist.

PCA from scratch, three ways, and a check that they agreePython · NumPy

import numpy as np
rng = np.random.default_rng(0)
np.set_printoptions(precision=4, suppress=True)

# ---------------- data with a known 3-dimensional latent structure -------------
n, d, k_true = 400, 12, 3
W = rng.normal(size=(k_true, d))
Zlat = rng.normal(size=(n, k_true)) * np.array([3.0, 2.0, 1.0])
X = Zlat @ W + 0.4*rng.normal(size=(n, d)) + np.array([5.0]*d)   # note the offset

Xc = X - X.mean(0)                      # CENTRE. always.
S = Xc.T @ Xc / n

# ---- method 1: eigendecomposition of the covariance ----
lam, Q = np.linalg.eigh(S)
o = np.argsort(lam)[::-1]; lam, Q = lam[o], Q[:, o]

# ---- method 2: SVD of the centred data ----
U, sv, Vt = np.linalg.svd(Xc, full_matrices=False)
lam_svd = sv**2 / n
V = Vt.T

print("eigenvalues (eigh):", lam[:6])
print("sigma^2/n  (svd) :", lam_svd[:6])
print("agree?", np.allclose(lam, lam_svd, atol=1e-9))
print("directions agree (up to sign)?",
      np.allclose(np.abs(Q[:, :5]), np.abs(V[:, :5]), atol=1e-8))
print()

# ---- method 3: maximise the Rayleigh quotient by brute force, in 2-D ----
print("--- BRUTE FORCE CHECK IN THE TOP 2-D SUBSPACE ---")
X2 = Xc @ Q[:, :2]                                   # work in the PC1-PC2 plane
S2 = X2.T @ X2 / n
best = max(((np.array([np.cos(t), np.sin(t)]) @ S2 @ np.array([np.cos(t), np.sin(t)]), t)
            for t in np.linspace(0, np.pi, 20001)))
print(f"  best projected variance over a fine angle sweep: {best[0]:.6f}")
print(f"  lambda_1 from the eigendecomposition:            {lam[0]:.6f}")
print(f"  angle found: {np.degrees(best[1]):.3f} deg  (0 deg = PC1, as expected)")
print()

# ---- reconstruction error IS the discarded eigenvalue sum ----
print("--- RECONSTRUCTION ERROR = SUM OF DISCARDED EIGENVALUES ---")
print(f"{'k':>3} {'explained':>11} {'measured err':>14} {'sum_{j>k} lam':>15}")
for kk in (1, 2, 3, 4, 6, 12):
    Uk = Q[:, :kk]
    recon = Xc @ Uk @ Uk.T
    err = float(np.mean(np.sum((Xc - recon)**2, axis=1)))
    print(f"{kk:3d} {lam[:kk].sum()/lam.sum():10.2%} {err:14.6f} {lam[kk:].sum():15.6f}")
print("  Exact, at every k. Eckart-Young in numbers.")
print()

# ---- PCA beats every random subspace ----
print("--- IS THE TOP-k SUBSPACE REALLY OPTIMAL? ---")
for kk in (1, 3, 5):
    Uk = Q[:, :kk]
    best_err = float(np.mean(np.sum((Xc - Xc@Uk@Uk.T)**2, axis=1)))
    rand_best = min(
        float(np.mean(np.sum((Xc - Xc@R@R.T)**2, axis=1)))
        for R in (np.linalg.qr(rng.normal(size=(d, kk)))[0] for _ in range(3000)))
    print(f"  k={kk}  PCA {best_err:10.5f}   best of 3000 random subspaces {rand_best:10.5f}"
          f"   {'PCA wins' if best_err <= rand_best + 1e-9 else 'PROBLEM'}")
print()

# ---- what happens if you forget to centre ----
print("--- FORGETTING TO CENTRE ---")
S_raw = X.T @ X / n
lr, Qr = np.linalg.eigh(S_raw); o = np.argsort(lr)[::-1]; lr, Qr = lr[o], Qr[:, o]
mean_dir = X.mean(0)/np.linalg.norm(X.mean(0))
print(f"  |cos angle between 'PC1' and the data mean| = {abs(Qr[:,0] @ mean_dir):.4f}")
print(f"  first eigenvalue: centred {lam[0]:.3f}   uncentred {lr[0]:.3f}")
print("  Uncentred, the top component just points at the mean. Useless.")
print()

# ---- whitening ----
print("--- WHITENING ---")
Zw = (Xc @ Q[:, :5]) / np.sqrt(lam[:5])
print("  covariance of the whitened scores (should be the identity):")
print(np.round(Zw.T @ Zw / n, 6))
print()

# ---- standardise or not? ----
print("--- SCALE SENSITIVITY ---")
Xs = X.copy(); Xs[:, 0] *= 1000.0
for tag, M in (("original", X), ("feature 0 x1000", Xs), ("standardised", (Xs-Xs.mean(0))/Xs.std(0))):
    Mc = M - M.mean(0)
    l, Qq = np.linalg.eigh(Mc.T@Mc/n); oo = np.argsort(l)[::-1]; l, Qq = l[oo], Qq[:, oo]
    print(f"  {tag:18s} PC1 loading on feature 0 = {abs(Qq[0,0]):.4f}"
          f"   explained by PC1 = {l[0]/l.sum():.2%}")
print("  Rescaling one feature hijacks PC1. Standardising restores balance.")
print()

# ---- PCA is blind to labels ----
print("--- VARIANCE IS NOT DISCRIMINATION ---")
lab = rng.integers(0, 2, n)
Xd = rng.normal(size=(n, 2)) * np.array([8.0, 1.0])       # huge variance on axis 0
Xd[:, 1] += 3.0*lab                                        # but the LABEL lives on axis 1
Xdc = Xd - Xd.mean(0)
ld, Qd = np.linalg.eigh(Xdc.T@Xdc/n); od = np.argsort(ld)[::-1]; ld, Qd = ld[od], Qd[:, od]
for j in (0, 1):
    s = Xdc @ Qd[:, j]
    sep = abs(s[lab==1].mean() - s[lab==0].mean()) / s.std()
    print(f"  PC{j+1}: {ld[j]/ld.sum():6.1%} of variance,"
          f"   class separation {sep:.3f} standard deviations")
print("  PC1 holds almost all the variance and almost none of the signal.")
print("  If you reduced to k=1 with PCA here, you would throw the label away.")

7Exercises

★

A covariance matrix has eigenvalues \(\{40, 25, 20, 10, 4, 1\}\).

What fraction of variance do the first two components explain?
How many components for 90%?
What is the reconstruction error at \(k=3\), and what is the total variance?
Would you call this a low-dimensional dataset?

1. Total \(=40+25+20+10+4+1=100\). First two: \(65/100 = \mathbf{65\%}\).

2. Cumulative: 40, 65, 85, 95, 99, 100. So \(k=4\) reaches 95%, and \(k=3\) only reaches 85%. Four components are needed for 90%.

3. Reconstruction error at \(k=3\) is \(\sum_{j>3}\lambda_j = 10+4+1 = \mathbf{15}\), which is 15% of the total variance \(\tr(\mS)=\mathbf{100}\). (Total variance is the trace, which equals the sum of the eigenvalues — Phase A Lecture 11.)

4. No, not really. The spectrum decays slowly: 40, 25, 20 are all comparable, and there is no sharp elbow. Genuinely low-dimensional data looks like \(\{80, 15, 2, 1, 1, 1\}\) — a big gap and then a flat floor. Here you need two-thirds of the components for 95% of the variance, which is barely a reduction.

What you can say: there is some structure (the top component is 40× the bottom one), and dropping the last two costs 5%, which may well be worth it for denoising. But do not describe this as "the data is really 2-dimensional".

★★

Prove that the two derivations give the same answer, by showing directly that

\[ \argmin_{\mU_k}\ \frac1n\sum_i\|\vx_i-\mU_k\mU_k^\top\vx_i\|^2 \;=\; \argmax_{\mU_k}\ \tr(\mU_k^\top\mS\mU_k). \]

Then explain why \(\tr(\mS)=\sum_j\lambda_j\) is the "total variance", and what that means for the interpretation of the explained-variance ratio.

Equivalence. Let \(\mP=\mU_k\mU_k^\top\), which is a projection: \(\mP^2=\mP\) and \(\mP^\top=\mP\). For each point, the projection and residual are orthogonal:

\[(\mP\vx)^\top(\vx-\mP\vx) = \vx^\top\mP\vx - \vx^\top\mP^\top\mP\vx = \vx^\top\mP\vx-\vx^\top\mP\vx = 0,\]

using \(\mP^\top\mP=\mP^2=\mP\). So Pythagoras gives \(\|\vx\|^2=\|\mP\vx\|^2+\|\vx-\mP\vx\|^2\). Averaging over \(i\):

\[\underbrace{\frac1n\sum_i\|\vx_i\|^2}_{=\,\tr(\mS),\ \text{fixed}} = \frac1n\sum_i\|\mP\vx_i\|^2 + J(\mU_k).\]

The left-hand side does not involve \(\mU_k\), so minimising \(J\) is exactly maximising the middle term. And \(\|\mP\vx\|^2 = \vx^\top\mU_k\mU_k^\top\vx = \|\mU_k^\top\vx\|^2\), whose average is \(\tr(\mU_k^\top\mS\mU_k)\) by the cyclic property of the trace. ∎

Why \(\tr(\mS)\) is the total variance. The diagonal entries \(S_{jj}\) are the individual feature variances, so \(\tr(\mS)=\sum_j\Var(x_j)\). And the trace equals the sum of the eigenvalues (Phase A Lecture 11's \(\lambda_1+\lambda_2=\tr\mA\), in general dimension). So

\[\sum_{j}\Var(x_j) = \sum_j \lambda_j.\]

The total variance is invariant under rotation — PCA redistributes it among orthogonal directions without creating or destroying any. That is what makes the ratio \(\sum_{j\le k}\lambda_j/\sum_j\lambda_j\) a meaningful "fraction explained": the denominator is a conserved quantity.

The interpretive caveat. It is a fraction of variance, not a fraction of information, and the two come apart badly in at least three ways. If you standardised, every feature contributes exactly 1 to the total, so \(\tr(\mS)=d\) and the ratio partly reflects your standardisation choice. If a feature is pure noise with large amplitude, it inflates the denominator and makes every ratio look worse. And if the signal you care about is a low-variance direction (the last figure in the code block), a component explaining 2% of the variance may carry 100% of what matters. "95% of variance explained" is a statement about geometry, never about utility.

★★

Whitening. Define \(\vz_i = \Lambda_k^{-1/2}\mQ_k^\top\vx_i\).

Show the whitened data has identity covariance.
Why is whitening useful before another algorithm? Give two concrete reasons using Phase A results.
Why is whitening dangerous when the spectrum has very small eigenvalues?

1. The covariance of the whitened scores is

\[\frac1n\sum_i \vz_i\vz_i^\top = \Lambda_k^{-1/2}\mQ_k^\top\left(\frac1n\sum_i\vx_i\vx_i^\top\right)\mQ_k\Lambda_k^{-1/2} = \Lambda_k^{-1/2}\big(\mQ_k^\top\mS\mQ_k\big)\Lambda_k^{-1/2}.\]

Since the columns of \(\mQ_k\) are eigenvectors, \(\mS\mQ_k=\mQ_k\Lambda_k\), so \(\mQ_k^\top\mS\mQ_k=\Lambda_k\). Hence the covariance is \(\Lambda_k^{-1/2}\Lambda_k\Lambda_k^{-1/2}=\mI_k\). ∎ (\(\Lambda_k^{-1/2}\) exists because the retained eigenvalues are positive.)

2. Two reasons, both from Phase A.

Conditioning. After whitening the covariance is \(\mI\), so \(\kappa=1\). By Lecture 03, gradient descent then converges in \(O(1)\) iterations instead of \(O(\kappa)\), and the zig-zag disappears entirely. Whitening is the strongest possible preconditioner for a quadratic problem — it is the Newton step, applied to the data instead of the gradient.

Isotropic distance. Euclidean distance in whitened space equals Mahalanobis distance in the original space (Phase A Lecture 14, exercise 3). Any algorithm that uses raw distances — \(k\)-means, \(k\)-NN, RBF kernels — is implicitly assuming isotropy, and whitening makes that assumption true.

3. The danger. Whitening divides by \(\sqrt{\lambda_j}\). A direction with \(\lambda_j\approx0\) is one in which the data barely varies at all — which usually means it is pure measurement noise. Dividing by a tiny number amplifies that noise to the same scale as the genuine signal, and can amplify it by orders of magnitude.

Concretely, with \(\lambda_{\min}=10^{-6}\) and \(\lambda_{\max}=1\), whitening multiplies the noisiest direction by 1000 relative to the strongest. You have "improved" the conditioning from \(10^6\) to \(1\) by promoting noise to first-class status.

The standard fixes: whiten only the top \(k\) components (discarding the rest, as the formula above does); or use \(\Lambda^{-1/2}\) with a floor, \((\Lambda+\epsilon\mI)^{-1/2}\), which is ridge regularisation applied to whitening. This is exactly the same trade-off as Lecture 06's shrinkage factors, viewed from the other side — there we shrank small-eigenvalue directions, here we must be careful not to inflate them.

★★★

Coding. Investigate when PCA helps a downstream classifier and when it hurts.

Build a dataset in \(\RR^{50}\) where the class signal lies in a high-variance direction. Fit logistic regression on the top \(k\) components for \(k=1,\dots,50\) and plot test accuracy against \(k\).
Now build one where the signal lies in a low-variance direction. Repeat.
Compare against LDA, which uses the labels, on both.
State the rule for when unsupervised dimensionality reduction is safe.

import numpy as np
rng = np.random.default_rng(0)
sig = lambda z: 1/(1+np.exp(-np.clip(z,-500,500)))

def logistic(X, y, T=3000, lr=0.3, lam=1e-3):
    Xb = np.column_stack([np.ones(len(y)), X]); w = np.zeros(Xb.shape[1])
    for _ in range(T):
        g = Xb.T @ (sig(Xb@w) - y) / len(y); g[1:] += lam*w[1:]
        w -= lr*g
    return w
def acc(X, y, w):
    return float((( sig(np.column_stack([np.ones(len(y)),X])@w) >= .5) == y).mean())

def make(kind, n=1500, d=50, seed=0):
    r = np.random.default_rng(seed)
    y = r.integers(0, 2, n).astype(float)
    scales = np.logspace(1, -1, d)                 # variance 10 down to 0.1
    X = r.normal(size=(n, d)) * scales
    j = 0 if kind == "high" else d-1               # which direction carries the label
    X[:, j] += 2.0*scales[j]*(2*y-1)
    return X, y

for kind in ("high", "low"):
    Xtr, ytr = make(kind, seed=1); Xte, yte = make(kind, seed=2)
    mu = Xtr.mean(0); Zc = Xtr - mu
    l, Q = np.linalg.eigh(Zc.T@Zc/len(Zc)); o = np.argsort(l)[::-1]; l, Q = l[o], Q[:,o]
    row = []
    for k in (1, 2, 5, 10, 25, 50):
        w = logistic(Zc@Q[:,:k], ytr)
        row.append((k, acc((Xte-mu)@Q[:,:k], yte, w), l[:k].sum()/l.sum()))
    print(f"signal in a {kind}-variance direction:")
    for k, a, ev in row:
        print(f"   k={k:3d}  test acc {a:.4f}   variance explained {ev:6.1%}")
    # LDA direction: Sigma^-1 (mu1 - mu0)
    m1, m0 = Xtr[ytr==1].mean(0), Xtr[ytr==0].mean(0)
    Sw = np.cov(Xtr.T) + 1e-6*np.eye(Xtr.shape[1])
    v = np.linalg.solve(Sw, m1-m0); v /= np.linalg.norm(v)
    w = logistic((Xtr@v)[:,None], ytr)
    print(f"   LDA, k=1               test acc {acc((Xte@v)[:,None], yte, w):.4f}\\n")

1. Signal in a high-variance direction. Accuracy is already near its ceiling at \(k=1\) and stays flat as \(k\) grows. PCA has done exactly what you hoped: reduced 50 dimensions to 1 with no loss. Note that \(k=1\) explains only a modest share of total variance and yet captures all the signal — variance explained and accuracy are not the same axis.

2. Signal in a low-variance direction. Accuracy at \(k=1\) is around chance, and stays near chance until \(k\) approaches 50, at which point it jumps. The discriminative direction is the last principal component, so any aggressive PCA truncation destroys the model. If you had followed "keep 95% of the variance" you would have discarded the only useful direction and concluded the problem was unlearnable.

3. LDA gets near-ceiling accuracy with \(k=1\) in both cases, because \(\mS_w^{-1}(\vm_1-\vm_0)\) is constructed from the labels: it finds the direction maximising between-class separation relative to within-class scatter, wherever that direction happens to sit in the spectrum. This is not magic — it is simply using information PCA is not given.

4. The rule. Unsupervised dimensionality reduction is safe exactly when the variance structure of the inputs is aligned with the structure of the target — and you cannot check that without looking at the target.

Practical consequences:

If you have labels, choose \(k\) by cross-validating the downstream metric, never by a variance threshold. The variance threshold is a heuristic for the case where you have no other criterion.
PCA's real strengths are noise reduction, decorrelation, conditioning, visualisation and compression — all of which are about the inputs alone and are not compromised by label-blindness.
If you specifically want a supervised projection, use LDA, partial least squares, or simply let a neural network learn its own first layer, which is a supervised dimensionality reduction.
Whatever you choose, fit PCA on the training split only — fitting it on all the data is Lecture 01's preprocessing leak.

Next: PCA finds directions. Lecture 09 finds groups — k-means, its probabilistic upgrade to Gaussian mixtures, and the EM algorithm that fits them, which is one of the most reused ideas in the subject.

Previous7 · Bias, Variance, and Honest Evaluation Next9 · k-means, Gaussian Mixtures, and EM