Sketching, Johnson–Lindenstrauss and random matrix theory

Before this lesson: lesson 09 (PCA and randomized SVD), lesson 12 (QR and least squares) and lesson 17 (perturbation theory: Weyl and Davis–Kahan). Time: about 90 minutes, plus the notebook.

---

Why this matters #

Some matrices are too big to factor, or even to read many times. A sketch multiplies the data by a short random matrix and works with the small result. The randomized SVD of lesson 09 was the first example. The same idea solves huge least-squares problems and compresses embeddings and gradients.

Randomness also explains noise. With $n = 50$ ensemble members and $p = 10^6$ grid points, the sample covariance of an ensemble forecast is mostly sampling noise. Random matrix theory says exactly how pure noise looks, so you can tell which eigenvalues or singular values carry signal. This replaces "look for the elbow" in a PCA or EOF analysis by a rule with a reason.

1. Predict first #

2. The picture #

One vector first. Fix $x \in \mathbb{R}^D$ and draw $g \in \mathbb{R}^D$ with independent $N(0,1)$ entries. Then $g^\top x \sim N(0, \|x\|^2)$, so $(g^\top x)^2$ is an unbiased estimate of $\|x\|^2$. It is a bad estimate: its relative standard deviation is $\sqrt 2$. Now use $k$ independent rows, stacked as $G \in \mathbb{R}^{k\times D}$. The average $\|Gx\|^2/k$ is still unbiased, and its relative standard deviation is $\sqrt{2/k}$: 0.1 for $k = 200$.

So a random map $x \mapsto Gx/\sqrt k$ keeps one length with error about $\sqrt{2/k}$, whatever $D$ is. The tails are close to Gaussian, so the probability of a large error falls exponentially in $k$. This is why we can keep $n^2$ distances at once with only $k \propto \log n$: a union bound over $n^2$ events multiplies the failure probability by $n^2$, and in the exponent this is only an extra term $2\ln n$.

Interactive figure: loading…

A whole subspace. A $d$-dimensional subspace holds infinitely many vectors. But its unit sphere is covered by a finite net of $e^{O(d)}$ points (every unit vector is close to one of them), so it behaves like $e^{O(d)}$ vectors. Then "$\log(\text{number of points})$" becomes $d$, and a sketch needs $k \propto d$ rows. Section 3 shows this with a figure.

Noise has a shape. Estimate a covariance from $n$ samples in $p$ dimensions. If $p/n$ is not small, the sample eigenvalues spread out, even when all true eigenvalues are equal. The spread is deterministic as $n, p \to \infty$: the Marchenko–Pastur law. A true signal direction is visible only if its eigenvalue sticks out of this "bulk".

Interactive figure: loading…

3. Definitions and main results #

Subspace embeddings #

Start with a picture. Take a plane (a 2-dimensional subspace) inside $\mathbb{R}^{5000}$ and its unit circle: all vectors of the plane with length 1. A sketch $S$ maps each of them to a vector $Sy$. If we draw each $y$ at the length $\|Sy\|$, the circle becomes an ellipse. A good sketch gives an ellipse that is almost round.

Interactive figure: loading…

From the picture to the formula, step by step. Let $U \in \mathbb{R}^{n\times d}$ have orthonormal columns that span the subspace $\mathcal U$.

1. Every $y \in \mathcal U$ is $y = Uz$ for one $z \in \mathbb{R}^d$, and $\|y\| = \|z\|$ because $U^\top U = I$.
2. Then $\|Sy\|^2 = z^\top (U^\top S^\top S U) z = z^\top M z$, with the small $d\times d$ matrix $M = U^\top S^\top SU$.
3. "All lengths are kept within $1 \pm \varepsilon$" means $(1-\varepsilon)\|z\|^2 \le z^\top Mz \le (1+\varepsilon)\|z\|^2$ for every $z$.
4. By the Rayleigh quotient (lesson 06), this holds exactly when every eigenvalue of $M$ lies in $[1-\varepsilon, 1+\varepsilon]$, that is, $\|I_d - M\|_2 \le \varepsilon$.
5. The eigenvalues of $M = (SU)^\top(SU)$ are the squared singular values of $SU$. So every singular value of $SU$ lies in $[\sqrt{1-\varepsilon}, \sqrt{1+\varepsilon}]$: the ellipse lies in the ring.

Definition. Let $\mathcal{U} \subseteq \mathbb{R}^n$ be a subspace and $0 < \varepsilon < 1$. A matrix $S \in \mathbb{R}^{k\times n}$ is an $\varepsilon$-subspace embedding for $\mathcal{U}$ if, for all $y \in \mathcal{U}$,

By steps 1–5, this is the same as $\|I_d - U^\top S^\top S U\|_2 \le \varepsilon$. A random $S$ is oblivious if it works with probability $\ge 1-\delta$ for any fixed $\mathcal U$ of dimension $d$, without looking at the data.

Three standard sketches for a $d$-dimensional subspace (rows up to constant factors, for a fixed failure probability $\delta$; $\operatorname{nnz}(A)$ is the number of nonzero entries of $A$):

Sketch	What $S$ is	Rows $k$ (order)	Cost of $SA$, $A \in \mathbb{R}^{n\times d}$
Gaussian	$G/\sqrt k$, $G_{ij} \sim N(0,1)$	$d/\varepsilon^2$	$O(k\cdot\operatorname{nnz}(A))$, dense: $O(knd)$
SRHT	random signs, then a fast Hadamard transform, then $k$ random rows	$d\log d/\varepsilon^2$	$O(nd\log k)$
CountSketch	each row of $A$ is added, with a random sign, to one random bucket	$d^2/\varepsilon^2$	$O(\operatorname{nnz}(A))$

