Here is one of the most audacious results in machine learning: you can slash a dataset from thousands of dimensions down to a few dozen by multiplying it with a matrix of random numbers — and the geometry survives. No looking at the data, no eigenvectors, no training. Just draw a random linear map R and send every point x to (1/√k)·R·x. The pairwise distances between points barely change.
The Johnson–Lindenstrauss lemma
The licence for this is the Johnson–Lindenstrauss lemma. For any set of n points, there exists a random linear map into k ≈ O(log n / ε²) dimensions that preserves every pairwise distance within a factor of (1 ± ε). The astonishing part is what's absent from that formula: the original dimension D. A 10,000-dimensional bag-of-words and a 1,000,000-dimensional one collapse to the same target size — k depends only on how many points you have and how much distortion you'll tolerate.
Why does a random direction preserve length? Because the projected squared length is an average of k independent contributions, each with mean equal to the true squared length. By concentration of measure it clusters tightly around the truth, and averaging more of them (bigger k) makes the estimate tighter — the distortion falls like 1/√k.
The entire "model" is a random matrix
// standard normal via Box–Muller
function randn(rng) {
let u = 0, v = 0;
while (u === 0) u = rng();
while (v === 0) v = rng();
return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}
// k×D Gaussian projection matrix, row-major, scaled 1/√k
function gaussianMatrix(k, D, rng) {
const R = new Float64Array(k * D);
const s = 1 / Math.sqrt(k);
for (let i = 0; i < k * D; i++) R[i] = s * randn(rng);
return R; // E‖Rx‖² = ‖x‖²
}
Nothing about your data is ever consulted — the map is chosen before the data is seen. The 1/√k factor is the calibration that makes the projection preserve lengths in expectation rather than inflating them by √k. Reducing a point is then a single matrix–vector multiply; for a whole dataset it's one (N×D)·(D×k) product. No fitting, no SVD, no iteration.
You don't even need real Gaussians
Achlioptas showed that entries drawn as √3 × {+1 with prob 1/6, 0 with prob 2/3, −1 with prob 1/6} give the same JL guarantee — yet two-thirds of the entries are zero (so you skip those multiplies) and the rest are just additions and subtractions.
function sparseMatrix(k, D, rng) {
const R = new Float64Array(k * D);
const s = Math.sqrt(3 / k); // √3 · (1/√k)
for (let i = 0; i < k * D; i++) {
const u = rng();
R[i] = u < 1/6 ? s : (u < 1/3 ? -s : 0); // ⅔ are zeros
}
return R;
}
Sizing k, and where it sits
The exact JL bound (the one scikit-learn ships) makes the constant explicit:
import numpy as np
def jl_min_dim(n, eps):
return int(np.ceil((4 * np.log(n)) / (eps**2/2 - eps**3/3)))
jl_min_dim(1_000_000, 0.1) # ~ 15,000 — independent of D
jl_min_dim(1_000_000, 0.3) # ~ 1,700
That bound is a worst-case guarantee holding for every pair simultaneously (a union bound over all n(n−1)/2 of them), so it's famously loose — the typical RMS distortion drops below ε at a far smaller k. In practice you reach for scikit-learn's SparseRandomProjection and johnson_lindenstrauss_min_dim, and the whole thing is a single blazing matmul.
Against the data-dependent reducers — PCA computes an SVD and keeps variance; t-SNE/UMAP preserve local neighbourhoods for visualisation — random projection makes the opposite bet: don't look at the data at all, just multiply by random noise. You lose PCA's optimality but gain a distance-preserving guarantee, a target size that ignores D, and speed. It even degrades gracefully — pick k too small and distances just blur, they never break. That's why it's the go-to front-end for kNN, k-means, and locality-sensitive hashing on enormous feature spaces.
Slide k, swap the matrix type, and watch the distortion histogram tighten around 1.0, live at: https://dev48v.infy.uk/ml/day54-random-projections.html
Top comments (0)