Symptom

You are building semantic search. Ten million documents, each a 1536-dimensional embedding, and queries need the twenty nearest neighbours in under 50 milliseconds.

The index is 1536 floats × 10M documents = 61 GB, which does not fit in the memory you have. Exact search means 10M distance computations of 1536 multiply-adds each, which is 15 billion operations per query and nowhere near 50 ms.

The obvious move is dimensionality reduction, and the obvious objection is that throwing away 90% of the dimensions must destroy the distances you are trying to rank by. It does not, the number of dimensions you need does not depend on 1536 at all, and the projection that works is a random matrix requiring no knowledge of your data whatsoever.

Statement

Johnson–Lindenstrauss (1984). Let $X$ be a set of $n$ points in $\mathbb{R}^d$ and let $0 < \epsilon < 1$. For

$$k = O\!\left(\frac{\log n}{\epsilon^2}\right)$$

there exists a linear map $F : \mathbb{R}^d \to \mathbb{R}^k$ such that for all $u, v \in X$:

$$(1-\epsilon)\|u - v\|^2 \le \|F(u) - F(v)\|^2 \le (1+\epsilon)\|u - v\|^2.$$

Moreover $F$ can be taken to be a random matrix with i.i.d. Gaussian entries, and it succeeds with probability at least $1 - 1/n$.

The target dimension $k$ depends on $\log n$ and on $\epsilon$. It does not depend on $d$. Ten million points at $\epsilon = 0.1$ need $k \approx 8\ln(10^7)/0.01 \approx 12{,}900$ with the textbook constant — worse than 1536, which is the honest reading of the bound at small $\epsilon$ — but at $\epsilon = 0.3$ it drops to about 1,400, and empirically the constant is far better than the proof’s, with $k = 128$ or 256 working well in practice.

Argument

The proof is Chernoff plus a union bound, which is why it belongs after the concentration post.

Let $F = \frac{1}{\sqrt{k}} A$ where $A$ is $k \times d$ with i.i.d. $\mathcal{N}(0,1)$ entries. Fix a vector $w = u - v$, and assume $|w| = 1$ by scaling.

Step 1: the expectation is right. Each coordinate of $Aw$ is $\sum_j A_{ij} w_j$, a sum of independent Gaussians, hence $\mathcal{N}(0, |w|^2) = \mathcal{N}(0,1)$. So $|Fw|^2 = \frac{1}{k}\sum_{i=1}^k Z_i^2$ with $Z_i$ i.i.d. standard normal — a chi-squared with $k$ degrees of freedom, scaled — and $\mathbb{E}[|Fw|^2] = 1$. The projection is unbiased for every vector, regardless of $d$.

Step 2: concentration. $|Fw|^2$ is an average of $k$ i.i.d. terms, so it concentrates. The chi-squared tail bound gives

$$\Pr\big[\big|\|Fw\|^2 - 1\big| > \epsilon\big] \le 2\exp\!\left(-\frac{k(\epsilon^2 - \epsilon^3)}{4}\right).$$

The exponent is linear in $k$ and quadratic in $\epsilon$, which is where $k = O(\log n/\epsilon^2)$ comes from.

Step 3: union bound. There are $\binom{n}{2} < n^2$ pairs. Taking $k \ge \frac{8\ln n}{\epsilon^2 - \epsilon^3}$ makes each failure probability at most $2n^{-4}$, so the total failure probability is at most $2n^2 \cdot n^{-4} = 2/n^2 < 1/n$. All pairs survive simultaneously with high probability, so such an $F$ exists. $\blacksquare$

The $\log n$ is where the leverage is. Going from a thousand points to ten million multiplies $\log n$ by about 3.2. Ten million to ten billion adds another factor of 1.5. The dimension you need grows almost not at all with the size of your dataset, which is the property that makes this deployable.

The $1/\epsilon^2$ is where the pain is. Halving the distortion quadruples the dimension. This is the binding constraint in practice and it is tight: Larsen and Nelson (2017) proved the bound optimal — for some point sets, no map into fewer than $\Omega(\log n / \epsilon^2)$ dimensions preserves all distances. There is no better projection waiting to be found.

Back to the search index. Project to $k = 128$ with a random Gaussian matrix. The index shrinks from 61 GB to 5 GB, which fits in memory. Each distance is 128 multiply-adds instead of 1536, a 12× reduction. Projecting a query costs one $128 \times 1536$ matrix-vector product, about 200k operations, negligible next to the search.