The Gaussian sketch needs the fewest rows but is the most expensive to apply. CountSketch is the cheapest to apply (one pass over the nonzeros) but needs more rows.

Exact row counts and the SRHT

Woodruff (2014), Theorems 6, 7 and 9, give these sizes for an $\varepsilon$-subspace embedding of a $d$-dimensional subspace with probability at least $1-\delta$:

Sketch	What $S$ is	Rows $k$
Gaussian	$G/\sqrt k$, $G_{ij} \sim N(0,1)$	$O((d + \log(1/\delta))/\varepsilon^2)$
SRHT	$\sqrt{n/k}\,RHD$: random signs $D$, normalized Walsh–Hadamard $H$ (entries $\pm1/\sqrt n$), $k$ sampled rows $R$	$O(\varepsilon^{-2}\log d\,(\sqrt d + \sqrt{\log n})^2)$
CountSketch	one $\pm1$ per column, in a random row	$O(d^2/(\delta\varepsilon^2))$

Why the SRHT needs $D$ and $H$: uniform row sampling alone fails when a few rows carry most of the energy (a vector like $e_1$). The random signs and the Hadamard transform spread the energy of every vector over all coordinates, so that sampling $k$ rows uniformly becomes safe.

Two ways to use a sketch for least squares #

We want $x^* = \arg\min_x \|Ax - b\|$ with $A \in \mathbb{R}^{n\times d}$ and $n \gg d$. There are two ways to use a sketch:

	Sketch-and-solve	Sketch-and-precondition
What it does	solves the small problem $\min_x \|S(Ax - b)\|$	uses $SA = QR$ only to build $R$; then runs LSQR on the full problem with the preconditioner $R^{-1}$
Accuracy	residual within a factor $\approx 1+\varepsilon$ of the optimum; $x$ itself can be far off	full accuracy (machine precision), like a direct solver
Cost	one sketch + one small QR	one sketch + one small QR + about 30 cheap iterations, each one product with $A$ and $A^\top$
When to use	you need a good fit (a small residual), not exact coefficients	you need accurate coefficients $x$, or $A$ is ill-conditioned

The next two theorems explain both rows of "Accuracy".

Theorem 1 (sketch-and-solve least squares) #

Let $A \in \mathbb{R}^{n\times d}$, $b \in \mathbb{R}^n$, $x^* = \arg\min_x\|Ax - b\|$ and $r^* = Ax^* - b$. Let $S$ be an $\varepsilon$-subspace embedding for $\mathcal U = \operatorname{range}([A\;\; b])$ (dimension at most $d+1$), and $\tilde x = \arg\min_x \|S(Ax - b)\|$. Then

Proof.

1. For every $x$, the vector $Ax - b$ lies in $\mathcal U$, so the embedding bounds apply to it.
2. Lower bound at $\tilde x$: $\|A\tilde x - b\|^2 \le \frac{1}{1-\varepsilon}\|S(A\tilde x - b)\|^2$.
3. $\tilde x$ minimizes the sketched problem: $\|S(A\tilde x - b)\|^2 \le \|S(Ax^* - b)\|^2$.
4. Upper bound at $x^*$: $\|S(Ax^* - b)\|^2 \le (1+\varepsilon)\|Ax^* - b\|^2$.
5. Chain 2–4 and take square roots. $\blacksquare$

For small $\varepsilon$ the factor is about $1+\varepsilon$. The guarantee is about the residual, not about $x$. Since $r^* \perp \operatorname{range}(A)$, Pythagoras gives $\|A(\tilde x - x^*)\|^2 = \|A\tilde x - b\|^2 - \|r^*\|^2 \le \frac{2\varepsilon}{1-\varepsilon}\|r^*\|^2$, so

With a large residual and an ill-conditioned $A$ (small $\sigma_{\min}$), $\tilde x$ can be far from $x^*$: section 5 shows a 170% error. For a Gaussian sketch the average is known exactly: $\mathbb{E}\|A\tilde x - b\|^2 = (1 + \frac{d}{k-d-1})\|r^*\|^2$.

Theorem 2 (sketch-and-precondition) #

Let $S$ be an $\varepsilon$-subspace embedding for $\operatorname{range}(A)$ and $SA = QR$ a thin QR ($R \in \mathbb{R}^{d\times d}$). Then $\kappa(AR^{-1}) \le \sqrt{(1+\varepsilon)/(1-\varepsilon)}$, where $\kappa$ is the 2-norm condition number (lesson 13).

In words: $R$ undoes almost all of the bad scaling of $A$. The matrix $AR^{-1}$ has nearly orthonormal columns, whatever $\kappa(A)$ was.

Proof of Theorem 2

1. $SAR^{-1} = Q$ has orthonormal columns, so $\|SAR^{-1}z\| = \|z\|$ for every $z \in \mathbb{R}^d$.
2. $AR^{-1}z \in \operatorname{range}(A)$, so the embedding gives $(1-\varepsilon)\|AR^{-1}z\|^2 \le \|SAR^{-1}z\|^2 = \|z\|^2 \le (1+\varepsilon)\|AR^{-1}z\|^2$.
3. Divide by $\|z\|^2$: $\|AR^{-1}z\|/\|z\|$ lies in $[1/\sqrt{1+\varepsilon},\, 1/\sqrt{1-\varepsilon}]$ for every $z \ne 0$. So every singular value of $AR^{-1}$ lies in this interval, and their ratio is at most $\sqrt{(1+\varepsilon)/(1-\varepsilon)}$. $\blacksquare$

