Bulk Johnson-Lindenstrauss Lemmas

The Johnson-Lindenstrauss lemma [JL84] concerns the approximate preservation of distances in a finite point set in Euclidean space. Specifically, for a subset $X\subset\mathbb{R}^{D}$ of $N$ points and a tolerance $\epsilon\in(0,1)$, there exist $k\times D$ matrices $Z$ and a constant $\gamma(\epsilon)$ for which

$$(1-\epsilon)\left\lVert x-x^{\prime}\right\rVert_{2}\leq\sqrt{\frac{\gamma(\epsilon)}{k}}\left\lVert Z(x-x^{\prime})\right\rVert_{2}\leq(1+\epsilon)\left\lVert x-x^{\prime}\right\rVert_{2}$$

(JL)

holds for all pairs of points $x,x^{\prime}\in X$ simultaneously, with probability at least $1-\delta$, provided

$$k=O(\epsilon^{-2}\log(N^{2}/\delta))=:D_{JL}(N)=:D_{JL}.$$

The matrices $Z$ are drawn randomly, and much work has been done since the original paper to equip $Z$ with special properties, such as allowing fast matrix multiplication, preserving sparsity, restricting the matrix entries to discrete distributions, and so forth; see [Nel20] for a recent review. The matrix $Z$ provides a linear method for dimension reduction, which, at the very least, reduces the amount of space needed to store the dataset $X$ on the computer, provided one can work with approximate versions of the pairwise distances. One would expect that “robust” downstream algorithms that depend on distance data, now working on the pointset $ZX$, should still return answers that approximate those found on the original pointset $X$, but now in shorter time, with less memory needed, etc. People appreciate this lemma, in theory [Vem04], for the above reasons, and if the algorithm scales linearly in the ambient dimension, in time or in space, then processing $ZX$ instead of $X$ will yield a proportional improvement.

However, certain algorithms, including many nearest-neighbor algorithms, scale exponentially in the dimension, especially if they only use space linear in $N$, due to packing arguments. For comparison to brute force, the all pairs nearest-neighbor problem may already be solved in $O(DN^{2})$ time by scanning the points of $X$ with respect to their distance to each query $x$. If the algorithm scales like $DN\exp(D)$, this scaling is only an improvement for $D=\log(N)$, and even then one would really prefer $D=o(\log N)$, as $N^{2}$ is too expensive when $N$ is large. Consequently, if we use multiplication by $Z$ as a preprocessing step, the target dimension $D_{JL}=O(\epsilon^{-2}\log(N^{2}))$ is much too large to be practical, as $\exp(D_{JL})$ is now polynomial in $N$ with exponent at least $2/\epsilon^{2}$ with $\epsilon<1$.

Apart from searching for better algorithms, it is then natural to ask whether there are matrices $Z$ with target dimension $k$ much smaller than $D_{JL}$ that would still satisfy equation (JL) for all pairs of points of $X$. However, Larsen and Nelson [LN14] showed no such matrix exists for general datasets – a union of the standard simplex and a sufficient number of Gaussian vectors forces at least one distance to be stretched or compressed too much. On the other hand, if further assumptions are made on the pointset, smaller $k$ is possible, for instance, when $X$ already lies in a low-dimensional subspace [Sar06] or when its unit difference set $\widehat{X-X}$ has low Gaussian complexity [KM05], even while allowing many more points in the dataset.

If one considers a given algorithm “robust”, one would hope that failing to preserve a few distances should not markedly change the final output; though, one would ideally have to prove such behavior for that algorithm. One is then led to ask: what happens to the other distances between the points of $ZX$ when $k$ is smaller than $D_{JL}$? To be concrete, can we approximately preserve $(1-\eta)$ of all the distances, for some fraction $\eta\in(0,1/2)$, ideally with $k=o(D_{JL}(N))$? The results in this paper show this is possible when $\eta$ is not too small and the data has high or even moderate “intrinsic dimension”, in a sense to be defined later. As a preview, for data sampled i.i.d. like a random vector $\xi\in\mathbb{R}^{D}$, corollary 4.1.8 shows that if $\xi$ has i.i.d. coordinates (No moment assumptions are made.), we can take

