Fast and Sparse Johnson-Lindenstrauss Transforms

Let $k = c\,\varepsilon^{-2}\log n$ for a suitable global constant $c$, and let $\Phi = PHD$, where:

• $D$ is a $d \times d$ diagonal matrix with independent uniform $\pm 1$ entries;
• $H$ is the normalized $d \times d$ Walsh-Hadamard matrix, $H_{ij} = d^{-1/2}(-1)^{\langle i-1,\, j-1\rangle}$, with the inner product of binary expansions taken mod 2;
• $P$ is $k \times d$ with independent entries, each $0$ with probability $1-q$ and $\mathcal{N}(0, q^{-1})$ with probability $q$, where $q = \min\{\Theta(\log^2 n / d),\ 1\}$.

With probability at least $2/3$ both of the following hold:

1. $(1-\varepsilon)\|x\|_2 \le k^{-1/2}\|\Phi x\|_2 \le (1+\varepsilon)\|x\|_2$ for every $x \in X$;
2. computing $\Phi x$ takes $O\big(d \log d + \min\{d\,\varepsilon^{-2}\log n,\ \varepsilon^{-2}\log^3 n\}\big)$ operations.

This is the Euclidean case $p = 2$ of Lemma 1 in Ailon and Chazelle (SIAM J. Comput. 2009), first presented at STOC 2006. The same lemma with $q = \min\{\Theta(\varepsilon^{-1}\log n / d), 1\}$ embeds into $\ell_1^k$ in $O(d\log d + \min\{d\,\varepsilon^{-2}\log n,\ \varepsilon^{-3}\log^2 n\})$ operations. The paper writes part 1 with a normalizing factor $\alpha_2$; the factor $k^{-1/2}$ used here is the one fixed by $\mathbb{E}\|Pu\|_2^2 = k\|u\|_2^2$. Section 2 of the paper assumes $n \ge d$ throughout; without it, the flattening step below gives the level $\sqrt{\log(nd)/d}$ in place of $\sqrt{\log n/d}$. When $X$ is known, success can be tested, and $O(\log(1/\delta))$ independent repetitions raise the success probability to $1 - \delta$.

$HD$ is orthogonal, so it changes no norm. A vector and its Walsh-Hadamard transform cannot both be concentrated on few coordinates, and the random signs in $D$ make it unlikely that $HDx$ is concentrated for any of the $n$ fixed inputs. After this rotation every input looks like the flat vector from the previous section, which is the good case for the sparse projection $P$.

Whenever $\varepsilon^{-2}\log^3 n \le d\log d$, one point costs $O(d \log d)$ instead of $O(\varepsilon^{-2} d\log n)$. Ailon and Chazelle used it for approximate nearest neighbors: for $n$ points in $\ell_2^d$ and any $\varepsilon > 0$ they build a randomized structure of size $n^{O(\varepsilon^{-2})}$ that answers $\varepsilon$-approximate nearest neighbor queries in $O(d\log d + \varepsilon^{-3}\log^2 n)$ time with high probability (Theorem 1 of the journal version).

The $O(d\log d)$ cost holds only while the second term is smaller. At fixed $\varepsilon$, $\varepsilon^{-2}\log^3 n \le d\log d$ means $\log n$ is at most of order $(d\log d)^{1/3}$; Krahmer and Ward summarize the construction as fast only for about $e^{d^{1/3}}$ points. The $d\log d$ term is paid in full even when $\mathrm{nnz}(x) = 1$, which rules the construction out for streaming updates. The guarantee is for a set fixed before the randomness is drawn, with success probability $2/3$. The Walsh-Hadamard matrix needs $d$ to be a power of two, so inputs are padded with zeros.

Fix a constant $\gamma > 0$. For all $d$ and $k$ with $k \le d^{1/2-\gamma}$ there is a distribution over $k \times d$ matrices $A$ with this Johnson-Lindenstrauss property such that $x \mapsto Ax$ can be computed in $O(d\log k)$ time for every $x \in \mathbb{R}^d$. Sampling $A$ uses $O(d)$ random bits, and the same holds for embeddings into $\ell_1$. The constants depend on $\gamma$ (Ailon and Liberty, Discrete Comput. Geom. 2009; statement and proof constants read in the author-hosted conference version).

With $k = C\varepsilon^{-2}\log n$ and a union bound, this gives an $\varepsilon$-JL map for a fixed $n$-point set at cost $O(d \log k)$ per point, provided $C\varepsilon^{-2}\log n \le d^{1/2-\gamma}$.

Two substitutions. The random sparse $P$ becomes a deterministic code matrix that is 4-wise independent: in any 4 rows, the entries of a uniformly random column are independent random signs, up to scaling. That makes its $2 \to 4$ operator norm small. The $\ell_\infty$ flatness of the input becomes $\ell_4$ flatness, which a few rounds of randomized Hadamard transforms deliver.