Now run LSQR (CG on the normal equations, done stably) on $\min_y \|AR^{-1}y - b\|$ and set $x = R^{-1}y$. With $\kappa = \kappa(AR^{-1})$, the error $\|A(x^{(m)} - x^*)\|$ after $m$ iterations is at most $2\big((\kappa - 1)/(\kappa + 1)\big)^m$ times its start value (Woodruff §2.6, Eq. (15)). Even a crude sketch with $\varepsilon = 1/2$ gives $\kappa \le \sqrt3$ and a factor $0.27$ per iteration: about 27 iterations to reach $10^{-15}$, whatever $\kappa(A)$ is. This reaches full accuracy because the sketch only builds the preconditioner; the iteration solves the original problem. Blendenpik and LSRN work this way.

Theorem 3 (Johnson–Lindenstrauss) #

Let $x_1, \dots, x_n \in \mathbb{R}^D$, $0 < \varepsilon < 1$, $0 < \delta < 1$. Let $G \in \mathbb{R}^{k\times D}$ have independent $N(0,1)$ entries and $f(x) = Gx/\sqrt k$. If

then with probability at least $1 - \delta$, for all pairs $i, j$:

So $k = O(\varepsilon^{-2}\log(n/\delta))$, independent of $D$.

Proof, step by step. Each step has the numbers for $n = 10^4$ points, $\varepsilon = 0.2$ and $\delta = 0.01$, where the bound asks for $k = 2657$.

1. Reduce to one vector. $f$ is linear, so $f(x_i) - f(x_j) = f(x_i - x_j)$. It suffices to show that for each fixed unit vector $u$ (the normalized difference), $Q = \|Gu\|^2$ satisfies $(1-\varepsilon)k \le Q \le (1+\varepsilon)k$, except with a small probability. Numbers: there are $n(n-1)/2 = 49{,}995{,}000 \approx 5\cdot10^7$ difference vectors.
2. Exact distribution. Row $i$ of $G$ gives $g_i^\top u \sim N(0, \|u\|^2) = N(0,1)$, and the rows are independent. So $Q$ is a sum of $k$ squared standard normals: $Q \sim \chi^2_k$, $\mathbb{E}Q = k$. Numbers: we want $Q$ in $[0.8k, 1.2k] = [2125.6, 3188.4]$.
3. Chernoff, upper tail. For $Z \sim N(0,1)$ and $\lambda < 1/2$, $\mathbb{E}e^{\lambda Z^2} = (1-2\lambda)^{-1/2}$. By Markov's inequality applied to $e^{\lambda Q}$, $P(Q \ge (1+\varepsilon)k) \le e^{-\lambda(1+\varepsilon)k}(1-2\lambda)^{-k/2}$. Choose $\lambda = \frac{\varepsilon}{2(1+\varepsilon)}$, so $1 - 2\lambda = \frac{1}{1+\varepsilon}$. The bound becomes $\big((1+\varepsilon)e^{-\varepsilon}\big)^{k/2} = \exp\!\big(-\tfrac k2(\varepsilon - \ln(1+\varepsilon))\big)$.
4. Use $\ln(1+\varepsilon) \le \varepsilon - \varepsilon^2/2 + \varepsilon^3/3$: $P(Q \ge (1+\varepsilon)k) \le \exp\!\big(-\tfrac k2(\varepsilon^2/2 - \varepsilon^3/3)\big)$. Numbers: $\varepsilon^2/2 - \varepsilon^3/3 = 0.02 - 0.00267 = 0.01733$, and $\exp(-\tfrac{2657}{2}\cdot 0.01733) = e^{-23.03} \approx 1.0\cdot10^{-10}$.
5. Lower tail, the same way with $e^{-\lambda Q}$ and $\lambda = \frac{\varepsilon}{2(1-\varepsilon)}$: $P(Q \le (1-\varepsilon)k) \le \big((1-\varepsilon)e^{\varepsilon}\big)^{k/2} \le e^{-k\varepsilon^2/4}$, using $\ln(1-\varepsilon) \le -\varepsilon - \varepsilon^2/2$. This is no larger than the bound of step 4. Numbers: one pair fails with probability at most $2 \cdot 1.0\cdot10^{-10} = 2.0\cdot10^{-10}$.
6. Union bound. There are fewer than $n^2/2$ pairs, each failing with probability at most $2\exp(-\tfrac k2(\varepsilon^2/2 - \varepsilon^3/3))$. The total is at most $n^2\exp(-\tfrac k2(\varepsilon^2/2 - \varepsilon^3/3))$, which is $\le \delta$ exactly when $k$ satisfies the bound above. Numbers: $10^8 \cdot e^{-23.03} \approx 0.0100 \le \delta$. $\blacksquare$

The bound is pessimistic. For $k = 2657$, the exact $\chi^2_k$ probability of leaving $[0.8k, 1.2k]$ is about $3.4\cdot10^{-12}$, about 60 times smaller than the Chernoff bound of step 5. This is why the histogram in section 2 almost always stays inside the aqua lines.

Remarks.

