Week 14: The Fourier Matrix, FFT, and PCA

The Fourier Matrix#

Week 13 closed with unitary matrices — complex analogues of orthogonal matrices that preserve length via $U^H U = I$. This week opens with the most important unitary family in signal processing and numerical linear algebra: the Fourier matrix.

Fix an integer $n \ge 1$ and let

be a primitive $n$-th root of unity (so $\omega^n = 1$ and $\omega^k \neq 1$ for $0 < k < n$). The Fourier matrix $F \in \mathbb{C}^{n \times n}$ has entries

(using $0$-based indexing):

Multiplying a vector $x$ by $F$ produces the discrete Fourier transform (DFTDiscrete Fourier Transform):

Each frequency bin $\hat{x}_j$ is the inner product of $x$ against a pure complex exponential of frequency $j$. The inverse transform recovers $x$ from $\hat{x}$ by conjugation (up to scaling): $x = \frac{1}{n} \overline{F}\,\hat{x}$.

Unitary (up to scaling). The columns of $F$ are the discrete complex exponentials

These columns are orthogonal with respect to the complex inner product: $f_j^H f_k = n$ when $j \equiv k \pmod n$ and $0$ otherwise. Equivalently,

Scaling by $1/\sqrt{n}$ produces a truly unitary matrix $\frac{1}{\sqrt{n}} F$ with $\bigl(\frac{1}{\sqrt{n}} F\bigr)^H\bigl(\frac{1}{\sqrt{n}} F\bigr) = I$. So $F$ is unitary up to that standard factor: a complex change of basis.

Circulant connection. A circulant matrix is constant on every wrap-around diagonal: each row is a right cyclic shift of the row above. Let $c_m$ be the entry on the diagonal at cyclic offset $m$, so $C_{jl} = c_{(j-l) \bmod n}$. The action on a Fourier column is then a scalar multiple:

a multiple of the $j$-th entry of $f_k$. So every $n \times n$ circulant matrix is diagonalized by the same $F$, with eigenvalue the DFT of its first row. That is why convolutions (which become multiplications after a DFT) and many shift-invariant filters are natural in the Fourier basis.

Tiny example. For $n = 2$, $\omega = e^{\pi i} = -1$, so

For $n = 4$, $\omega = i$ and the second column is $(1, i, -1, -i)^T$. Direct matrix–vector multiplication costs $O(n^2)$ arithmetic — fine for $n = 4$, painful for $n = 10^6$.

The Fast Fourier Transform#

A naive product $Fx$ walks every entry of $F$ and every component of $x$, which is $\Theta(n^2)$ work. The FFTFast Fourier Transform (Fast Fourier Transform) is a family of algorithms that compute the same DFT in $O(n\log n)$ operations when $n$ is a power of two (and with similar near-linear cost for highly composite $n$).

The key is a divide-and-conquer factorization of $F_n$. Split the input into even and odd indices:

Each of those half-length vectors has a DFT of size $n/2$. The Cooley–Tukey identity says the full DFT can be assembled from the two half-size DFTs with only $O(n)$ extra multiplications by “twiddle factors” $\omega^j$:

for $j = 0,\ldots,n/2-1$ (indices of the half-size transforms taken modulo $n/2$). In matrix language, $F_n$ factors into a product involving a permutation (even/odd reordering), a block-diagonal matrix with two copies of $F_{n/2}$, and a sparse diagonal of twiddles glued by an $n \times n$ “butterfly” step.

Complexity. Let $T(n)$ be the cost of an $n$-point FFT. The recurrence is

Unrolling gives $T(n) = O(n\log n)$ — each of $\log_2 n$ levels costs $O(n)$ work across the whole tree. Compared with the $O(n^2)$ direct matrix product, that is the difference between “impossible at audio or image scale” and “routine in real time.”

The butterflies are an implementation detail; what to carry away is the recursion, and libraries such as numpy.fft.fft implement it in carefully optimized variants — including more general factorizations when $n$ is not a power of two.

Principal Component Analysis#

