Week 14: PCA, the Fourier Matrix and the FFT

You have 240 measurements of three features, x, y and z, and want to keep two numbers each. Feature z has the smallest variance, 2.89 against 5.81 and 4.56, so you drop it. It throws away 21.8% of the variance. The cloud is thin, but not along any raw axis.

The amber line is PC1, the solid line PC2 and the dashed line PC3, each drawn two standard deviations either side of the centre. Slide $t$ from $0$ to $1$ and every point moves onto the PC1–PC2 plane, keeping 98.1% of the variance. From $1$ to $2$ the points move onto the PC1 line, which keeps 67.9%.

Principal components from the SVD#

PCAPrincipal Component Analysis is the statistical reading of Week 13's SVD. Collect $m$ observations of $n$ features into $X \in \mathbb{R}^{m \times n}$, rows as observations. Subtract each column's mean, so the cloud sits at the origin, and call the result $\tilde{X}$. Uncentred, the first direction is pulled toward the mean.

Take the SVD $\tilde{X} = U\Sigma V^T$. The principal components are the columns $v_i$ of $V$, ordered by decreasing $\sigma_i$.

Why maximum variance. For a unit vector $w$, the variance of the data along $w$ is

$C$ is the sample covariance. With $y = V^Tw$, $w^TCw = \sum_i \lambda_i y_i^2 \le \lambda_1$, so it is largest at the top eigenvector $v_1$, with value $\sigma_1^2/(m-1)$.

The scene's numbers#

The scene's cloud is built so that its covariance is exact:

with Week 13's axes $v_1 = (2, 2, 1)/3$, $v_2 = (-2, 1, 2)/3$ and $v_3 = (1, -2, 2)/3$. The total variance is the trace, and the trace does not change with the basis:

So PC1, PC2 and PC3 carry $9/13.25 = 67.9\%$, $30.2\%$ and $1.9\%$. The diagonal of $C$ gives the raw features $43.8\%$, $34.4\%$ and $21.8\%$.

Dropping a raw feature is not dropping a component. Each raw variance is a mix of all three eigenvalues. For $z$, $C_{33} = \tfrac{1}{9}(1 \cdot 9 + 4 \cdot 4 + 4 \cdot 0.25) = 2.89$, because $z$ leans on $v_2$ and $v_3$ alike. No choice of two raw features keeps more than $78.2\%$. The best plane keeps $98.1\%$, and it is oblique.

Truncation is projection. Keeping $k$ components replaces $\tilde{X}$ by $\tilde{X}_k = U_k\Sigma_kV_k^T$, which is the slider at $t = 1$ for $k = 2$. Each point moves to $\tilde{x} - (\tilde{x}\cdot v_3)v_3$, and the share of variance lost is $\sigma_3^2/\sum_j\sigma_j^2$. By Eckart–Young this is the best rank-$k$ approximation in both the Frobenius and the spectral norm, so no other plane keeps more.

The Fourier matrix and the FFT#

Week 13 ended with unitary matrices. The most used one is the Fourier matrix. With $\omega = e^{2\pi i/n}$, $F_{jk} = \omega^{jk}$ for $j, k = 0, \ldots, n-1$, and $\hat{x} = Fx$ is the discrete Fourier transform (DFTDiscrete Fourier Transform), $\hat{x}_j = \sum_k x_k\omega^{jk}$. NumPy and most libraries use the conjugate sign, $\overline{F}x$. For $n = 4$, $\omega = i$ and

Its columns are orthogonal: $F^HF = nI$, so $F/\sqrt{n}$ is unitary and $x = \tfrac{1}{n}F^H\hat{x}$.

The product $Fx$ costs $O(n^2)$. The FFTFast Fourier Transform splits $x$ into its even and odd entries, transforms each half at size $n/2$ (with $\omega^2$), and combines them, for $n$ a power of 2:

The cost obeys $T(n) = 2T(n/2) + O(n)$, so $T(n) = O(n\log n)$: $\log_2 n$ levels of $O(n)$ work.