$$k=\frac{C^{\prime}}{\epsilon^{2}}\left(\log(4e/\eta)\log(6De/\zeta)+\frac{\log(N^{2}/\delta)}{\eta D}\right)$$

for preserving $(1-\eta)(1-\zeta)$ of all pairwise distances, with probability at least $1-2\delta$ over the draw for $Z$ and $X$, provided $N=\Omega(\zeta^{-1}D\log(N/\delta))$ and $N=\Omega(D\log(6De/\zeta))$. The matrix $Z$ may have i.i.d. standard Gaussian, or in general, zero-mean unit-variance sub-gausssian coordinates. This estimate for $k$ is an “improvement” over the $D_{JL}$ target dimension as soon as $\eta D=\omega(1)$ or say a fractional power of $\log(N)$; “improvement” must be in quotes, as we have only guaranteed “most” distances are approximately preserved, not all. For general $\xi$, the $1/\eta D$ is replaced by a $1/\eta\widehat{r}$, with $\widehat{r}$ the intrinsic dimension $1/\left\lVert\mathbb{E}\widehat{y}\widehat{y}^{\top}\right\rVert$ of the unit normalized random vector $\widehat{y}=\widehat{\xi-\xi^{\prime}}$ for an independent copy $\xi^{\prime}$ of $\xi$. We have $1\leq\widehat{r}\leq D$ always, and we may estimate it using corollary 4.2.1.

Both of our main results – theorems 2.0.6 for general datasets and theorem 4.1.5 for i.i.d. samples – rely on a dual viewpoint for the dimension reduction problem, namely, instead of asking how $Z$ transforms the data $X$, we ask how $X^{\top}$ transforms the test matrix $Z^{\top}$; we can then exploit known properties of how general matrices act on $Z^{\top}$ or its columns. Standard results like the Hanson-Wright inequality may be viewed in this light, and we indeed do so in this paper. Treating $Z^{\top}$ as a test matrix with known properties has been done previously in randomized numerical linear algebra [HMT10]; however, unlike there, slow decay in singular values is actually a good case for us.

The rest of the paper is organized as follows. The main argument allowing us to quantify how many distances can be preserved is in section 2 and leads to our first theorem 2.0.6, which, by itself, only suggests smaller target dimensions $k$ may be possible by considering “small” batches of difference vectors. We then recall the Walecki construction in section 3, which gives us a way to choose these batches that works well for i.i.d. samples as well as the standard simplex. Section 4 then presents the rest of our results specializing to data sampled i.i.d. from a given distribution, which may just be a larger dataset. This section leads up to our second main theorem 4.1.5 allowing us to make $k$ depend on the intrinsic dimension $\widehat{r}$ mentioned above. We also show how to estimate $\widehat{r}$ from the data. The appendix contains proofs for the properties of $Z$ that we use in the paper, namely, a variant of the Hanson-Wright inequality, and a corresponding independent proof in the Gaussian case in order to have decent explicit constants for $k$.

In this paper, we only study dimension reduction matrices $Z:\mathbb{R}^{D}\to\mathbb{R}^{k}$ with isotropic rows, that is, $\mathbb{E}(Z_{i},y)^{2}=\left\lVert y\right\rVert_{2}^{2}$ for all $y$ and all rows $Z_{i}$. So, $\mathbb{E}\left\lVert Zy\right\rVert_{2}^{2}=k\left\lVert y\right\rVert_{2}^{2}$. We say more about isotropic random vectors in section 1.1.1. To guarantee equation (JL) holds for a difference vector $y\in Y$, the usual proof for the Johnson-Lindenstrauss lemma considers each vector individually, providing upper bounds for the failure probabilities