Then the crucial practical step: retrieve 200 candidates in the projected space, and re-rank those 200 with exact distances in the original space. The projection needs only to be good enough to get the true top-20 into the top-200, which is a far weaker requirement than $\epsilon$-preservation of every pair, and it is why $k = 128$ works when the bound demands thousands. Recall@20 in the high 90s is typical.

Cheaper projections. Achlioptas (2001) showed the Gaussian entries can be replaced by $\pm 1$ with probability $1/2$ each, or by ${-1, 0, +1}$ with probabilities ${1/6, 2/3, 1/6}$ — the same guarantee with a third of the arithmetic and no multiplications. The Fast JL transform (Ailon and Chazelle) uses a sparse matrix composed with a Hadamard transform to project in $O(d \log d)$ rather than $O(dk)$. The lemma is robust to how you generate the randomness, which is unusual and useful.

Where else this shows up. Random projection is the mechanism underneath locality-sensitive hashing for cosine similarity — SimHash is a sign-random-projection whose collision probability is a function of the angle. It gives the sketching step in streaming algorithms such as AMS and count-sketch. And it explains why random features approximate kernels: the same concentration argument bounds the error.

Forbids

Beating $\Omega(\log n / \epsilon^2)$ dimensions for worst-case point sets. Larsen–Nelson. The lemma is optimal, and a better general-purpose projection does not exist.

Preserving all pairwise distances in $o(\log n)$ dimensions. Even at constant distortion the $\log n$ is needed.

Expecting reduction to be free. $\epsilon$ distortion is real, and it means the ranking of two nearly-equidistant neighbours can flip. If exact ranking matters, re-ranking in the original space is not optional.

Using the bound’s constant as a design parameter. The proof’s $k$ is a sufficient worst-case value, often larger than $d$ itself for reasonable $\epsilon$. Choosing $k$ by measured recall is the correct engineering procedure.

Does not forbid

It does not require knowing the data, which is the property people find hardest to believe and which is the whole reason this is practical. The matrix is random and fixed before any point is seen. Unlike PCA there is nothing to fit, no pass over the data, and new points project with the same matrix — so the embedding is streaming-compatible and stable as the corpus grows, which PCA is not.

It does not preserve inner products or angles as tightly as norms, though it nearly does: the polarisation identity turns the norm guarantee into an additive $O(\epsilon)$ guarantee on inner products between unit vectors. Additive, not multiplicative, which matters when the true inner product is near zero.

It does not beat PCA on quality. For a fixed $k$, PCA minimises squared reconstruction error and typically preserves distances better on real data with low-rank structure. JL wins on cost, on being data-independent, and on giving a guarantee for every pair rather than in aggregate — and PCA offers no per-pair guarantee at all.

It does not solve nearest-neighbour search. It reduces the dimension; you still need an index. The curse of dimensionality bites at $k = 128$ too, which is why production systems use HNSW or IVF-PQ on top of the projection rather than instead of it.

It does not extend to $\ell_1$. The analogous result fails there (Brinkman–Charikar), and this is a genuine limitation rather than an open problem. JL is a Euclidean phenomenon, resting on the rotational invariance of the Gaussian.

Boundary

  • The guarantee is for a fixed point set chosen before the projection. Adversarially chosen new queries can be distorted more, which matters if an attacker influences your inputs.
  • The constant in the $O(\cdot)$ is large in the proof and small in practice. Take the theorem as licence to try small $k$, then measure recall. Do not size your index from the formula.
  • High probability, not certainty. A given random matrix can be bad, with probability $1/n$. Verify on a sample of pairs; regenerate if it fails.
  • Distances are preserved, not semantics. If the original embedding’s distances do not mean what you want, projecting preserves that faithfully.
  • Sparse inputs become dense. Projecting a sparse bag-of-words vector into $k$ dimensions destroys the sparsity, and for very sparse high-dimensional data that can increase memory use. Feature hashing is the alternative in that regime.

The idea to keep: high-dimensional geometry has far fewer degrees of freedom than it looks like. A cloud of $n$ points, however many coordinates you gave it, really lives in about $\log n$ dimensions as far as distances are concerned, and a random matrix finds those dimensions without looking.