Johnson-Lindenstrauss Proof & Sparse Embeddings

Random Projections That Preserve Distance

Imagine photographing a three-dimensional sculpture. One photograph loses depth, yet several carefully chosen views preserve enough geometry to compare shapes. A random projection applies the same idea to data in (\mathbb R^d): it maps every vector to (\mathbb R^m), with (m\ll d), while approximately preserving the distances among a finite collection of points.

For a linear map (A\in\mathbb R^{m\times d}), the desired event for points (x_1,\ldots,x_N) is \[ (1-\varepsilon)\lVert x_i-x_j\rVert_2^2 \le \lVert A(x_i-x_j)\rVert_2^2 \le (1+\varepsilon)\lVert x_i-x_j\rVert_2^2 \] for every pair. The Johnson–Lindenstrauss lemma says a random map can achieve this with (m=O(\varepsilon^{-2}\log N)). The original dimension (d) does not appear in the target dimension bound.

For a tiny example, take (u=(1,0,1,0)), (v=(0,1,0,1)), so (u-v=(1,-1,1,-1)) and (\lVert u-v\rVert_2^2=4). If a particular projection gives (A(u-v)=(1.2,-1.5)), its squared norm is (3.69), a relative distortion of (-7.75\%\). This one calculation verifies one pair for one realized map; it does not prove that every pair or future point is preserved.

Concentration says a random quantity stays near its mean. With a Gaussian map (A_{rj}\sim N(0,1/m)), (\lVert Ax\rVert_2^2/\lVert x\rVert_2^2) is a normalized chi-square variable for fixed (x). A sub-Gaussian map replaces Gaussian entries with other light-tailed variables such as scaled random signs. An epsilon net is a finite cover used to extend guarantees from finitely many directions to every direction in a subspace. A sparse embedding gives each input column few nonzeros, reducing multiplication work. Feature hashing is the extreme one-nonzero-per-column version. An FJLT first randomizes and spreads energy with a fast orthogonal transform, then samples or projects to fewer coordinates.

This deep dive supplies what a dimensionality-reduction overview usually omits: the fixed-vector tail bound, the all-pairs union bound, and the reasons a faster sparse transform needs more careful parameter choices than a dense Gaussian proof.

Concentration for One Vector

Fix a nonzero vector (x\in\mathbb R^d) before drawing (A), with independent entries (A_{rj}\sim N(0,1/m)). Rotational invariance gives each row dot product (\langle A_r,x\rangle\sim N(0,\lVert x\rVert_2^2/m)). Therefore \[ Y=\frac{m\lVert Ax\rVert_2^2}{\lVert x\rVert_2^2} =\sum_{r=1}^m g_r^2\sim\chi_m^2, \qquad \mathbb E\lVert Ax\rVert_2^2=\lVert x\rVert_2^2. \]

For the upper tail, Markov's inequality on $e^{\lambda Y}$, where $0\lt\lambda\lt1/2$, uses the moment generating function $\mathbb E e^{\lambda Y}=(1-2\lambda)^{-m/2}$. Optimizing $\lambda$ gives \[ \Pr[Y\ge m(1+\varepsilon)] \le \exp\!\left[-\frac m2\bigl(\varepsilon-\ln(1+\varepsilon)\bigr)\right]. \] A parallel calculation for the lower tail yields an exponential bound. For $0\lt\varepsilon\lt1$, a convenient combined form is \[ \Pr\!\left[\left|\lVert Ax\rVert_2^2-\lVert x\rVert_2^2\right|\gt\varepsilon\lVert x\rVert_2^2\right] \le 2\exp\!\left[-m\left(\frac{\varepsilon^2}{4}-\frac{\varepsilon^3}{6}\right)\right]. \] Thus (m=\Theta(\varepsilon^{-2}\log(1/\delta))) preserves one fixed vector with failure probability at most (delta).