• The map is oblivious: it does not look at the data, unlike PCA.
• Inner products survive too. By $4\langle u, v\rangle = \|u+v\|^2 - \|u-v\|^2$, if $f$ keeps the squared norms of $u \pm v$, then $|\langle f(u), f(v)\rangle - \langle u, v\rangle| \le \varepsilon$ for unit $u, v$.
• The $\log n$ cannot be improved, even by nonlinear maps (Vershynin, Remark 5.3.4), and Larsen and Nelson (2017) showed that $\varepsilon^{-2}\log n$ as a whole is optimal (for $\varepsilon$ not extremely small).
• Random signs ($\pm1$ entries), sparse JL matrices and the fast JL transform (the SRHT) give the same guarantee, up to constants, at lower cost.
• With an $e^{O(d)}$-point net on the unit sphere of a $d$-dimensional subspace, the same proof gives the Gaussian row of the table: $k = O((d + \log(1/\delta))/\varepsilon^2)$.

The Marchenko–Pastur law #

Setup. $X \in \mathbb{R}^{n\times p}$ has independent entries with mean 0 and variance $\sigma^2$: $n$ samples (rows), $p$ features. The sample covariance is $\hat\Sigma = X^\top X/n$ ($p\times p$), and the aspect ratio is $\gamma = p/n$.

Theorem 4 (Marchenko–Pastur 1967). As $n, p \to \infty$ with $p/n \to \gamma \in (0, \infty)$, the histogram of the $p$ eigenvalues of $\hat\Sigma$ converges to the density

on $[\lambda_-, \lambda_+]$, where $\lambda_\pm = \sigma^2(1 \pm \sqrt\gamma)^2$. For $\gamma \le 1$ this is the whole law. For $\gamma > 1$ this part has mass $1/\gamma$, and a point mass $1 - 1/\gamma$ sits at 0 ($\hat\Sigma$ has rank at most $n < p$).

What to read from it:

• The mean of the law is $\sigma^2$: the trace is right, but the eigenvalues spread around it.
• The spread is set by $\gamma$ alone. For $\gamma \to 0$ ($n \gg p$) the support shrinks to $\{\sigma^2\}$.
• Singular values of $X$ itself fill about $[\sigma(\sqrt n - \sqrt p),\, \sigma(\sqrt n + \sqrt p)]$. Non-asymptotic version for Gaussian entries: $\mathbb{E}\|X\|_2 \le \sigma(\sqrt n + \sqrt p)$ (Vershynin, Theorem 7.3.1).
• With a finite fourth moment, the largest eigenvalue converges to $\lambda_+$. Its fluctuations are of relative order $n^{-2/3}$ (the Tracy–Widom law), so the edge is sharp.
• Link to lesson 17: for $\Sigma = \sigma^2 I$, $\|\hat\Sigma - \Sigma\|_2 \to \sigma^2(2\sqrt\gamma + \gamma)$. Weyl's inequality then lets each sample eigenvalue move by up to this amount.

Theorem 5 (spiked covariance, BBP phase transition). Let the rows have covariance $\sigma^2(I + \theta vv^\top)$, with $\|v\| = 1$ and spike strength $\theta > 0$. As above, with $\gamma \le 1$:

• If $\theta > \sqrt\gamma$: the top sample eigenvalue converges to $\sigma^2(1+\theta)(1+\gamma/\theta) > \lambda_+$, and the top sample eigenvector $\hat v$ satisfies $\langle\hat v, v\rangle^2 \to \dfrac{1 - \gamma/\theta^2}{1 + \gamma/\theta}$.
• If $\theta \le \sqrt\gamma$: the top eigenvalue sticks to $\lambda_+$, and $\langle\hat v, v\rangle \to 0$.

Two lessons: a signal weaker than $\sqrt\gamma$ is undetectable by PCA, and a detected signal is biased: the sample eigenvalue is too large and the eigenvector is tilted (Davis–Kahan in action). In the figure of section 2, set $\gamma = 0.25$ and move $\theta$ across $0.5$.

More on the BBP transition

• Where the threshold comes from. Write $\ell(\theta) = (1+\theta)(1+\gamma/\theta) = 1 + \gamma + \theta + \gamma/\theta$. Its derivative $1 - \gamma/\theta^2$ is zero at $\theta = \sqrt\gamma$, where $\ell(\sqrt\gamma) = 1 + 2\sqrt\gamma + \gamma = (1+\sqrt\gamma)^2 = \lambda_+/\sigma^2$. So the formula for the top eigenvalue touches the bulk edge exactly at the threshold. Below it the formula would rise again (for example $\ell(0.4) = 2.275$ for $\gamma = 0.25$), but it no longer describes an eigenvalue: the top eigenvalue stays at the edge.
• The overlap $(1 - \gamma/\theta^2)/(1 + \gamma/\theta)$ is 0 at $\theta = \sqrt\gamma$ and tends to 1 as $\theta \to \infty$. Near the threshold the eigenvalue is visible but the eigenvector is almost useless.
• Sources. Baik, Ben Arous and Péché (2005) for the eigenvalue (complex Gaussian data), Baik and Silverstein (2006) and Paul (2007) for real data and the eigenvector.
• Finite samples. At finite $n$ the top noise eigenvalue fluctuates around $\lambda_+$ on the Tracy–Widom scale $n^{-2/3}$, and the transition is smooth instead of sharp.

How to decide if a top eigenvalue is signal. A step list for a sample covariance with $n$ samples and $p$ features:

1. Compute $\gamma = p/n$.
2. Estimate the noise variance $\sigma^2$: from data with no signal, or from the bulk (for example the median eigenvalue compared with the median of the MP law).
3. Compute the bulk edge $\lambda_+ = \sigma^2(1+\sqrt\gamma)^2$.
4. An eigenvalue clearly above $\lambda_+$ (with a margin for finite-$n$ fluctuations) is signal. Everything at or below $\lambda_+$ is consistent with pure noise.
5. For a detected eigenvalue $\hat\lambda$, remove the bias. Set $\ell = \hat\lambda/\sigma^2$ and solve $\ell = (1+\theta)(1+\gamma/\theta)$ for $\theta$: $\theta = \tfrac12\big[(\ell - 1 - \gamma) + \sqrt{(\ell - 1 - \gamma)^2 - 4\gamma}\,\big]$. The true variance in that direction is $\sigma^2(1+\theta)$, not $\hat\lambda$.
6. A direction you cannot see may still exist: it is only known that $\theta \le \sqrt\gamma$. More samples (smaller $\gamma$) lower the threshold.

Signal plus noise for one data matrix. Let $Y = X_0 + \sigma Z$ be $m\times n$ with $m \le n$, $\beta = m/n$, $X_0$ of low rank and $Z$ with independent $N(0,1)$ entries. The noise singular values fill $[\sqrt n\,\sigma(1-\sqrt\beta),\; \sqrt n\,\sigma(1+\sqrt\beta)]$. A rule to separate:

1. Estimate $\sigma$ (from a noise-only region, or from the median singular value).
2. Keep the singular values above the bulk edge $(1+\sqrt\beta)\sqrt n\,\sigma$ (plus a margin): below it, signal and noise cannot be told apart.
3. If the goal is the smallest mean squared error of $\hat X_0$, use the slightly higher Gavish–Donoho threshold $\tau^* = \lambda^*(\beta)\sqrt n\,\sigma$, with $\lambda^*(1) = 4/\sqrt3 \approx 2.309$ for square matrices (2.858 times the median singular value when $\sigma$ is unknown). Components just above the edge are real but so noisy that keeping them increases the error.

For a square matrix and a rank-one signal of size $x\sqrt n\,\sigma$, the observed top singular value is about $\sqrt n\,\sigma\,(1+x^2)/x$ if $x > 1$, and stays at the edge $2\sqrt n\,\sigma$ if $x \le 1$.

High-dimensional geometry and superposition #

• Norms concentrate. For $g \sim N(0, I_D)$, $\|g\| = \sqrt D + O(1)$ (Vershynin §3.1): random vectors live near a sphere.
• Random directions are almost orthogonal. For independent uniform unit vectors $u, v \in \mathbb{R}^d$, $\langle u, v\rangle$ has mean 0 and variance $1/d$, and $P(|\langle u, v\rangle| \ge t) \le 2e^{-dt^2/2}$. The union bound of the JL proof then gives $N \approx e^{d\varepsilon^2/4}$ directions with all $\lvert\cos\rvert \le \varepsilon$; equivalently, $N$ random directions have $\max\lvert\cos\rvert \lesssim \sqrt{4\ln N/d}$. In $d = 1000$ you get about 20,000 directions with $\lvert\cos\rvert \le 0.2$, but only 1000 exactly orthogonal ones.
• Superposition (Elhage et al. 2022): a network can represent more features than it has dimensions by giving each feature an almost-orthogonal direction. The price is interference $\langle w_i, w_j\rangle \ne 0$. It pays off when features are sparse (rarely active together), and a nonlinearity such as ReLU filters out the small interference. A consequence: features need not align with single neurons. This is why interpretability work uses sparse autoencoders to find them.

Interactive figure: loading…

Answers to section 1. (1) No. $k$ depends only on $n$ and $\varepsilon$: with $\varepsilon = 0.2$, $k \ge 4\ln(10^4)/(0.02 - 0.00267) \approx 2126$ (about 2,660 for a 99% guarantee), the same for $D = 10^6$ or $D = 10^{12}$. (2) $\gamma = 0.25$, so about $[(1-0.5)^2, (1+0.5)^2] = [0.25, 2.25]$: a factor 9 between the smallest and the largest, although all true eigenvalues are 1.

4. Worked example #

(a) JL dimension. $n = 10^4$ points, $\varepsilon = 0.2$. Then $\varepsilon^2/2 - \varepsilon^3/3 = 0.02 - 0.002667 = 0.017333$ and $4\ln(10^4) = 36.84$, so $k \ge 36.84/0.017333 = 2125.5$: $k = 2126$. This is the bound without the $\delta$ term (the form scikit-learn uses). With the exact number $n(n-1)/2$ of pairs it still guarantees success with probability at least $1/n$, and in practice the union bound is very pessimistic. For a real guarantee, $\delta = 0.01$ adds $2\ln 100 = 9.21$ to the numerator: $k \ge 2657$. With $\varepsilon = 0.1$ and $n = 10^6$: $k \ge 11{,}842$. The $\varepsilon^{-2}$ is expensive; the $\log n$ is cheap.

(b) A tiny CountSketch. Fit a constant: $A = (1,1,1,1)^\top$, $b = (1,2,3,6)^\top$. The exact answer is the mean, $x^* = 3$, with $\|r^*\|^2 = 4 + 1 + 0 + 9 = 14$. Hash rows 1, 3 to bucket 1 and rows 2, 4 to bucket 2, with signs $(+, +, -, +)$:

The sketched problem $\min_x\,(0\cdot x + 2)^2 + (2x - 8)^2$ gives $\tilde x = 4$, and $\|A\tilde x - b\|^2 = 9 + 4 + 1 + 4 = 18$: the residual is $\sqrt{18/14} = 1.13$ times the optimum. This $S$ distorts $\operatorname{range}([A\;b])$ by $\varepsilon = 0.79$, so Theorem 1 only promises $\sqrt{1.79/0.21} = 2.9$. Two buckets for a 2-dimensional subspace is far too few; real sketches use $k \gg d$.