PCAPrincipal Component Analysis (Principal Component Analysis) is the statistical reading of the SVD from Week 13: for a cloud of points in feature space it returns the orthonormal directions along which the cloud varies the most.

Setup. Collect $m$ observations of $n$ features into a data matrix $X \in \mathbb{R}^{m \times n}$, with rows as observations and columns as features. (Images, sensor readings, word counts — any tabular cloud of points works.)

Center the data. Subtract the sample mean of each column so the cloud sits at the origin:

Without centering, the first “component” often just points toward the mean rather than along the spread of the data.

SVD of the centered matrix. Compute

with orthogonal $U$ ($m \times m$ or thin), diagonal singular values $\sigma_1 \ge \sigma_2 \ge \cdots \ge 0$, and orthogonal $V$ ($n \times n$). The principal components are the columns of $V$ — the right singular vectors — ordered by decreasing singular value. Column $v_1$ is the direction of greatest variance; $v_2$ is the next, orthogonal to $v_1$; and so on.

Why “maximum variance”? For a unit vector $w \in \mathbb{R}^n$, the sample variance of the projected (centered) data is

The matrix $\tilde{X}^T\tilde{X}$ is symmetric positive semidefinite. By the spectral theorem (Week 12) and the SVD link $\tilde{X}^T\tilde{X} = V\Sigma^T\Sigma V^T$, its leading eigenvector is $v_1$ and the maximum of $\|\tilde{X} w\|^2$ over $\|w\|=1$ is $\sigma_1^2$. Each subsequent principal component solves the same problem in the orthogonal complement of the previous ones.

Scores. The coordinates of the data in the principal-component basis are the rows of $\tilde{X}V = U\Sigma$. Keeping only the first $k$ columns of $V$ projects every observation onto a $k$-dimensional subspace spanned by the top principal components.

Variance and Dimensionality Reduction#

Singular values quantify how much structure each principal component captures. Because

(with $r = \operatorname{rank}(\tilde{X})$), the fraction of variance explained by principal component $i$ is

The cumulative fraction for the first $k$ components is $\bigl(\sum_{i=1}^{k}\sigma_i^2\bigr)/\bigl(\sum_j\sigma_j^2\bigr)$. If that number is $0.95$, the top $k$ directions retain 95% of the total (centered) energy of the data.

Low-rank approximation. Truncating the SVD to the top $k$ singular values produces

the best rank-$k$ approximation of $\tilde{X}$ in both Frobenius and spectral norms (Eckart–Young, previewed in Week 13). PCA as dimensionality reduction is exactly this truncation: store the $n \times k$ matrix $V_k$ and the $m \times k$ scores $U_k\Sigma_k$, and reconstruct approximately when needed.

Practice sketch. Suppose centered 2-D data has singular values $\sigma_1 = 3$ and $\sigma_2 = 1$. Then

Ninety percent of the variance lies along PC1. Projecting onto that single axis compresses each point from two numbers to one while losing only 10% of the energy — often an excellent trade for visualization or noise filtering.

In high dimensions the same arithmetic applies: when the centered data matrix has rapidly decaying singular values, a handful of principal components can stand in for hundreds of raw features — and the Eckart–Young theorem (Week 13) guarantees that this truncation is the best rank-$k$ approximation of the data. The variance and covariance behind this calculation are developed in Probability & Statistics Week 3: Expectation, Variance, and Covariance, and PCA's role in learned representations continues in Generative Models Week 11: Representation Learning with Generative Models. The browser lab below runs PCA via np.linalg.svd on a synthetic correlated cloud and prints the explained-variance ratios.

Course Recap#

Fourteen weeks, one throughline: factorizations that turn hard linear problems into easy ones.