The proof tracks four obligations: isotropy makes the expectation correct; independence factors the moment generating function; light tails make it finite; and $x$ must be independent of the sampled map. If $x$ is chosen afterward, any nonzero nullspace vector of a dimension-reducing map has $Ax=0$, a decisive counterexample.

Scaled random-sign matrices satisfy comparable concentration through sub-Gaussian tools without rotational invariance. Sparse embeddings need separate analyses because collisions create dependencies and higher moments. An epsilon net is unnecessary for a finite set of pair differences, but it becomes essential when preserving every vector in a (k)-dimensional subspace; a net of size exponential in (k) turns the logarithm of its size into an (O(k)) dimension term.

flowchart TD
    A[Fix x independently of random map] --> B[Choose isotropic Gaussian rows]
    B --> C[Each row dot x is scaled normal]
    C --> D[Projected energy is normalized chi square]
    D --> E[Apply exponential moment bound]
    E --> F[Optimize lambda for upper and lower tails]
    F --> G[Failure decays as exp minus c epsilon squared m]
    G --> H[Choose m for target delta]
    H --> I[Union bound finite differences]
    I --> J[All pairwise distances preserved]
    A --> K{Vector chosen after seeing A?}
    K -->|yes| L[Nullspace vector destroys guarantee]
    K -->|no| B
    classDef navy fill:#17365D,color:#fff,stroke:#17365D;
    classDef teal fill:#147D84,color:#fff,stroke:#0D565B;
    classDef purple fill:#6C3BAA,color:#fff,stroke:#472675;
    classDef crimson fill:#C6284D,color:#fff,stroke:#861B35;
    class A,B,C navy;
    class D,E,F,G,H purple;
    class I,J teal;
    class K,L crimson;

From One Pair to Every Pair

Distances reduce to norms of differences: preserving $x_i-x_j$ preserves the distance between $x_i$ and $x_j$. There are $M=\binom N2$, fewer than $N^2/2$, nonzero pair differences. If the fixed-vector failure probability is at most $\delta/M$, the union bound gives total failure at most $\delta$, regardless of dependence among pair events. Consequently, \[ m\ge C\varepsilon^{-2}\bigl(\log N+\log(1/\delta)\bigr) \] rows suffice for an absolute constant (C). A conservative implementation may use a larger explicit constant derived from the selected tail inequality; asymptotic notation must not be mistaken for a production parameter.

An epsilon net solves a different quantifier problem. To preserve an entire (k)-dimensional subspace, cover its unit sphere by a finite net, apply concentration to each net point, then extend to all unit vectors via continuity and an operator-norm argument. The net has roughly ((1+2/\eta)^k) points, so its logarithm contributes (O(k)).

import math
import random

def gaussian_matrix(rows, cols, seed):
    if rows <= 0 or cols <= 0:
        raise ValueError("matrix dimensions must be positive")
    rng = random.Random(seed)
    scale = 1.0 / math.sqrt(rows)
    return [[rng.gauss(0.0, 1.0) * scale for _ in range(cols)]
            for _ in range(rows)]

def matvec(matrix, vector):
    if not matrix or any(len(row) != len(vector) for row in matrix):
        raise ValueError("incompatible matrix and vector")
    return [sum(a * b for a, b in zip(row, vector)) for row in matrix]

def squared_distance(a, b):
    if len(a) != len(b):
        raise ValueError("dimension mismatch")
    return sum((x - y) ** 2 for x, y in zip(a, b))

def distortion_report(points, projected):
    worst = 0.0
    ratios = []
    for i in range(len(points)):
        for j in range(i):
            before = squared_distance(points[i], points[j])
            after = squared_distance(projected[i], projected[j])
            if before == 0.0:
                assert after < 1e-20
                continue
            ratio = after / before
            ratios.append(ratio)
            worst = max(worst, abs(ratio - 1.0))
    return worst, ratios