(c) Noise or signal? $n = 1000$ samples, $p = 250$ features, $\sigma^2 = 1$, so $\gamma = 0.25$ and $\sqrt\gamma = 0.5$.

• Noise bulk: eigenvalues of $\hat\Sigma$ in $[0.25, 2.25]$; singular values of $X$ in $\sqrt{1000}\,[0.5, 1.5] = [15.8, 47.4]$.
• Spike $\theta = 3$ (true variance 4 along $v$): $3 > 0.5$, so it is detected. The top sample eigenvalue is about $4(1 + 0.25/3) = 4.33$, not 4, and $\langle\hat v, v\rangle^2 \approx (1 - 0.25/9)/(1 + 0.25/3) = 0.897$, an angle of about 19°. Going back (step 5 of the step list): $\ell = 4.33$, $\ell - 1 - \gamma = 3.083$, and $\theta = \tfrac12\big(3.083 + \sqrt{3.083^2 - 1}\big) = \tfrac12(3.083 + 2.917) = 3$. ✓
• Spike $\theta = 0.4$ (true variance 1.4): $0.4 < 0.5$, so it is invisible. The top eigenvalue stays at about 2.25, and $\hat v$ is unrelated to $v$.
• So with 1000 samples, a direction with 40% more variance than the others cannot be found by PCA at all. Four times more samples ($\gamma = 0.0625$, threshold 0.25) would reveal it.

5. In code #

NumPy, SciPy and scikit-learn check the claims of sections 3 and 4. The notebook for this lesson has these cells (at smaller sizes) and the exercises.

import numpy as np
from scipy.spatial.distance import pdist
from scipy.sparse.linalg import LinearOperator, lsqr
from scipy.linalg import solve_triangular
from sklearn.random_projection import GaussianRandomProjection, johnson_lindenstrauss_min_dim

rng = np.random.default_rng(0)

# 1) JL: 500 points in D = 10_000 dimensions, eps = 0.25.
n, D, eps = 500, 10_000, 0.25
k = johnson_lindenstrauss_min_dim(n, eps=eps)     # the bound of Theorem 3 with delta = 1
X = rng.standard_normal((n, D))
Y = GaussianRandomProjection(n_components=k, random_state=0).fit_transform(X)
ratio = pdist(Y, "sqeuclidean") / pdist(X, "sqeuclidean")
print(k, ratio.min().round(3), ratio.max().round(3))   # 954 0.814 1.246: inside [0.75, 1.25]

# 2) Sketch-and-solve vs sketch-and-precondition, kappa(A) = 1e6.
m, d, ks = 20_000, 50, 500
A = rng.standard_normal((m, d)) @ np.diag(np.logspace(0, -6, d))
b = A @ rng.standard_normal(d) + 1e-3 * rng.standard_normal(m)
x_opt = np.linalg.lstsq(A, b, rcond=None)[0]
S = rng.standard_normal((ks, m)) / np.sqrt(ks)          # dense Gaussian sketch
x_sk = np.linalg.lstsq(S @ A, S @ b, rcond=None)[0]
res = lambda x: np.linalg.norm(A @ x - b)
rel = lambda x: np.linalg.norm(x - x_opt) / np.linalg.norm(x_opt)
print(res(x_sk) / res(x_opt), rel(x_sk))   # 1.044 (theory: E ratio^2 = 1 + d/(k-d-1) = 1.11), 1.73
R = np.linalg.qr(S @ A, mode="r")
AR = LinearOperator((m, d), matvec=lambda y: A @ solve_triangular(R, y),
                    rmatvec=lambda z: solve_triangular(R, A.T @ z, trans="T"))
print(np.linalg.cond(A), np.linalg.cond(solve_triangular(R, A.T, trans="T").T))  # 1e6 -> 1.9
y, istop, iters = lsqr(AR, b, atol=1e-14, btol=1e-14, iter_lim=500)[:3]
print(iters, rel(solve_triangular(R, y)))                   # 28 iterations, error ~3e-12
print(lsqr(A, b, atol=1e-14, btol=1e-14, iter_lim=500)[2])  # 500: no preconditioner, no convergence

# 3) Marchenko-Pastur: pure noise, then one spike of strength theta.
n, p = 1000, 250
gamma = p / n
Z = rng.standard_normal((n, p))
ev = np.linalg.eigvalsh(Z.T @ Z / n)
print(ev.min().round(3), ev.max().round(3), (1 - gamma**0.5)**2, (1 + gamma**0.5)**2)
v = np.zeros(p); v[0] = 1.0
for theta in [3.0, 0.4]:
    Xs = Z + (np.sqrt(1 + theta) - 1) * np.outer(Z @ v, v)    # covariance I + theta v v^T
    w, V = np.linalg.eigh(Xs.T @ Xs / n)
    print(theta, w[-1].round(3), (V[:, -1] @ v)**2)   # 3.0: ~4.4, ~0.91;  0.4: ~2.24, ~0

# 4) Rank-3 signal + noise (500 x 500, sigma = 1): which singular values are signal?
N = 500
U0 = np.linalg.qr(rng.standard_normal((N, 3)))[0]
V0 = np.linalg.qr(rng.standard_normal((N, 3)))[0]
Yn = U0 @ np.diag([60.0, 30.0, 18.0]) @ V0.T + rng.standard_normal((N, N))
s = np.linalg.svd(Yn, compute_uv=False)
print(s[:5].round(1), round(2 * N**0.5, 1), round(4 / 3**0.5 * N**0.5, 1))
# [67.1 47.2 44.  43.7 43.4]  edge 44.7  Gavish-Donoho 51.6