Up to $k \le d^{1/2-\gamma}$ this removes the $k$-dependent term of the original fast transform and beats both it and the dense $O(kd)$ map. For very small $k$ the $\log k$ factor replaces $\log d$.

The cap $k \le d^{1/2-\gamma}$ is binding. In the worked example below, $k \approx 1382$ exceeds $d^{1/2} = 1024$, so the theorem does not apply. The constant in the $\ell_4$ bound grows exponentially in $1/\gamma$.

Fix $\eta > 0$, $\varepsilon \in (0,1)$, and a finite set $E \subset \mathbb{R}^d$ with $|E| = n$. Let $r \ge 40\log(4n/\eta)$ and let $\Phi \in \mathbb{R}^{k \times d}$ have the restricted isometry property of order $r$ with constant $\delta_r \le \varepsilon/4$. Let $\xi$ be uniform on $\{-1,1\}^d$ and $D_\xi = \mathrm{diag}(\xi)$. With probability exceeding $1 - \eta$,

(Krahmer and Ward, SIAM J. Math. Anal. 2011, Theorem 3.1.)

Sort the coordinates of $x$ by magnitude. The few large coordinates form a sparse vector, which RIP handles deterministically. The remaining coordinates are small, and their interaction with the large block and with each other is a sum with random signs, which concentrates.

Every RIP bound becomes a JL bound. A random partial Hadamard or partial Fourier matrix has RIP of order $r$ and constant $\delta$ once $k \gtrsim \delta^{-2} r \log^4 d$. With random signs and $E$ the set of pairwise differences, $k \gtrsim \varepsilon^{-2}\log(n/\eta)\log^4 d$ rows give an $\varepsilon$-JL map computable in $O(d\log d)$ time. Krahmer and Ward call the construction unrestricted: $n$ has no cap of the form $\exp(d^{1/3})$ or $\exp(d^{1/2})$. More than $d$ rows saves nothing over the identity map, so $k \le d$ still limits $n$ to $\exp(\tilde O(\varepsilon^2 d))$. This improves the $\varepsilon^{-4}$ of Ailon and Liberty to $\varepsilon^{-2}$, the optimal dependence up to the $\log^4 d$ factor.

The signs must be independent of $E$. A fixed $\Phi$ with $k < d$ has a nontrivial null space and maps those vectors to $0$, so RIP alone never implies JL for an arbitrary set. The polylogarithmic factor in $d$ inherited from structured RIP bounds remains.

For $0 < \varepsilon, \delta < 1/2$ there are $k = \Theta(\varepsilon^{-2}\log(1/\delta))$, $s = \Theta(\varepsilon^{-1}\log(1/\delta))$ and a distribution over $k \times d$ matrices $S$ with exactly $s$ nonzero entries in every column, each equal to $\pm s^{-1/2}$, such that for every $x \in \mathbb{R}^d$

Two placements work. In the graph construction each column puts its $s$ nonzeros in $s$ distinct random rows. In the block construction the $k$ rows are split into $s$ blocks of $k/s$ rows and each column picks one random row per block. Hash locations and signs need only $O(\log(1/\delta))$-wise independence (for the graph construction, almost $O(\log(1/\delta))$-wise independent permutations), and one of the two can be sampled with $O(\log(1/\delta)\log d)$ random bits. Computing $Sx$ costs $O(s \cdot \mathrm{nnz}(x))$ (Kane and Nelson, J. ACM 2014; conference version SODA 2012).

With $\delta < 1/\binom{n}{2}$: $k = O(\varepsilon^{-2}\log n)$, $s = O(\varepsilon^{-1}\log n)$, and one point costs $O(\varepsilon^{-1}\log n \cdot \mathrm{nnz}(x))$.

Write $S_{ri} = \eta_{ri}\sigma_{ri}/\sqrt{s}$ with $\eta_{ri} \in \{0,1\}$ the nonzero pattern and $\sigma_{ri}$ the signs. Because every column has exactly $s$ entries of size $s^{-1/2}$, the diagonal part of $\|Sx\|_2^2$ equals $\|x\|_2^2$ exactly. All error comes from collisions, pairs $i \ne j$ that share a row:

In a turnstile stream, an update $x_j \leftarrow x_j + \Delta$ changes the sketch by $\Delta \cdot Se_j$, which costs $s = O(\varepsilon^{-1}\log(1/\delta))$ instead of $k$ for a dense matrix or $\Omega(d\log d)$ for a Hadamard-based one. Kane and Nelson use this to speed up the streaming regression and low-rank approximation algorithms of Clarkson and Woodruff (their Section 6). The output lives in $\ell_2$, not behind a median estimator, so it can feed nearest-neighbor search or gradient-based training on hashed features.