• Elimination ($LU$) — Weeks 2–4. Gaussian elimination factors $A = LU$ (with pivots and permutations when needed) so that $Ax = b$ becomes two triangular solves. The algebra of row operations becomes a product of elementary matrices.
• Orthogonalization ($QR$) — Weeks 7–8. Gram–Schmidt and Householder build $A = QR$ with orthonormal columns in $Q$. Least squares $\min\|Ax-b\|$ stabilizes into solving $Rx = Q^Tb$ instead of forming the often-ill-conditioned normal equations alone.
• Eigendecomposition — Weeks 9–12. $A = S\Lambda S^{-1}$ (and $A = Q\Lambda Q^T$ when $A$ is symmetric) explains powers $A^k$, differential systems, quadratic forms, and positive definiteness. Dynamics become independent modes along eigenvectors.
• The SVD — Week 13. $A = U\Sigma V^T$ works for every matrix, square or not. It unifies and generalizes the spectral picture: singular values replace eigenvalues when symmetry or squareness is missing.
• This week. The Fourier matrix is a structured unitary (up to scale) change of basis that diagonalizes every circulant operator; the FFT makes that change cheap. PCA reads the SVD of centered data as variance-maximizing coordinates — the same $V$ and $\sigma_i$ you already know, reinterpreted for statistics and compression.

If you remember only one sentence from the course: choose the right factorization, and the geometry of the problem becomes visible. $LU$ for systems, $QR$ for least squares, eigenvalues for dynamics, SVD for almost everything else — including the Fourier and PCA tools that close the loop into applications.

Knowledge check#

Browser lab: FFT and PCA#

Run an FFT on a synthetic signal, then PCA on synthetic 2D data via the SVD.

import numpy as np
import matplotlib.pyplot as plt

# Part 1: FFT
t = np.linspace(0, 1, 256, endpoint=False)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
spectrum = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(t), d=t[1] - t[0])

# Part 2: PCA via SVD
rng = np.random.default_rng(0)
mean = [0, 0]
cov = [[3, 1.5], [1.5, 1]]
data = rng.multivariate_normal(mean, cov, size=200)
data_centered = data - data.mean(axis=0)
U, S, Vt = np.linalg.svd(data_centered, full_matrices=False)
explained_variance_ratio = (S ** 2) / np.sum(S ** 2)

print("explained variance ratio =", explained_variance_ratio)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
axes[0].plot(freqs[freqs >= 0], np.abs(spectrum[freqs >= 0]))
axes[0].set_title("FFT magnitude spectrum")
axes[0].set_xlabel("frequency (Hz)")

axes[1].scatter(data_centered[:, 0], data_centered[:, 1], alpha=0.4, s=10)
for i in range(2):
    v = Vt[i] * S[i] / np.sqrt(len(data))
    axes[1].plot([0, v[0]], [0, v[1]], linewidth=2,
                 label=f"PC{i+1} ({explained_variance_ratio[i]*100:.1f}%)")
axes[1].set_aspect("equal")
axes[1].legend()
axes[1].set_title("Principal components")
plt.tight_layout()
plt.show()

print("Notice: the spectrum peaks at 5 Hz and 20 Hz, and PC1 carries most of the explained variance.")

Part 1 builds a sum of 5 Hz and 20 Hz sines and plots the magnitude of np.fft.fft against non-negative frequencies — expect clear peaks at those two bins. Part 2 draws a correlated Gaussian cloud, centers it, runs np.linalg.svd, and overlays the principal directions scaled by singular values. The printed explained_variance_ratio should sum to $1$ and put most of the mass on PC1, matching the elongated shape of the scatter.

Done when#

You can construct the Fourier matrix, explain the FFT's $O(n\log n)$ cost, and perform PCA via the SVD to identify variance-maximizing directions. Check it against the browser lab: the magnitude spectrum peaks at 5 Hz and 20 Hz, the printed explained_variance_ratio sums to 1 with PC1 above 0.9 (roughly [0.95, 0.05]), and the PC1 arrow lies along the long axis of the scatter.

Further Reading#

• Gilbert Strang, Linear Algebra and Its Applications — Chapters on the discrete Fourier transform, FFT, and principal component analysis.
• L. N. Trefethen and D. Bau, Numerical Linear Algebra — SVD, Eckart–Young theorem, and the link to PCA as low-rank approximation.
• S. L. Brunton and J. N. Kutz, Data-Driven Science and Engineering — FFT for signal processing and PCA for dimensionality reduction in modern data science.