# 5) Near-orthogonality: 4000 random unit vectors in R^1000.
W = rng.standard_normal((4000, 1000))
W /= np.linalg.norm(W, axis=1, keepdims=True)
G = W @ W.T
np.fill_diagonal(G, 0)
print(np.abs(G).max().round(3), np.sqrt(4 * np.log(4000) / 1000).round(3))  # 0.157 <= 0.182

Part 2: the residual is near-optimal but $x$ is off by 170%, because $\kappa(A) = 10^6$; the same sketch as a preconditioner gives full accuracy. Part 4: the signal singular values $60, 30, 18$ are $2.68, 1.34, 0.80$ in units of $\sqrt N\sigma = 22.4$, so theory predicts observed values $68.3$, $46.7$ and "at the edge" ($0.80 < 1$). The bulk-edge rule keeps two components, Gavish–Donoho keeps one, and the third is undetectable.

6. Where it shows up #

• Large least squares and low-rank approximation. Sketch-and-precondition solvers (Blendenpik, LSRN) for $n \gg d$; the randomized SVD (lesson 09) sketches the range with $A\Omega$; streaming PCA keeps only a sketch.
• Embeddings and gradients. Random projection before nearest-neighbor search; TRAK (Park et al. 2023) compresses per-example gradients with a JL projection for data attribution.
• Ensemble data assimilation. With $N = 50$ members for a state of size $10^6$, the sample covariance has rank at most 49 and spurious long-range correlations. Localization, inflation and shrinkage (Ledoit–Wolf) are the fixes.
• Choosing the number of EOFs or POD modes. Compare the singular values with the noise bulk edge, or use the Gavish–Donoho threshold, instead of an "elbow" by eye.
• Neural networks. For a random initial weight matrix $W$, the eigenvalues of $W^\top W$ (the squared singular values) follow the MP law; trained layers show outliers above the bulk and heavy tails, which carry learned structure (Martin and Mahoney, JMLR 2021). Superposition motivates sparse autoencoders in interpretability.

7. Common mistakes #

• Using sketch-and-solve when you need $x$. It controls the residual only. For an accurate $x$, sketch-and-precondition.
• A dense Gaussian sketch on dense data. Forming $SA$ costs $O(knd)$, more than the $O(nd^2)$ of Householder QR when $k > d$. Use SRHT or CountSketch for speed; Gaussian is for analysis and for small problems.
• Forgetting the $1/\sqrt k$ scaling, or applying the JL bound for $n$ points to a whole subspace (infinitely many points): a subspace needs $k \propto d$.
• Believing $k$ is small for small $\varepsilon$. $\varepsilon = 0.01$ costs $\approx 8\cdot10^4\ln n$ dimensions. JL helps when $D$ is huge; for 768-dimensional embeddings, PCA or quantization is often better.
• Reading sample eigenvalues as population eigenvalues when $p/n$ is not small, or mixing up $\gamma = p/n$ with $n/p$.
• Treating a value just above $\lambda_+$ as certain signal. Finite-$n$ (Tracy–Widom) fluctuations, a wrong $\sigma$, or correlated noise move the edge.

8. Check yourself #

1. Use the bound of Theorem 3 for $n = 10^6$ points and $\varepsilon = 0.1$ ($\delta$ not counted). Does the answer change for $D = 10^4$?
   Answer

   $k \ge 4\ln(10^6)/(0.005 - 0.000333) = 55.26/0.004667 \approx 11{,}842$. It does not depend on $D$, so for $D = 10^4$ JL gives no reduction at all: use a larger $\varepsilon$, or PCA.
2. A sketch has $\|I - U^\top S^\top SU\|_2 = 0.2$ on $\operatorname{range}([A\;b])$. Bound the sketch-and-solve residual, and $\kappa(AR^{-1})$ for $SA = QR$.
   Answer

   Both by $\sqrt{1.2/0.8} = \sqrt{1.5} \approx 1.225$: the residual is at most 22.5% above the optimum, and $\kappa(AR^{-1}) \le 1.225$.
3. Why does the JL dimension grow like $\log n$ and not like $n$?
   Answer

   One distance fails with probability $\le 2e^{-ck\varepsilon^2}$ (exponential tail of $\chi^2_k$). A union bound over fewer than $n^2/2$ pairs multiplies this by $n^2$, which is absorbed when $ck\varepsilon^2 \ge 2\ln n + \ln(1/\delta)$.
4. $n = 2000$ samples, $p = 500$ features, noise variance 1. Where is the noise bulk of $X^\top X/n$? Is a spike $\theta = 1$ visible? Give the predicted top eigenvalue and squared overlap.
   Answer

   $\gamma = 0.25$: bulk $[0.25, 2.25]$. $\theta = 1 > \sqrt\gamma = 0.5$, so yes: top eigenvalue $\to 2(1 + 0.25) = 2.5$, and $\langle\hat v, v\rangle^2 \to (1 - 0.25)/(1 + 0.25) = 0.6$.
5. $Y = X_0 + \sigma Z$ is $1000\times1000$ with $\sigma = 0.1$. Give the bulk edge and the Gavish–Donoho threshold. Do you keep an observed singular value of 7.0?
   Answer

   Edge $2\sqrt{1000}\cdot0.1 = 6.32$; $\tau^* = 2.309\cdot\sqrt{1000}\cdot0.1 = 7.30$. The value 7.0 is clearly above the edge, so it is almost surely real signal, but it is below $\tau^*$: for the smallest mean squared error, drop it (its singular vectors are too noisy to help).
