PCA as SVD of centered data, and randomized SVD

Before this lesson: lesson 08 (the SVD, Eckart–Young–Mirsky, polar decomposition and the pseudoinverse) and lesson 06 (the spectral theorem and the Rayleigh quotient). Time: about 70 minutes, plus the notebook.

---

Why this matters #

PCA is the first thing to try on any high-dimensional data set: it finds the few directions in which the data varies most. In climate science the same computation is called EOF analysis; in fluid dynamics it is called POD (proper orthogonal decomposition). It is also an interview classic, because it connects statistics (variance), optimization (a constrained maximum) and linear algebra (eigenvectors and singular vectors).

When the data matrix is huge (millions of grid points, many snapshots), an exact SVD is too expensive. The randomized SVD finds the top singular directions with a few passes over the matrix. Both ideas rest on the SVD of lesson 08.

1. Predict first #

2. The picture #

Start with a cloud of points in the plane, centered at the origin and shaped like a long thin ellipse. Pick a unit vector $w$ and project every point onto the line through the origin in direction $w$. The projections are numbers ("shadows" on the line).

• If $w$ points along the long axis, the shadows are spread out: large variance.
• If $w$ points along the short axis, the shadows are bunched together: small variance.

PC1 (the first principal direction) is the $w$ with the largest spread of shadows. There is a second, equivalent view. For each point, its squared length splits into the squared shadow plus the squared perpendicular distance to the line (Pythagoras). The total squared length is fixed. So maximizing the shadows = minimizing the perpendicular distances. PCA finds the best-fitting line (or $k$-dimensional subspace) in this sense.

Interactive figure: loading…

The randomized SVD rests on one picture too. With $A = \sum_j \sigma_j u_j v_j^\top$, a random vector $\omega$ gives

The numbers $v_j^\top \omega$ are random and of size about 1, so the large $\sigma_j$ dominate: a random sample of the range of $A$ points mostly along the top singular directions. A few such samples span (almost) the top-$k$ subspace. "Mostly" is the weak point: if $v_1^\top\omega$ happens to be small, the sample misses $u_1$. Two cheap fixes help. Oversampling takes a few extra random vectors, so one bad draw does not matter. Power iterations multiply by $AA^\top$ again, which turns the weights $\sigma_j$ into $\sigma_j^3, \sigma_j^5, \dots$, so the top direction wins by much more.

Interactive figure: loading…

3. Definitions and main results #

Setup. $X \in \mathbb{R}^{N \times d}$ holds $N$ samples as rows (the convention of scikit-learn), with $d$ features. The mean is $\mu = \frac{1}{N}X^\top \mathbf{1}$. Centering subtracts it from every row: $X_c = X - \mathbf{1}\mu^\top$. The sample covariance matrix is

$S$ is symmetric PSD, so $S = Q\Lambda Q^\top$ with $\lambda_1 \ge \dots \ge \lambda_d \ge 0$ (lesson 06).