$$p_{+}(y):=\mathbb{P}_{Z}\left\{\left\lVert Zy\right\rVert_{2}^{2}>(1+\widetilde{\epsilon}_{+})k\left\lVert y\right\rVert_{2}^{2}\right\}$$

and

$$p_{-}(y):=\mathbb{P}_{Z}\left\{\left\lVert Zy\right\rVert_{2}^{2}<(1-\widetilde{\epsilon}_{-})k\left\lVert y\right\rVert_{2}^{2}\right\},$$

implicitly working with the formulation

$$(1-\widetilde{\epsilon}_{-})k\left\lVert x-x^{\prime}\right\rVert_{2}^{2}\leq\left\lVert Z(x-x^{\prime})\right\rVert_{2}^{2}\leq(1+\widetilde{\epsilon}_{+})k\left\lVert x-x^{\prime}\right\rVert_{2}^{2}$$

($JL^{2}$)

which is often easier to manage. In lemma 5.0.5 of the appendix, we show how to choose $\widetilde{\epsilon}_{\pm}$ in line with the original equation (JL).

Suppose $\widetilde{\epsilon}_{-}\leq\widetilde{\epsilon}_{+}$, as it will for this paper. If the distribution of $Z$ is sufficiently nice, sub-gaussian for example, then one may show $p_{+}+p_{-}\leq 2\exp(-\widetilde{\epsilon}_{-}^{2}k/C)$ for each fixed $y$, with $C$ a constant depending on the distribution of $Z$. By the union bound (Boole’s inequality), the probability that a random $Z$ fails for some $y\in Y$ is at most

$$\binom{N}{2}2\exp(-\widetilde{\epsilon}_{-}^{2}k/C)\leq N^{2}\exp(-\widetilde{\epsilon}_{-}^{2}k/C)<\delta,$$

provided $k>C\widetilde{\epsilon}_{-}^{-2}\log(N^{2}/\delta)=D_{JL}$. The probability there is some $Z$ that succeeds for all $y\in Y$ is then at least $1-\delta$, so at least one such $Z$ exists.

The argument above considers the vectors $y\in Y$ one at a time, making sure $Z$ succeeds for each $y$; if we consider the unit normalized vectors $\widehat{y}\in\widehat{Y}$, we may view this argument as controlling the extreme order statistics of the random variables $\left\lVert Z\widehat{y}_{i}\right\rVert_{2}^{2}$, induced by $Z$, namely

$$(1-\widetilde{\epsilon}_{-})k\leq\left\lVert Z\widehat{y}_{(0)}\right\rVert_{2}^{2}\leq\ldots\leq\left\lVert Z\widehat{y}_{(j)}\right\rVert_{2}^{2}\leq\ldots\leq\left\lVert Z\widehat{y}_{\big{(}\binom{N}{2}-1\big{)}}\right\rVert_{2}^{2}\leq(1+\widetilde{\epsilon}_{+})k.$$

If we only wish to preserve a fraction of the distances, say $(1-\eta)$ with $\eta$ hopefully small, we can consider controlling the intermediate or central order statistics of the $\left\lVert Z\widehat{y}_{i}\right\rVert_{2}^{2}$ instead. We do so as follows.

If we divide the difference vectors into batches of size $M$, and preserve $(1-\eta)M$ distances there, then we still recover

$$(1-\eta)M\binom{N}{2}/M=(1-\eta)\binom{N}{2}\quad\text{ of the total distances.}\quad$$

We assume $\eta M$ is a strictly positive integer here, and for simplicity of discussion, we also assume $M$ divides $N$; we shall return to this point later. Let $\Upsilon\subset Y$ be a given batch of size $M$. Each given matrix $Z$ also induces an ordering on the points of $\Upsilon$: with $\widehat{y}=y/\left\lVert y\right\rVert_{2}$,

$$\left\lVert Z\widehat{y}_{(0)}\right\rVert_{2}^{2}\leq\ldots\leq\left\lVert Z\widehat{y}_{(i)}\right\rVert_{2}^{2}\leq\ldots\leq\left\lVert Z\widehat{y}_{(M-1)}\right\rVert_{2}^{2}.$$