6. How can a 1000-dimensional residual stream hold $10^5$ features with little interference?
   Answer

   Random directions are almost orthogonal: $10^5$ random unit vectors in $\mathbb{R}^{1000}$ have $\max\lvert\cos\rvert \lesssim \sqrt{4\ln 10^5/1000} \approx 0.21$. If the features are sparse, few are active at once, so the interference stays small, and a nonlinearity can remove it.

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. Define an $\varepsilon$-subspace embedding $S$ for the subspace with orthonormal basis $U$.

   Show answer

2. What does sketch-and-solve least squares guarantee, and what does it not guarantee?

   Show answer

3. What is the idea of sketch-and-precondition?

   Show answer

4. State the Johnson–Lindenstrauss lemma (the size of $k$).

   Show answer

5. Why does the JL dimension grow like $\log n$?

   Show answer

6. Pure noise, variance $\sigma^2$, $\gamma = p/n \le 1$: where are the eigenvalues of $X^\top X/n$?

   Show answer

7. When can PCA detect a spike $\theta$ (covariance $\sigma^2(I + \theta vv^\top)$), and where is the top eigenvalue then?

   Show answer

8. How long is a random Gaussian vector $g \sim N(0, I_d)$?

   Show answer

9. Compute the typical spread of the angle between two random directions in $\mathbb{R}^{10\,000}$.

   Show answer

10. Rank the Gaussian, SRHT and CountSketch by the cost of computing $SA$ for $A \in \mathbb{R}^{n\times d}$.

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "State and prove JL." — $\chi^2_k$ concentration for one vector by Chernoff, then a union bound over pairs: $k = O(\varepsilon^{-2}\log n)$, independent of $D$.
• "How do you solve least squares with $10^9$ rows and 100 columns?" — Sketch (SRHT or CountSketch) to a few times $d$ rows, QR of $SA$ as a preconditioner, then LSQR: about $O(\operatorname{nnz}(A)\log(1/\text{tol}))$ work plus the small QR, and full accuracy.
• "Random projection or PCA?" — JL is oblivious, cheap and keeps every pairwise distance up to $1\pm\varepsilon$, with $k$ set by $n$. PCA adapts to the data and is optimal on average for rank $k$, but gives no worst-case guarantee for single pairs.
• "Is this eigenvalue of my covariance signal?" — Compare with the MP edge $\sigma^2(1+\sqrt{p/n})^2$; spikes below $\sqrt{p/n}$ are undetectable, and detected ones are biased upward.

10. Practice #

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

1. Write jl_min_dim(n, eps, delta) from Theorem 3 and compare it with scikit-learn.
2. Write a Gaussian random projection and measure the distortion of all pairwise distances against $\sqrt{2/k}$ and the JL $\varepsilon$.
3. Sketch-and-solve least squares: check the residual factor of Theorem 1 and see how far $\tilde x$ is from $x^*$ for an ill-conditioned $A$.
4. Sketch-and-precondition: build $R$ from the QR of $SA$, check $\kappa(AR^{-1})$, and solve to full accuracy with LSQR.
5. Signal or noise: compute the Marchenko–Pastur edges and the BBP prediction, and count the singular values of a noisy low-rank matrix above the bulk edge.

11. Watch and read #

Class 8, Video 1: Johnson-Lindenstrauss LemmaMary Wootters · Stanford, statement and proof of the JL lemma · 20 min

Follow the proof: one vector, the chi-squared tail, then the union bound over pairs (Theorem 3).

What is Random Matrix Theory, Visually ExplainedJingyuan Hu · noise spectra, Marchenko–Pastur and the BBP transition · 9 min

Watch the noise bulk appear and a spike leave it; compare with the figure in section 2.

RM+ML: 11. What Does a Large Random Covariance Matrix Look Like? — Marchenko–Pastur LawRoland Speicher · Random matrices and machine learning, lecture 11 · 59 min

A full lecture on Theorem 4: where the Marchenko–Pastur density and its edges come from.

Randomized Numerical Linear Algebra: OverviewPetros Drineas · Simons Institute (sketching for matrix multiplication, least squares and low rank) · 31 min

Sketching for least squares and low rank, the topic of Theorems 1 and 2.

How might LLMs store facts | Deep Learning Chapter 73Blue1Brown · from 17:04, superposition and almost-orthogonal directions · 22 min

From 17:04, almost-orthogonal directions in high dimension and superposition.

• Roman Vershynin, High-Dimensional Probability, 2nd ed. (free PDF): §3.1 (concentration of the norm), §5.3 (Theorem 5.3.1, the JL lemma, and Remark 5.3.4 on optimality), §4.6 and §7.3 (Theorems 4.6.1 and 7.3.1, singular values of random matrices).
• David Woodruff, Sketching as a Tool for Numerical Linear Algebra (2014): §2.1 (Definitions 1–2, Theorems 6, 7, 9), §2.5 (Theorem 21, regression) and §2.6 (sketch-and-precondition).
• Gavish and Donoho, The Optimal Hard Threshold for Singular Values is 4/√3 (IEEE Trans. Inf. Theory 2014): §I (the bulk edge and the threshold), Eq. (11) and Table I for $\lambda^*(\beta)$.
• Elhage et al., Toy Models of Superposition (Anthropic, 2022): "Definitions and Motivation" (the superposition hypothesis) and "Demonstrating Superposition".