Randomized & Approximate Linear Algebra
How PCA/SVD actually scale to millions of rows, via random projection.
On this page
Beginner: exact SVD (section 1.10) of a huge matrix — millions of rows or columns — can simply be too slow to compute, even though you usually only care about the top few singular values and vectors anyway. Randomized numerical linear algebra trades a small, controllable amount of accuracy for an enormous speedup, by first compressing the problem using random projections.
Intermediate: the core enabling fact is the Johnson-Lindenstrauss lemma: projecting high-dimensional points onto a random, much lower-dimensional subspace approximately preserves the distances between them — the number of dimensions you need to keep depends only on how many points you have and how much distortion you'll tolerate, not on the original dimensionality at all.
Advanced: randomized SVD exploits this directly: multiply the original matrix by a small random matrix to project it into a much lower-dimensional sketch, compute an exact (cheap) SVD of that small sketch, then project the result back — producing an excellent approximation to the top singular values/vectors of the original huge matrix, at a tiny fraction of the cost of an exact full SVD.
A is projected through a random matrix Ω into a much smaller matrix Y — an orthonormal basis (via QR, section 1.14) for Y's columns then captures most of A's dominant structure, at a fraction of the cost of factoring A directly.
Suppose A truly has rank k (only k independent directions in its column space) and you sketch it with a random Ω of exactly k columns: Y = AΩ. Write A's compact SVD as A = UΣVᵀ with U, Σ, V all of rank k. Then:
The parenthesized term ΣVᵀΩ is a k×k matrix. Since Ω is random (its columns aren't specially aligned to anything), this k×k matrix is invertible with probability 1 — random matrices are singular only on a measure-zero set of unlucky draws. So Y equals U times an invertible k×k matrix, which means Y's columns span exactly the same k-dimensional space as U's — i.e., exactly A's column space, recovered from a matrix k columns wide instead of A's original width. When A is only approximately low-rank (the realistic case — singular values decay but never hit exactly zero), this argument degrades gracefully rather than breaking outright: the sketch captures the dominant directions well and the small residual is proportional to the singular values being discarded, which is why fast-decaying spectra sketch so well in practice.
Where this is used: this is the exact justification behind sklearn.decomposition.TruncatedSVD and randomized PCA — the random sketch isn't a heuristic approximation of a vague notion, it's provably recovering A's dominant subspace with a concrete, computable failure probability.
This 2D-to-1D toy version is the same idea randomized SVD relies on at massive scale — regenerate to see it hold up across different random directions.
On a matrix this size the speedup is already noticeable; on real production-scale data (millions of rows) it's the difference between feasible and not.
import numpy as np
import time
A = np.random.rand(3000, 500)
t0 = time.time()
U, S, Vt = np.linalg.svd(A, full_matrices=False)
exact_time = time.time() - t0
def randomized_svd(A, k, oversample=10):
n = A.shape[1]
omega = np.random.randn(n, k + oversample)
Y = A @ omega
Q, _ = np.linalg.qr(Y) # section 1.14
B = Q.T @ A
Ub, Sb, Vtb = np.linalg.svd(B, full_matrices=False)
return Q @ Ub[:, :k], Sb[:k], Vtb[:k, :]
t0 = time.time()
Uk, Sk, Vtk = randomized_svd(A, k=10)
rand_time = time.time() - t0
print(f"exact: {exact_time:.3f}s, randomized (top 10): {rand_time:.3f}s")
print("top singular values close?", np.allclose(S[:10], Sk, atol=1e-1))- scikit-learn's
TruncatedSVDandPCA(svd_solver='randomized')use exactly this technique by default once a dataset gets large. - Large-scale recommender systems factor enormous, sparse user-item matrices (section 1.22) using randomized or streaming SVD variants — an exact SVD would be computationally infeasible at that scale.
- Random projection is also used directly as a fast, simple dimensionality reduction technique in its own right, without any SVD step at all, when approximate distance preservation is all that's needed.
- Using too few random projection dimensions relative to how many top singular values you actually need — the "oversampling" parameter in the code above exists specifically to improve accuracy, and skipping it is a common source of poor approximations.
- Assuming randomized SVD approximates the smallest singular values well — it's specifically designed to accurately capture the largest ones; the tail is intentionally sacrificed for speed.
Going deeper
The accuracy of randomized SVD depends on how quickly a matrix's singular values decay — a matrix whose singular values drop off fast (most real-world data matrices) is approximated extremely well; a matrix with a long flat tail of similarly-sized singular values is a much harder case for randomized methods.
At the master level: this entire family of techniques belongs to the broader field of matrix sketching, which also includes methods like CUR decomposition (selecting actual rows and columns of the original matrix rather than random linear combinations) — an active, ongoing area of numerical linear algebra research directly driven by the scale of modern ML data.