As $Z$ is random, this ordering is too, treating $Y$ as fixed. If we could guarantee that $\left\lVert Z\widehat{y}_{((1-\eta)M)}\right\rVert_{2}^{2}\leq(1+\epsilon)k$, then all $\left\lVert Z\widehat{y}_{(i)}\right\rVert_{2}^{2}$ with $i\leq(1-\eta)M$ also have this upper bound, with an analogous statement for a lower bound of $(1-\epsilon)k$ on $\left\lVert Z\widehat{y}_{(\eta M-1)}\right\rVert_{2}^{2}$, altogether guaranteeing $(1-2\eta)M+2>(1-2\eta)M$ of the vectors of $\Upsilon$ have

$$(1-\widetilde{\epsilon}_{-})k\leq\left\lVert Z\widehat{y}\right\rVert_{2}^{2}\leq(1+\widetilde{\epsilon}_{+})k.$$

We need to control the probabilities

	$\displaystyle p_{+}(\Upsilon):$	$\displaystyle=\mathbb{P}_{Z}\left\{\left\lVert Z\widehat{y}_{((1-\eta)M)}\right\rVert_{2}^{2}>(1+\widetilde{\epsilon}_{+})k\right\}$		(2.1)
and
	$\displaystyle p_{-}(\Upsilon):$	$\displaystyle=\mathbb{P}_{Z}\left\{\left\lVert Z\widehat{y}_{(\eta M-1)}\right\rVert_{2}^{2}<(1-\widetilde{\epsilon}_{-})k\right\}.$		(2.2)

We can recast control of $p_{\pm}(\Upsilon)$ using the following lemma. Let $\widehat{\Upsilon}=\left\{\widehat{y}_{0},\ldots,\widehat{y}_{M-1}\right\}$, that is, the unit normalized version of $\Upsilon$.

We now make assumptions on the matrix $Z$, allowing us to make use of the stable rank of the minibatches $\Upsilon(\eta)$.

One point we want to stress even now is the presence of the product $k\,r_{\infty}$ in the exponent. If one wants a fixed failure probability, $k$ need not be as large when $r_{\infty}$ is sizable. We shall give several examples in this paper where $r_{\infty}$ is large.

Consider the term $\left\lVert Z\Upsilon(\eta)\right\rVert_{F}^{2}$ in lemma 2.0.1, taking $\Lambda_{\Upsilon(\eta)}$ as the identity. (The following discussion will also hold for other scalings, such as $\widehat{\Upsilon}(\eta)$.) Labeling the rows of $Z$ as $Z_{i}$, we can interchange the sums implicit in $\left\lVert Z\Upsilon(\eta)\right\rVert_{F}^{2}$ to our advantage:

$\displaystyle\sum_{y\in\Upsilon(\eta)}\left\lVert Zy\right\rVert_{2}^{2}=\sum_{y\in\Upsilon(\eta)}\sum_{i=1}^{k}(Z_{i},y)^{2}=\sum_{i=1}^{k}\sum_{y\in\Upsilon(\eta)}(y,Z_{i})^{2}=\sum_{i=1}^{k}\left\lVert\Upsilon(\eta)^{\top}Z_{i}\right\rVert_{2}^{2}.$

Because we assume the rows $Z_{i}$ are isotropic, $\mathbb{E}\left\lVert\Upsilon(\eta)^{\top}Z_{i}\right\rVert_{2}^{2}=\left\lVert\Upsilon(\eta)^{\top}\right\rVert_{F}^{2}=\left\lVert\Upsilon(\eta)\right\rVert_{F}^{2}$, and we have transformed the bound on $p_{+}(\Upsilon)$ to involve the probabilities

$$\mathbb{P}_{Z}\left\{\sum_{i=1}^{k}\left\lVert\Upsilon(\eta)^{\top}Z_{i}\right\rVert_{2}^{2}>(1+\widetilde{\epsilon}_{+})\sum_{i=1}^{k}\left\lVert\Upsilon(\eta)\right\rVert_{F}^{2}\right\},$$