rng = random.Random(41)
n_points, source_dim, target_dim = 24, 80, 48
points = [[rng.gauss(0.0, 1.0) for _ in range(source_dim)]
          for _ in range(n_points)]
points.append(points[0][:])
projection = gaussian_matrix(target_dim, source_dim, seed=7)
embedded = [matvec(projection, point) for point in points]
worst, ratios = distortion_report(points, embedded)
assert len(ratios) == len(points) * (len(points) - 1) // 2 - 1
assert min(ratios) > 0.0

zero = [0.0] * source_dim
assert matvec(projection, zero) == [0.0] * target_dim

try:
    matvec(projection, [1.0, 2.0])
    raise AssertionError("dimension mismatch accepted")
except ValueError:
    pass

bad_projection = [[1.0] + [0.0] * (source_dim - 1)]
bad_points = [[0.0] * source_dim, [0.0, 1.0] + [0.0] * (source_dim - 2)]
bad_embedded = [matvec(bad_projection, point) for point in bad_points]
bad_worst, _ = distortion_report(bad_points, bad_embedded)
assert bad_worst == 1.0

print("pairs checked:", len(ratios))
print("worst squared-distance distortion:", round(worst, 4))
print("mean ratio:", round(sum(ratios) / len(ratios), 4))
print("adversarial deterministic projection distortion:", bad_worst)

Sparse Embeddings and Fast Geometry

Sparse projections replace dense matrix-vector multiplication with a few signed writes per nonzero input feature. In a common construction, each input coordinate chooses (s) output rows and contributes a scaled signed value. The scaling preserves expected energy; multiple destinations reduce the chance that one collision dominates. Feature hashing sets (s=1), making application linear in the input's number of nonzeros, but arbitrary dense or spiky vectors can require more rows or stronger assumptions than dense JL.

FJLT takes another route: multiply coordinates by random signs, apply a normalized Hadamard or Fourier transform to flatten leverage, then subsample or apply a sparse projection. The spreading step prevents a single source coordinate from controlling one sampled output. It is fast for dense vectors and structured dimensions, though padding, memory layout, batching, and transform constants matter in practice.

In a teaching case, project a fixed point cloud with Gaussian, sign, and one-sparse hashing matrices, then plot the full distribution of pairwise distortion. In a professional nearest-neighbor pipeline, freeze and version the projection, validate recall after projection, and retain exact re-ranking on original vectors. In a scientific least-squares problem, use a subspace-embedding theorem rather than an all-pairs finite-set statement. Reject the method if queries adapt to a public map, if rare spiky vectors dominate, if downstream margins are smaller than allowed distortion, or if dimension reduction saves less than transform overhead.

flowchart TD
    A[Define protected points pairs or subspace] --> B[Set distortion epsilon and failure delta]
    B --> C{Dense or sparse inputs?}
    C -->|dense| D[Benchmark Gaussian sign and FJLT]
    C -->|sparse| E[Benchmark sparse JL and feature hashing]
    D --> F[Validate maximum distortion and task quality]
    E --> F
    F --> G{Guarantee and latency pass?}
    G -->|yes| H[Freeze seed scale dimension and version]
    G -->|no| I[Increase dimension change transform or reject]
    H --> J[Monitor drift recall and adaptive exposure]
    J --> K{Assumptions remain valid?}
    K -->|yes| J
    K -->|no| I
    classDef navy fill:#17365D,color:#fff,stroke:#17365D;
    classDef teal fill:#147D84,color:#fff,stroke:#0D565B;
    classDef purple fill:#6C3BAA,color:#fff,stroke:#472675;
    classDef crimson fill:#C6284D,color:#fff,stroke:#861B35;
    class A,B,C navy;
    class D,E,F purple;
    class G,H,J teal;
    class I,K crimson;

Prove and Test a JL Dimension

The following lab moves from the fixed-vector proof to production judgment. Keep squared norms and ordinary distances distinct: a (1\pm\varepsilon) guarantee on squared distance translates to square-root factors on distance.