The sparsity cannot drop much further. Kane and Nelson show that, for $k = C\varepsilon^{-2}\log(1/\delta)$ and $\delta$ below a constant depending on $C$, their graph and block constructions need $s = \Omega(\varepsilon^{-1}\log(1/\delta))$ even with fully random hash functions (Theorems 18 and 19), and that the Dasgupta-Kumar-Sarlós scheme needs $s = \Omega(\varepsilon^{-1}\log^2(1/\delta)/\log^2(1/\varepsilon))$ (Theorem 16). Nelson and Nguyen (STOC 2013) show that some $n$-point set forces many nonzeros per column in every JL matrix on it. Their Theorem 10 concerns $n$ unit vectors in $\mathbb{R}^k$ with at most $s$ nonzeros each and pairwise inner products at most $\varepsilon$ in absolute value. For $1/\sqrt{n} < \varepsilon < \varepsilon_0$, with $\varepsilon_0 \le 1/2$ a fixed constant, and $k < O(n/\log(1/\varepsilon))$, it gives $s = \Omega(\varepsilon^{-1}\log n/\log(k/\log n))$. When also $k \le \mathrm{poly}(1/\varepsilon)\log n$, this reads $s = \Omega(\varepsilon^{-1}\log n/\log(1/\varepsilon))$ (their Corollary 11), which matches Kane-Nelson up to $\log(1/\varepsilon)$. On a dense input, $\mathrm{nnz}(x) = d$ and the cost $O(\varepsilon^{-1}d\log n)$ is only a factor $\varepsilon^{-1}$ below the dense matrix.

There is an absolute constant $C > 0$ such that for every $d > 1$, every $0 < \varepsilon < 1/2$ and every $N > d^C$, some $N$-point set $X \subset \mathbb{R}^d$ has the following property: every linear map $A: \mathbb{R}^d \to \mathbb{R}^k$ with the $\varepsilon$-JL guarantee on $X$ has

(Larsen and Nelson, ICALP 2016, LIPIcs 55, 82:1-82:11.)

Include the basis vectors, which force every column of $A$ to have norm near $1$, so $\mathrm{tr}(A^\top A) \approx d$. Include many Gaussian vectors $g$. For any fixed matrix, $\|Ag\|_2^2 - \mathrm{tr}(A^\top A)$ fluctuates at scale $\|A^\top A\|_F$, and some $g$ hits a large fluctuation. The JL guarantee caps the fluctuation, which caps $\|A^\top A\|_F$, and a trace-to-Frobenius ratio bounds the rank from below.

Every known JL construction is linear, and linearity is required in turnstile streaming and compressed sensing. This theorem rules out a better linear map even with full knowledge of $X$, for example one chosen after a first pass over the data.

The hard instance here is nearly orthogonal, and nearly orthogonal sets admit nonlinear embeddings into $o(\varepsilon^{-2}\log n)$ dimensions when $\varepsilon < 2^{-\omega(\sqrt{\log n})}$ (Larsen and Nelson 2017, Section 1.1). So this theorem alone leaves room for nonlinear maps.

For all integers $n, d \ge 2$ and every $\varepsilon$ with

there is a set $X \subset \mathbb{R}^d$ with $|X| = n$ such that every map $f: X \to \mathbb{R}^k$ with the $\varepsilon$-JL guarantee on $X$ has

In particular, if also $\varepsilon \ge n^{-0.4999}$, then $\log(\varepsilon^2 n) \ge 0.0002\log n$ and $k = \Omega(\varepsilon^{-2}\log n)$ (Larsen and Nelson, FOCS 2017, pp. 633-638, Theorem 2).

Build many point sets that differ in which coordinates each point touches, in a way any JL map must preserve. If $k$ were small, the embedded sets could be written down in too few bits to tell all the original sets apart.

Distributional JL plus a union bound is optimal: no nonlinear, data-dependent embedding improves the dimension for worst-case inputs. The necessary condition $\varepsilon = \Omega(1/\sqrt{\min\{n,d\}})$ for any bound of this form follows from the upper bound $O(\min\{n,d,\varepsilon^{-2}\log n\})$, and the theorem matches it up to the $(\log n)^{0.5001}$ factor. Larsen and Nelson conjecture the optimal dimension is $\Theta(\min\{n, d, \varepsilon^{-2}\log(2 + \varepsilon^2 n)\})$ for all $n, d > 1$ and $0 < \varepsilon < 1$. The revised arXiv version (Section 1.2) reports that Alon and Klartag proved the lower bound over the full parameter range, removing the $(\log n)^{0.5001}$ condition.

The bound is about the worst $n$-point set. Sets of low dimension, or the simplex at very small $\varepsilon$, embed into fewer dimensions. The earlier hard instances were nearly orthogonal and could not have given this result, because of the Reed-Solomon embeddings noted above; the encoding argument needs the different family.