Course recap#

• $A = LU$ (Weeks 3–4) turns $Ax = b$ into two triangular solves.
• $A = QR$ (Weeks 7–8) solves least squares without the normal equations.
• $A = S\Lambda S^{-1}$ and $A = Q\Lambda Q^T$ (Weeks 9–12) give powers and quadratic forms.
• $A = U\Sigma V^T$ (Weeks 13–14) works for every matrix, and for centred data its $V$ and $\sigma_i$ are PCA.
• $F/\sqrt{n}$ is a unitary change of basis, made cheap by the FFT.

Knowledge Check#

Browser lab: PCA of the scene's cloud#

Rebuild the scene's cloud, run PCA, and check the FFT.

import numpy as np
import matplotlib.pyplot as plt

V = np.array([[2, 2, 1], [-2, 1, 2], [1, -2, 2]]).T / 3   # columns v1, v2, v3
sigma = np.array([3.0, 2.0, 0.5])

rng = np.random.default_rng(14)
Z = rng.normal(size=(240, 3))
Z -= Z.mean(axis=0)
L = np.linalg.cholesky(np.cov(Z, rowvar=False))
Z = Z @ np.linalg.inv(L).T            # sample covariance exactly I
X = (Z * sigma) @ V.T                 # sample covariance exactly V diag(9, 4, 0.25) Vᵀ

Xc = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
var = S**2 / (len(X) - 1)
total = var.sum()
print("PC variances  ", var.round(4), "shares", (var / total).round(3))
print("raw variances ", Xc.var(axis=0, ddof=1).round(4))
print("keep x, y      ", round(1 - Xc[:, 2].var(ddof=1) / total, 3))
print("keep PC1, PC2  ", round((var[0] + var[1]) / total, 3))
print("Vᵀ rows       ", Vt.round(3).tolist())

scores = Xc @ Vt.T
fig, axes = plt.subplots(1, 2, figsize=(10, 4.2))
axes[0].scatter(Xc[:, 0], Xc[:, 1], s=8)
axes[0].set_title("raw features x, y")
axes[1].scatter(scores[:, 0], scores[:, 1], s=8)
axes[1].set_title("scores on PC1, PC2")
for ax in axes:
    ax.set_aspect("equal")
plt.tight_layout()
plt.show()

n = 8
F = np.exp(2j * np.pi * np.outer(np.arange(n), np.arange(n)) / n)
x = rng.normal(size=n)
print("F^H F = nI:", np.allclose(F.conj().T @ F, n * np.eye(n)))
print("fft matches F̄x:", np.allclose(np.fft.fft(x), F.conj() @ x))

The PC variances print as $9$, $4$ and $0.25$. Keeping $x$ and $y$ keeps $0.782$, against $0.981$. The printed rows are $v_1$, $v_2$ and $v_3$, possibly with a sign flipped, as in Week 13. The FFT check matches $\overline{F}x$, NumPy's sign. Change sigma to [3, 2, 1] and predict the gap before you run it.

Key takeaways#

PCA is the SVD of the centred data: the components are the $v_i$, and the variance along $v_i$ is $\sigma_i^2/(m-1)$. A thin cloud need not be thin along a raw axis: in the scene, dropping $z$ keeps 78.2% and dropping PC3 keeps 98.1%. The Fourier matrix is unitary up to $\sqrt{n}$, and the FFT applies it in $O(n\log n)$.

Further reading#

• Gilbert Strang, Introduction to Linear Algebra (6th ed., Wellesley-Cambridge Press, 2023). Covers PCA through the SVD, and the Fourier matrix and the FFT.
• Lloyd N. Trefethen & David Bau, Numerical Linear Algebra (SIAM, 1997). The SVD and the Eckart–Young theorem behind PCA as low-rank approximation.
• Steven L. Brunton & J. Nathan Kutz, Data-Driven Science and Engineering (Cambridge University Press, 2019). PCA for dimensionality reduction and the FFT for signal processing.