Why $N - 1$ and not $N$? The deviations are measured from the sample mean $\mu$, which is fitted to the same data, so they are a little too small on average. Dividing by $N - 1$ instead of $N$ corrects this exactly: $S$ is then an unbiased estimate of the population covariance (Bessel's correction). The principal directions are the same with either choice; only the variances change, by the factor $N/(N-1)$.

Theorem 1 (PCA as an eigenproblem) #

The unit vector $w$ that maximizes the sample variance of the projections $X_c w$ is $w = q_1$, and the maximum variance is $\lambda_1$. Under the extra constraint $w \perp q_1, \dots, q_{k-1}$, the maximizer is $q_k$, with variance $\lambda_k$.

Proof.

1. The projections are $z_i = w^\top \tilde x_i$, where $\tilde x_i$ is row $i$ of $X_c$. Their mean is $w^\top \frac1N \sum_i \tilde x_i = 0$, because the data are centered.
2. So their sample variance is $\frac{1}{N-1}\sum_i (w^\top \tilde x_i)^2 = \frac{1}{N-1}\|X_c w\|^2 = w^\top S w$.
3. Write $w = Qc$. Then $\|c\| = \|w\| = 1$ and $w^\top S w = \sum_j \lambda_j c_j^2 \le \lambda_1 \sum_j c_j^2 = \lambda_1$, with equality for $c = e_1$, that is $w = q_1$.
4. The constraint $w \perp q_1, \dots, q_{k-1}$ means $c_1 = \dots = c_{k-1} = 0$. The same bound gives $w^\top S w \le \lambda_k$, with equality at $w = q_k$. $\blacksquare$

The Lagrange route gives the same result: the gradient of $w^\top S w - \lambda(w^\top w - 1)$ is $2Sw - 2\lambda w$, so a maximizer satisfies $Sw = \lambda w$, and then the objective equals $\lambda$. If $\lambda_1 = \lambda_2$, PC1 is not unique: every unit vector in that eigenspace is a maximizer.

Theorem 2 (PCA from the SVD of centered data) #

Let $X_c = U\Sigma V^\top$ (full SVD: $V$ is $d\times d$, $\Sigma$ is $N\times d$). Then:

• the principal directions are the columns of $V$ (the right singular vectors);
• the variance along direction $i$ is $\lambda_i = \sigma_i^2/(N-1)$ (zero beyond the rank);
• the scores (coordinates of the samples on the principal directions) are $X_c V = U\Sigma$.

Proof.

1. $X_c^\top X_c = V\Sigma^\top U^\top U \Sigma V^\top = V (\Sigma^\top\Sigma) V^\top$, because $U^\top U = I$; $\Sigma^\top\Sigma = \operatorname{diag}(\sigma_1^2, \sigma_2^2, \dots)$ is $d\times d$.
2. Divide by $N-1$: $S = V \,\frac{\Sigma^\top\Sigma}{N-1}\, V^\top$. This is a spectral decomposition of $S$ (orthogonal $V$, diagonal non-negative middle factor).
3. The eigenvalues of $S$ are unique, so $\lambda_i = \sigma_i^2/(N-1)$, and the columns of $V$ are eigenvectors of $S$.
4. $X_c V = U\Sigma V^\top V = U\Sigma$. In practice keep only the top $k$ columns (thin SVD). $\blacksquare$

The explained variance ratio of PC $i$ is $\lambda_i / \sum_j \lambda_j = \sigma_i^2/\|X_c\|_F^2$. If the samples are columns instead ($A = X_c^\top$, size $d \times N$, common in science), the principal directions are the left singular vectors of $A$, and $S = AA^\top/(N-1)$.

At most $N - 1$ components have nonzero variance #

This surprises many people when $d$ is much larger than $N$. Here is why, step by step.

1. Centering makes the rows add up to zero. Each column of $X_c$ has mean 0, so $\mathbf{1}^\top X_c = 0$: the sum of the $N$ rows is the zero vector.
2. So the rows are dependent. The last row equals minus the sum of the others. The row space is spanned by $N - 1$ rows, so $\operatorname{rank}(X_c) \le N - 1$ (and of course $\le d$).
3. Variance needs the row space. The variance along $w$ is $\|X_c w\|^2/(N-1)$. If $w$ is orthogonal to every row, then $X_c w = 0$ and the variance is 0.
4. Count. By Theorem 2, the number of nonzero variances is the number of nonzero $\sigma_i$, which is $\operatorname{rank}(X_c) \le N - 1$. (Equivalently, $\operatorname{rank}(S) = \operatorname{rank}(X_c^\top X_c) = \operatorname{rank}(X_c)$, lesson 02.)

Small cases make it concrete. Two points, after centering, are $\tilde x$ and $-\tilde x$: they lie on one line through 0, so only one component has variance. Three points, after centering, lie in a plane through 0 (the plane of their triangle, moved to the origin), so at most two components have variance, even in $\mathbb{R}^{1000}$. A numerical check: a centered random $50 \times 1000$ matrix has rank 49.

Theorem 3 (maximum variance = minimum reconstruction error) #

For $W \in \mathbb{R}^{d\times k}$ with orthonormal columns, the projection of $\tilde x_i$ onto the span of $W$ is $WW^\top \tilde x_i$, and

Proof.

1. $WW^\top$ is an orthogonal projector, so $\tilde x_i - WW^\top \tilde x_i \perp WW^\top \tilde x_i$.
2. Pythagoras, with $\|WW^\top \tilde x_i\| = \|W^\top \tilde x_i\|$ (orthonormal columns): $\|\tilde x_i\|^2 = \|W^\top \tilde x_i\|^2 + \|\tilde x_i - WW^\top \tilde x_i\|^2$.
3. Sum over $i$: $\sum_i \|W^\top \tilde x_i\|^2 = \|X_c W\|_F^2 = (N-1)\operatorname{tr}(W^\top S W)$. $\blacksquare$

So the subspace that captures the most variance is the one with the smallest reconstruction error. By Eckart–Young (lesson 08), the best choice is $W = V_k$ (top $k$ right singular vectors), with error $\sum_{j>k}\sigma_j^2 = (N-1)\sum_{j>k}\lambda_j$. Eckart–Young compares with every matrix of rank at most $k$, so no other $k$-dimensional subspace does better.

Why center? #

Without centering you analyze the second-moment matrix. Expanding $x_i = \tilde x_i + \mu$ and using $\sum_i \tilde x_i = 0$:

The extra rank-one term pulls the top eigenvector toward the mean $\mu$. PC1 then describes where the data are, not how they spread.

Interactive figure: loading…

Units matter too. PCA is not scale invariant. If one feature is a pressure in pascals (spread about $10^3$) and another a temperature in kelvin (spread about 5), the pressure variance is about $4 \times 10^4$ times larger, and PC1 is almost the pressure axis only because of the units. Then standardize first: divide each centered feature by its standard deviation. This is PCA of the correlation matrix.

Covariance or SVD? #

Forming $X_c^\top X_c$ squares the singular values, so it squares the condition number. In float64 (machine epsilon $\varepsilon = 2^{-52} \approx 2.2 \times 10^{-16}$), directions with $\sigma_j/\sigma_1$ below about $\sqrt{\varepsilon} \approx 1.5 \times 10^{-8}$ are lost in rounding. The SVD works on $X_c$ directly and keeps them. Cost decides the rest:

• $N \gg d$, small $d$: the $d \times d$ covariance is cheap, and the lost directions are usually noise. scikit-learn (checked with version 1.9) uses its covariance_eigh solver by default when $d \le 1000$ and $N \ge 10d$.
• $d \gg N$ ($10^6$ grid points, $N = 50$ members): never form the $d\times d$ matrix. Use the thin SVD of $X_c$ ($O(N^2 d)$ work), or the $N\times N$ Gram matrix $X_c X_c^\top = U\Sigma^2U^\top$ and then $V = X_c^\top U\Sigma^{-1}$ (the "method of snapshots" of POD).
• Both large: randomized SVD.

The randomized SVD (Halko, Martinsson, Tropp 2011; "HMT") #

Input: $A \in \mathbb{R}^{m\times n}$, target rank $k$, oversampling $p$ (5 or 10 is typical), power iterations $q$ (1 or 2 is typical). Let $\ell = k + p$.

1. Draw a Gaussian test matrix $\Omega \in \mathbb{R}^{n\times \ell}$ (independent standard normal entries).
2. $Y = A\Omega$ ($m\times\ell$); orthonormalize, $Q = \operatorname{orth}(Y)$ (thin QR). This is the range finder: $Q$ is an orthonormal basis for (almost) the top-$\ell$ part of the range of $A$.
3. Repeat $q$ times: $\tilde Q = \operatorname{orth}(A^\top Q)$ ($n\times\ell$), then $Q = \operatorname{orth}(A\tilde Q)$ ($m\times\ell$). In exact arithmetic the final $Q$ spans the range of $(AA^\top)^q A\Omega$.
4. $B = Q^\top A$ ($\ell\times n$). Then $A \approx QB = QQ^\top A$.
5. Exact SVD of the small matrix: $B = \tilde U \Sigma V^\top$ ($\tilde U$ is $\ell\times\ell$).
6. $U = Q\tilde U$ ($m\times\ell$). Keep the top $k$ triplets.

The shapes, for $A$ of size $m \times n$:

Matrix	Shape	Made from
$\Omega$	$n \times \ell$	random numbers
$Y = A\Omega$, $Q$	$m \times \ell$	one pass over $A$
$A^\top Q$, $\tilde Q$ (each power iteration)	$n \times \ell$	one pass over $A$
$A\tilde Q$, new $Q$ (each power iteration)	$m \times \ell$	one pass over $A$
$B = Q^\top A$	$\ell \times n$	one pass over $A$
$\tilde U$, $\Sigma$	$\ell \times \ell$	small SVD, no pass
$V^\top$	$\ell \times n$	small SVD, no pass
$U = Q\tilde U$	$m \times \ell$, keep $m \times k$	no pass

So the algorithm touches $A$ in $2q + 2$ passes. The cost is $O(mn\ell)$ per pass plus $O((m+n)\ell^2)$ for the QR and the small SVD, instead of $O(mn\min(m,n))$ for a full SVD.

Why oversampling. The error depends on how well the samples "see" each of the top $k$ directions. Split every random vector $\omega$ into its parts along $v_1, \dots, v_k$ and along the other right singular vectors. The top parts form a $k \times \ell$ matrix $\Omega_1 = V_k^\top\Omega$. A Gaussian matrix stays Gaussian under rotation, so $\Omega_1$ is again a Gaussian matrix.

• With $p = 0$ it is square ($k \times k$), and a random square matrix is sometimes close to singular. Then some top direction gets almost no weight in every sample: the figure in section 2 shows this with one bad $\omega$.
• With a few extra columns ($p = 5$ or $10$), $\Omega_1$ is wide, and a wide Gaussian matrix is very rarely close to rank-deficient. The error is then close to the best possible one.

The error bounds of HMT

Write $V_\perp$ for the other right singular vectors, $\Omega_2 = V_\perp^\top\Omega$, and let $\Sigma_2$ hold $\sigma_{k+1}, \sigma_{k+2}, \dots$. If $\Omega_1$ has full row rank, the deterministic bound of HMT (Theorem 9.1) is

$$\|A - QQ^\top A\|^2 \le \|\Sigma_2\|^2 + \|\Sigma_2\Omega_2\Omega_1^{+}\|^2 \quad(\text{spectral or Frobenius norm}).$$

So the error is small when $\Omega_1$ is well conditioned. For a $k\times(k+p)$ Gaussian matrix, $\mathbb{E}\|\Omega_1^{+}\|_F^2 = k/(p-1)$ for $p \ge 2$. This gives (HMT Theorem 10.5, $k, p \ge 2$, $k + p \le \min(m, n)$)

$$\mathbb{E}\,\|A - QQ^\top A\|_F \le \Big(1 + \frac{k}{p-1}\Big)^{1/2}\Big(\sum_{j>k}\sigma_j^2\Big)^{1/2}.$$

The right factor is the optimal rank-$k$ error (Eckart–Young). The error is also small with high probability: a spectral-norm tail bound (HMT Corollary 10.9, $p \ge 4$) fails with probability at most $6p^{-p}$ (about $0.002$ for $p = 5$).

For power iterations, $(AA^\top)^q A = U\Sigma^{2q+1}V^\top$. In the spectral-norm bound (HMT Corollary 10.10) the factor in front of $\sigma_{k+1}$ is raised to the power $1/(2q+1)$, so it goes to 1 quickly as $q$ grows.

Why power iterations. $(AA^\top)^q A = U\Sigma^{2q+1}V^\top$: same singular vectors, singular values $\sigma_j^{2q+1}$. If the spectrum decays slowly (noisy data), the tail $\sigma_{k+1}, \sigma_{k+2}, \dots$ is large in total and pollutes the samples. Raising to the power $2q+1$ makes it decay fast. For example, $\sigma_j = 1/\sqrt j$ becomes $j^{-3/2}$ with $q = 1$.

Interactive figure: loading…

Why re-orthonormalize. Without the QR in step 3, all columns of $(AA^\top)^qA\Omega$ turn toward $u_1$, as in power iteration. Information about $\sigma_j$ below about $\varepsilon^{1/(2q+1)}\sigma_1$ is lost to rounding (HMT Remark 4.3). Orthonormalizing after every product (subspace iteration) keeps the basis well conditioned.

Answers to section 1. (1) At most 49: $\mathbf{1}^\top X_c = 0$, so the 50 rows of $X_c$ are dependent and $\operatorname{rank}(X_c) \le N - 1$. (2) With exactly $k$ samples the $k\times k$ block $\Omega_1$ can be nearly singular, so one top direction can be badly missed. A few extra columns make that very unlikely.

4. Worked example #

PCA by hand in five steps. Five points in $\mathbb{R}^2$ (rows), so $N = 5$ and $d = 2$:

Step 1. Center. The column means are $\mu = (15/5,\ 10/5) = (3, 2)$. Subtract $\mu$ from every row:

Check: each column of $X_c$ adds up to 0.

Step 2. Covariance (or SVD). Multiply column by column: $(X_c^\top X_c)_{11} = 4 + 1 + 0 + 1 + 4 = 10$, $(X_c^\top X_c)_{12} = 4 + 1 + 0 + 1 + 0 = 6$, $(X_c^\top X_c)_{22} = 4 + 1 + 4 + 1 + 0 = 10$. Divide by $N - 1 = 4$:

Step 3. Directions. A matrix $\begin{pmatrix} a & b \\ b & a \end{pmatrix}$ has eigenvalues $a \pm b$ with eigenvectors $(1, 1)/\sqrt2$ and $(1, -1)/\sqrt2$. So $\lambda_1 = 4$ with $q_1 = (1, 1)/\sqrt2$, and $\lambda_2 = 1$ with $q_2 = (1, -1)/\sqrt2$. The SVD route gives the same: the singular values of $X_c$ are $\sqrt{16} = 4$ and $\sqrt4 = 2$, and $4^2/4 = 4$, $2^2/4 = 1$. ✓

Step 4. Scores. Project the centered rows on the directions:

• on PC1: $X_c q_1 = (-4, -2, 2, 2, 2)/\sqrt2 = (-2\sqrt2, -\sqrt2, \sqrt2, \sqrt2, \sqrt2)$;
• on PC2: $X_c q_2 = (0, 0, -2, 0, 2)/\sqrt2 = (0, 0, -\sqrt2, 0, \sqrt2)$.

Check: the PC1 scores have mean 0 and variance $(8 + 2 + 2 + 2 + 2)/4 = 4 = \lambda_1$. ✓

Step 5. Explained variance. PC1 explains $\lambda_1/(\lambda_1 + \lambda_2) = 4/5 = 80\%$ of the variance, PC2 the other $20\%$. If you keep only PC1, the reconstruction error is the sum of the squared PC2 scores: $2 + 2 = 4 = \sigma_2^2 = (N - 1)\lambda_2$, as Theorem 3 says. ✓

Without centering: $X^\top X = \begin{pmatrix} 55 & 36 \\ 36 & 30\end{pmatrix} = X_c^\top X_c + 5\mu\mu^\top$. Its top eigenvector is about $(0.815, 0.580)$: only $1.7°$ from the mean direction $(3, 2)/\sqrt{13}$, and $9.6°$ from the true PC1.

Randomized SVD by hand. $A = \operatorname{diag}(3, 1)$, $k = 1$, $p = 0$, and the "random" vector $\Omega = (1, 1)^\top$. The optimal rank-1 error is $\sigma_2 = 1$.

$q$	$Y = (AA^\top)^q A\Omega$	estimate $\lVert Q^\top A\rVert$ of $\sigma_1 = 3$	error $\lVert A - QQ^\top A\rVert_2$
0	$(3, 1)$	$\sqrt{8.2} \approx 2.864$	$\sqrt{1.8} \approx 1.342$
1	$(27, 1)$	$\approx 2.998$	$\approx 1.005$

The weight ratio of $u_1$ to $u_2$ in $Y$ grows from $3:1$ to $27:1 = 3^3:1$. That is the $\sigma^{2q+1}$ effect.

5. In code #

Libraries check the identities of section 3. The notebook asks you to write pca and randomized_svd yourself.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.utils.extmath import randomized_svd

# 1) The worked example: PCA variances are sigma_i^2 / (N - 1).
X = np.array([[1, 0], [2, 1], [3, 4], [4, 3], [5, 2]], dtype=float)
pca = PCA(n_components=2).fit(X)
s = np.linalg.svd(X - X.mean(axis=0), compute_uv=False)
print(s)                             # [4. 2.]
print(pca.explained_variance_)       # [4. 1.]
print(s**2 / (len(X) - 1))           # [4. 1.]
print(pca.explained_variance_ratio_) # [0.8 0.2]
print(np.abs(pca.components_))       # rows: (1,1)/√2 and (1,-1)/√2, up to sign

# 2) Forming Z^T Z squares the condition number: a singular value 1e-9 is lost.
rng = np.random.default_rng(0)
P, _ = np.linalg.qr(rng.standard_normal((1000, 2)))
R = np.array([[1.0, -1.0], [1.0, 1.0]]) / np.sqrt(2)   # rotate the input axes by 45°
Z = P @ np.diag([1.0, 1e-9]) @ R.T                     # singular values 1 and 1e-9
print(np.linalg.svd(Z, compute_uv=False))      # [1.e+00 1.e-09]: correct
print(np.linalg.eigvalsh(Z.T @ Z))             # small one: rounding noise (about -3e-17), not 1e-18

# 3) Randomized SVD on a slowly decaying spectrum sigma_j = 1/sqrt(j).
m = n = 500
U, _ = np.linalg.qr(rng.standard_normal((m, m)))
V, _ = np.linalg.qr(rng.standard_normal((n, n)))
sig = 1 / np.sqrt(np.arange(1, n + 1))
A = (U * sig) @ V.T
k = 10
for q in [0, 1, 2]:
    Uk, sk, Vtk = randomized_svd(A, k, n_oversamples=10, n_iter=q,
                                 power_iteration_normalizer="QR", random_state=0)
    err = np.linalg.norm(A - (Uk * sk) @ Vtk, 2)
    print(q, round(err / sig[k], 3))  # ratio to the optimum sigma_{k+1}: about 1.9, 1.02, 1.00

In part 2 the small eigenvalue comes out negative: pure rounding noise. In part 3, without power iterations the error is almost twice the optimum; two power iterations reach it. (For this matrix, raising $p$ from 10 to 40 with $q = 0$ helps less than one power iteration.)

6. Where it shows up #

• Climate and weather (EOF analysis). Rows are times or ensemble members, columns are grid points. The principal directions are spatial patterns (EOFs), the scores are their time series (PCs). On a latitude–longitude grid, multiply each column by $\sqrt{\cos(\text{latitude})}$ first, so that each grid point counts with its area.
• Ensemble forecasts. The leading PCs of the ensemble anomalies give a low-rank model of the forecast uncertainty. Without localization, the ensemble Kalman filter update lies in the span of these anomalies.
• Reduced-order models (POD). The leading modes of simulation snapshots give a small basis for a fast surrogate model.
• ML. Embedding visualization, whitening, feature compression. A linear autoencoder trained with squared error learns the PCA subspace (not necessarily the individual directions).
• Libraries. sklearn.decomposition.PCA(svd_solver="randomized"), torch.svd_lowrank and torch.pca_lowrank (based on HMT Algorithm 5.1). For small dense matrices a full SVD is faster.

7. Common mistakes #

• Not centering. PC1 then points toward the mean. (torch.svd_lowrank does not center; torch.pca_lowrank does by default.)
• Mixing conventions. Samples as rows: directions = right singular vectors $V$. Samples as columns: directions = left singular vectors $U$.
• $1/N$ vs $1/(N-1)$. The directions are the same; the variances differ by the factor $N/(N-1)$. scikit-learn uses $1/(N-1)$.
• Comparing signs. Each principal direction is defined only up to sign.
• Not scaling features with different units. PCA is not scale invariant. If one variable is in pascals and another in kelvin, standardize first (PCA of the correlation matrix).
• Expecting more than $N - 1$ components. With $N$ samples, at most $N - 1$ components have nonzero variance, however large $d$ is.
• Randomized SVD with $p = 0$, or power iterations without re-orthonormalization.

8. Check yourself #

1. $X_c$ ($N = 10$) has singular values $6, 3, 0$. Give the variances along the PCs and the explained variance ratio of PC1.
   Answer

   $\lambda = 36/9 = 4$, $9/9 = 1$, $0$. Ratio $4/5 = 80\%$.
2. Data with mean $\mu = (10, 0)$ and a small spread. What does "PCA" without centering give as PC1?
   Answer

   About $e_1$, the mean direction: $X^\top X/(N-1) = S + \frac{N}{N-1}\mu\mu^\top$ is dominated by $\mu\mu^\top$ (largest eigenvalue about $100N/(N-1)$).
3. Randomized SVD of a $10^6 \times 10^4$ matrix with $k = 20$, $p = 10$, $q = 2$: shapes of $\Omega$, $Q$, $B$, and the number of passes over $A$?
   Answer

   $\Omega$: $10^4\times 30$; $Q$: $10^6 \times 30$; $B = Q^\top A$: $30\times 10^4$. Passes: $2q+2 = 6$.
4. The singular values decay like $\sigma_j = 1/\sqrt j$ and the error is twice the optimum. Increase $p$ or $q$? Why?
   Answer

   Mainly $q$. The tail is flat, so many tail directions leak into the samples. Power iterations turn $\sigma_j$ into $\sigma_j^{2q+1}$ (here $j^{-5/2}$ for $q = 2$), which suppresses the tail. Oversampling alone helps more slowly.
5. With $q = 3$ power iterations and no re-orthonormalization in float64, roughly which singular values are lost?
   Answer

   Those below about $\varepsilon^{1/(2q+1)}\sigma_1 = (2^{-52})^{1/7}\sigma_1 \approx 6\times10^{-3}\,\sigma_1$.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What does the first principal direction maximize, and which eigenproblem gives it?

   Show answer

2. How do you get the principal directions and variances from the SVD of the centered data $X_c$ (samples as rows)?

   Show answer

3. Why is "maximum variance" the same as "minimum reconstruction error" in PCA?

   Show answer

4. Why must you center the data before PCA?

   Show answer

5. $N$ samples in $\mathbb{R}^d$ with $d > N$. At most how many principal components have nonzero variance, and why?

   Show answer

6. Why does the sample covariance divide by $N - 1$ and not by $N$?

   Show answer

7. Randomized SVD: what are the first two steps (the range finder)?

   Show answer

8. Randomized SVD: what happens after the range finder?

   Show answer

9. Why does the randomized SVD oversample ($\ell = k + p$ with $p \approx 5$–$10$)?

   Show answer

10. Why do power iterations help the randomized SVD, and what must be done after each product?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Derive PCA." — Variance of projections $= w^\top S w$; maximize on the unit sphere; get $Sw = \lambda w$; then link to the SVD: $\lambda_i = \sigma_i^2/(N-1)$.
• "Eigendecomposition of the covariance or SVD of the data?" — SVD for accuracy (no squared condition number) and when $d$ is large; covariance when $N \gg d$ and $d$ is small.
• "How do you run PCA on a $10^7 \times 10^5$ matrix?" — Randomized SVD (or a Krylov method), with $k + p$ samples and 1–2 power iterations: $2q+2$ passes over the data.
• "Why does a random projection find the top singular subspace?" — $A\omega$ weights $u_j$ by $\sigma_j$; oversampling makes $\Omega_1$ well conditioned; the error is close to $\sigma_{k+1}$.

10. Practice #

The notebook for this lesson runs in Colab. It has four exercises, each with a check cell and a folded solution:

1. pca: principal directions, explained variances and scores from the SVD of the centered data; compare with scikit-learn on the digits data set.
2. reconstruction_error: project the digits onto the top $k$ components, check Theorem 3 (error $= \sum_{j>k}\sigma_j^2$), and find how many components keep 90% of the variance.
3. uncentered_pc1: see how far the top direction of $X^\top X$ moves when you shift the data.
4. randomized_svd: the HMT algorithm with oversampling, power iterations and QR after every product; compare its error with $\sigma_{k+1}$ and with a full SVD on a $1000 \times 400$ matrix.

11. Watch and read #

Principal Component Analysis (PCA)Steve Brunton · Data-Driven Science and Engineering · 14 min

StatQuest: Principal Component Analysis (PCA), Step-by-StepJosh Starmer · StatQuest · 22 min

Randomized Singular Value Decomposition (SVD)Steve Brunton · Data-Driven Science and Engineering · 13 min

7. Eckart-Young: The Closest Rank k Matrix to AGilbert Strang · MIT 18.065, lecture 7 (includes PCA) · 47 min

• Halko, Martinsson, Tropp, Finding structure with randomness (SIAM Review 2011): §1.6 (prototype randomized SVD), §4.1–4.5 (Algorithms 4.1, 4.3, 4.4 and Remark 4.3), §9.2 (Theorem 9.1), §10.2–10.4 (Theorem 10.5, Corollaries 10.9 and 10.10).
• Brunton and Kutz, Data-Driven Science and Engineering, 2nd edition (free PDF), chapter 1: §1.5 (PCA), §1.6 (eigenfaces), §1.8 (randomized SVD).
• scikit-learn user guide, Principal component analysis: the solvers (full, covariance_eigh, randomized) and when each is used.