and $\Upsilon(\eta)^{\top}$ is now viewed as a fixed matrix acting on the random vectors $Z_{i}$. We may thus take a dual viewpoint on the dimension reduction problem: instead of considering how the matrix $Z$ acts on each difference vector, we consider how the transposed batch of difference vectors $\Upsilon(\eta)^{\top}$ acts on the matrix $Z^{\top}$, as mentioned in the introduction. If we take the $Z_{i}$ to be independent, with i.i.d. mean-zero unit-variance sub-gaussian entries, we can use lemma 5.0.3 in the appendix to harness both the independence of the $Z_{i}$ and the Hanson-Wright inequality:

We now can control the probabilities in equations (2.1) for a given batch $\Upsilon$ of size $M$. The control is in terms of $r_{\infty}(\Upsilon(\eta))$ or $r_{4}(\Upsilon(\eta))$ for subsets of size $\eta M$. Because the target dimension $k$ is a global parameter, it needs to be in terms of global quantitities, but the stable ranks above vary over minibatches. To make this transition and to help with bookkeeping, recall that a decomposition of a graph is a partition of its edges.

Let $\mathcal{P}$ be a decomposition of the complete graph on $N$ vertices, into batches $\Upsilon$ of size $M$. If $M$ does not divide $N$, that is, with $N=jM+n$, we can expand several of the batches to $M+s\geq M$ with $s=\left\lceil n/j\right\rceil<M$. For those batches, $\eta(M+s)$ need not be an integer, so take $\tilde{\eta}$ as

$$\tilde{\eta}=\eta\frac{M}{M+s}\quad\text{in order for}\quad\tilde{\eta}(M+s)=\eta M,$$

(2.3)

and set $\Upsilon(\eta)=\Upsilon(\tilde{\eta})$ when the batch size is not $M$. Note smaller $\eta$ values only help us ensure a total fraction $(1-\eta)$ of distances is preserved. For this decomposition, let

$$R_{\infty}(\eta M):=R_{\infty}(\eta M;\mathcal{P}):=\min_{\Upsilon\in\mathcal{P}}\min_{\Upsilon(\eta)\subset\Upsilon}r_{\infty}(\Upsilon(\eta))$$

be the minimum stable rank of the $\eta M$ sized minibatches from such batches. We then have our first theorem, written in terms of $\epsilon$ in the original equation (JL).

To make sense of the above, suppose $\eta M=D$, and then recall $1\leq R_{\infty}(\eta M)\leq D$ as the stable rank is always bounded above by the ambient dimension. If $R_{\infty}(\eta M)=cD$, for $c$ bounded away from 0, the term attached to $\log(2e/\eta)$ becomes constant, while the remaining becomes constant as soon as $D\gtrsim\log(2N^{2}/(D\delta))$.

Note $R_{\infty}(\eta M)$ depends on the decomposition $\mathcal{P}$, which we are free to choose at will. The remainder of the paper will address how to choose $\mathcal{P}$ (and the batch size $M$).

From lemma 2.0.1 we are free to scale vectors in each batch $\Upsilon$ to have unit norm. So we can also define for a given decomposition

$$\widehat{R}_{\infty}(\eta M;\mathcal{P})=\min_{\widehat{\Upsilon}\in\mathcal{P}}\min_{\widehat{\Upsilon}(\eta)\subset\widehat{\Upsilon}}r_{\infty}(\widehat{\Upsilon}(\eta)).$$

This normalization allows us to use linear algebra to control $r_{\infty}(\widehat{\Upsilon}(\eta))$ deterministically, for any of the $\eta M$-sized subsets of $\widehat{\Upsilon}$, once we have control of $\sigma_{1}(\widehat{\Upsilon})$. When the underlying data is random, we can thus avoid the need to take union bounds over the minibatches.