Almost Optimal Unrestricted Fast Johnson-Lindenstrauss Transform

Introduction

Designing computationally efficient transformations that reduce dimensionality of data while approximately preserving its metric information lies at the heart of many problems. While in compressed sensing such techniques are sought for sparse data in a real or complex metric space (with respect to some basis), in random projections, following the seminal work of Johnson and Lindenstrauss, one seeks to reduce dimension of any set of finite data.The term ”random projections” describes Johnson and Lindenstrauss’ s original construction and became synonymous with the process of approximate metric preserving dimension reduction using randomized linear mappings. However, these linear mappings need not be (and indeed are usually not) projections in the linear algebraic sense of the word. In both applications, random matrices of a suitable size result in optimal construction in the parameters $n$ (the original dimension), $k$ (the target dimension), $N$ (the number of input vectors) and $\delta$ (the distortion). However, these constructions’ resulting running time complexity, measured as number of operations needed in order to map a vector, is suboptimal.

The transformation we derive here is a composition of two random matrices: A random sign matrix and a random selection of a suitable number $k$ of rows from a Fourier matrix, where $k=O(\delta^{-4}(\log N)\operatorname{polylog}(n))$, and $\delta$ is the tolerated distortion level. The result, for constant $\delta$, is believed to be suboptimal within the $\operatorname{polylog}(n)$ factor in the target dimension $k$. The running time of performing the transformation on a vector is dominated by the $O(n\log n)$ of the Fast Fourier Transform, and is believed to be optimal. The possibility of obtaining such a running time for fixed distortion was left as an open problem in Ailon and Chazelle and Ailon and Liberty’s work, and here we resolve it up to a factor of $\operatorname{polylog}(n)$. The dependence on the constant $\delta$ is also believed to be suboptimal, and the “correct” dependence shoould be $\delta^{-2}$. The question of improving this dependence is left as an open problem.

The use of a combination of random sign matrices and various forms of subsampled Fourier matrices was also used in the work of Ailon and Chazelle and later Ailon and Liberty , as well as that of Matousek . Here we obtain improved analysis using recent work by Rudelson and Vershynin for sparse reconstruction .

An underlying idea common to both random projections and sparse reconstruction is the preservation of metric information under a dimension reducing transformation. In sparse reconstruction theory, this property is known as restricted isometry . A matrix $\Phi$ is a restricted isometry with sparseness paramater $r$ if for some $\delta>0$,

In , Rudelson and Vershynin construct a distribution over $k\times n$ matrices $\Phi$ such that, with high probability, $\Phi$ has the restricted isometry property with sparseness parameter $r$ and arbitrarily small $\delta>0$.Their analysis is done over the complex field, but we restrict the discussion to the reals here. In their analysis, $k=O(\delta^{-2}r\log(n)\cdot\log^{2}(r)\log(r\log n))$ and $\Phi$ can be applied (to a given vector $x$) in running time $O(n\log n)$. Assuming $r$ polynomial in $n$, this takes the simpler form of $k=O(\delta^{-2}r\log^{4}n)$.In their work, the dependence of $k$ on $\delta$ is not analyzed because $\delta$ is assumed to be fixed (for sparse signal reconstruction purposes, this dependence is not important). It is not hard to derive the quadratic dependence of $k$ in $\delta^{-1}$ from their work. In fact, $\Phi$ is (up to a constant) nothing other than a random choice of $k$ rows from the (unnormalized) Hadamard matrix, defined as $\Psi_{\omega,t}=(-1)^{\langle\omega,t\rangle}$, where $\langle\cdot,\cdot\rangle$ is the dot product over the binary field, $n$ is assumed to be a power of $2$ and $\omega,t$ are thought of as $\log n$ dimensional vectors over the binary field in an obvious way.Rudelson and Vershynin use the complex Discrete Fourier Transform matrix, but their analysis does not change when using the Hadamard matrix. As a corollary of the result, one obtains a universal matrix for reconstructing sparse signals, which can be applied to a vector in time $O(n\log n)$. The conjecture is that the same distribution with $k=O(\delta^{-2}r\log n)$ should work as well, but this is a major open question beyond the scope of this work. For an excellent survey explaining how restricted isometry can be used for sparse reconstruction, and why designing such matrices with good computational properties is important we refer the readers to and to references therein.

with constant probability. Additionally, the number of steps required for applying $\Phi$ on any given $x$ is $O(n\log n)$. In their result $k$ was taken as $O(\delta^{-2}\log N)$, which is also essentially the best possible . Unfortunately, both results break down when $k=\Omega(n^{1/2})$.Ailon and Chazelle and Ailon and Liberty used $d$ to denote the data dimension, $n$ its cardinality and $\varepsilon$ the sought distortion bound. Here we follow Rudelson and Vershynin’s convention using $n$ to denote the dimension and $\delta$ the distortion bound. We now use $N$ to denote the data cardinality. Assuming the tolerance parameter $\delta$ fixed, this limitation can be rephrased as follows: The techniques fail when the number of vectors $N$ is in $\exp\{\Omega(n^{1/2})\}$.

In both Ailon and Chazelle and Ailon and Liberty’s results, as well as in previous work the bounds (1.2) are obtained by proving strong tail bounds on the distribution of the estimator $\|\Phi y\|_{2}$, and then applying a simple union bound on the finite collection $Y$. It is worth a moment’s thought to realize that Ailon and Chazelle’s result as well as that of Ailon and Liberty can be used for restricted isometry as well. Indeed, a simple epsilon-net argument for the set of $r$-sparse vectors can turn that set into a finite set of $\exp\{O(r\log n)\}$ vectors, on which a union bound can be applied. However, the current limitation of random projections mentioned above will limit $r$ to be in $n^{O(1/2-\mu)}$ (for arbitrarily small $\mu$). Interestingly, Rudelson and Vershynin’s result does not break down for $r$ polynomial in $n$. A careful inspection of their techniques reveals that instead of union bounding on a finite set of strongly concentrated random variables, they use a result due to Dudley to bound extreme values of Gaussian processes. Can this idea be used to improve and ? Intuitively there is no reason why a result which is designed for preserving the metric of sparse vectors should help with preserving the metric of any finite set of vectors. It turns out, luckily, that such a reduction can be done, though not in an immediate way. A suitable generalization of Rudelson and Vershynin’s result (Section 2), combined with Ailon and Chazelle and Ailon and Liberty’s method of random sign matrix preconditioning achieves this in Section 3.

2. Notation

Now let $\Phi$ be a random $k\times n$ matrix obtained by picking $k$ random rows from the unnormalized $n\times n$ Hadamard matrix (the Euclidean norm of each column of $\Phi$ is $\sqrt{k}$). Let $\Omega$ denote the probability space for the choice of $\Phi$.

Restricted isometry result generalization

We follow the main path of Rudelson et al. in to prove a more general formulation of their main theorem which is more suitable for us here.

[Derived from Rudelson and Vershynin] Let $\alpha>0$ be any real number. Define $E_{\alpha}$ as

In particular, if $\frac{(\log^{3/2}n)(\log^{1/2}k)}{\sqrt{k}}=O(\alpha)$, then

If we also assume that $k=\Theta(r\log^{4}n)$, then (2.3) will hold, from which we conclude that

Now we notice that $D_{y}=\frac{1}{\sqrt{r}}\operatorname{Id}_{\operatorname{supp}{y}}$, where for a set of indexes $T$ the diagonal matrix $\operatorname{Id}_{T}$ (as defined in ) has $1$ in diagonal position $i$ if and only if $i\in T$. Using this observation and multiplying (2.4) by $r$ we conclude that

which is exactly the main result of Rudelson and Vershynin in for restricted isometry.

The proof of Theorem 2.1 below points out the necessary changes to the proof of Theorem 3.4 in . The difference between the theorems is that in our case, the supremum in the definition of $E_{\alpha}$ is taken not only over the set of sparse vectors, but over a richer set. It turns out however that uses sparsity in a very limited way: In fact, the dominating effect of sparsity there is obtained using the fact that the $L_{1}$ norm of a sparse vector is small, compared to its $L_{2}$ norm. These arguments appear at the very end of their proof. For the sake of contributing to the self containment of the paper we walk through the main milestones of the proof of Theorem 3.4 in , and point out the changes necessary for our purposes. The reader is nevertheless encouraged to refer to the enlightening exposition in first.

Clearly $E[\frac{1}{k}D_{y}\Phi^{t}\Phi D_{y}]=D_{y}^{2}$. We define new independent random i.i.d. variables $\epsilon_{1},\dots,\epsilon_{n}$ obtaining each the values $+1,-1$ with equal probability. Let $\Pi$ denote the probability space for $\epsilon_{1},\dots,\epsilon_{n}$. It suffices to prove (using a symmetrization argument, see Lemma 6.3 in ) that

where $x_{i}$ is the (random) $i$’th row of $\Phi_{k}$. To that end, as claimed in (Lemma 3.5), if we can show that for any fixed choice of $\Phi$,

for some number $k_{1}$, then by taking $E_{\Omega}$ on both sides and using Jensen’s inequality (to swap $(\cdot)^{1/2}$ on the RHS with $E_{\Omega}$) and the triangle inequality, the conclusion would be that

Since $\|D_{y}^{2}\|=\|y\|^{2}_{\infty}\leq\alpha$, we would get the stated result. It thus suffices to prove (2.6) with $k_{1}=O((\log^{3/2}n)(\log^{1/2}k))$. To do so, continue by replacing the $k$ binary random variables $\epsilon_{1},\dots,\epsilon_{k}$ in (2.6) with $k$ Gaussian random variables $g_{1},\dots,g_{k}$ using a comparison principle (inequality (4.8) in ), reducing the problem to that of bounding the expected extreme value of a Gaussian process. Using Dudley’s inequality (Theorem 11.17 in ), as Rudelson and Vershynin do, one concludes that (2.6) will hold with $k_{1}$ taken as:

For a norm $\|\cdot\|_{\star}$, a set $S$ and number $u$, $\mathcal{N}(S,\|\cdot\|_{\star},u)$ denotes the minimal number of balls of radius $u$ in norm $\|\cdot\|_{\star}$ centered in points of $S$ needed to cover the set $S$,

$B$ is defined as $\cup_{y\in B_{2}\cap\alpha B_{\infty}}B_{y}$, where $B_{y}=\{D_{y}z:z\in B_{2}\}$, and

$\|x\|_{X}=\max_{i\leq k}|\langle x_{i},x\rangle|$, where we remind the reader that $x_{i}$ is the $i^{\prime}th$ row of $\Phi$.

Rudelson and Vershynin derive bounds on $\mathcal{N}(B_{RV},\|\cdot\|_{X},u)$ for small $u$ and for large $u$ separately, where in their case $B_{RV}$ was the set of $r$-sparse vectors of Euclidean norm $1$ (denoted by $D_{2}^{r,n}$ in ). The sparsity of the vectors in the set $B_{RV}$ is used in both derivations, as follows:

For large $u$, they use containment argument (11) in , asserting that $B_{RV}\subseteq\sqrt{r}B_{1}$. Note that by Cauchy Schwartz and the definition of $B$, $B\subseteq B_{1}$ hence we ”gain” a factor of $\sqrt{r}$ when deriving $k_{1}$.

For small $u$, inequality (13) in asserts that $\mathcal{N}(B_{RV},\|\cdot\|_{X},u)\leq d(n,r)(1+2/u)^{r}$, where $d(n,r)$ is the number of ways to choose $r$ elements from a set of $n$ elements. Since the best sparseness we can assume for vectors in $B$ here is trivially $n$, we replace the expression $d(n,r)$ with $d(n,n)=1$, and $(1+2/u)^{r}$ with $(1+2/u)^{n}$.To be exact, in they use the expression $(1+2K/u)^{r}$ and not $(1+2/u)^{r}$, but the parameter $K$ in their work can be taken as $1$ for our purposes.

Rudelson and Vershynin then derive a bound for $\int_{0}^{\infty}\mathcal{N}^{1/2}(B_{RV},\|\cdot\|_{X},u)du$ by balancing the two bounds at $u=1/\sqrt{r}$. In our case we balance at $u=1/\sqrt{n}$. The net result will lead to a $k_{1}$ which is as the one in the statement of Lemma 3.5 , except that the $\sqrt{r}$ will disappear and $\log r$ will be replaced by $\log n$. The conclusion is that we can take $k_{1}$ to be

Random Projections

Our main result claims that the same construction used by Rudelson et al. also gives improved bounds for random projections. In what follows, we fix $r$ to be $\lceil\delta^{-2}\log N\rceil$ and $\alpha$ to be $1/\sqrt{r}$. Additionally, we assume that $\Phi$ is such that

Indeed, Theorem 2.1 guarantees that this holds with probability at least $0.99$ in $\Omega$.

Let $Y\subseteq B_{2}$ denote a set of cardinality $N$, and let $\Phi$ satisfy (3.1). With probability at least $0.98$ (in $\Gamma$) we have the following uniform bound for all $y\in Y$:

Let $r$ and $\alpha$ be defined as in Section 2. For each $y\in Y$ we write $y={\hat{y}}+{\check{y}}$, where ${\hat{y}}$ is the restriction of $y$ to its $r$ largest (in absolute value) coordinates and ${\check{y}}$ is the restriction to its remaining coordinates. Note that $\|y\|^{2}=\|{\hat{y}}\|^{2}+\|{\check{y}}\|^{2}$ and that ${\hat{y}}$ is $r$-sparse and that $\|{\check{y}}\|_{\infty}\leq\alpha$.

For the first term we have $\left\|\frac{1}{\sqrt{k}}\Phi D_{{\hat{y}}}b\right\|^{2}=\|{\hat{y}}\|^{2}+O(\delta)$ from Theorem 2.1 and the fact that ${\hat{y}}$ is $r$-sparse.

In what follows we will use the bound on $\|{\check{y}}\|_{\infty}$ to show that with high probability, for all $y\in Y$, $\left\|\frac{1}{\sqrt{k}}\Phi D_{{\check{y}}}b\right\|^{2}=\|{\check{y}}\|^{2}+O(\delta)$. A similar argument will bound the cross product $\frac{2}{k}b^{t}D_{{\hat{y}}}\Phi^{t}\Phi D_{{\check{y}}}b$. Combining the three gives the desired result that $\left\|\frac{1}{\sqrt{k}}\Phi D_{y}b\right\|^{2}=\|y\|^{2}+O(\delta)$.

We start by analyzing the measure concentration properties of $\left\|\frac{1}{\sqrt{k}}\Phi D_{{\check{y}}}b\right\|^{2}$. Let $X_{{\check{y}}}$ be the Rademacher random variable defined by

Let $\mu_{\check{y}}$ denote a median of $X_{\check{y}}$. By Talagrand , we have that for all $t>0$,

for some global $C_{2}$, where $\sigma_{\check{y}}=\left\|\frac{1}{\sqrt{k}}\Phi D_{{\check{y}}}\right\|$. By the triangle inequality and Equation (3.1) we have $\sigma_{\check{y}}^{2}=\|\frac{1}{k}D_{\check{y}}\Phi^{t}\Phi D_{\check{y}}-D_{\check{y}}^{2}+D_{\check{y}}^{2}\|\leq\alpha^{2}+\|D_{\check{y}}^{2}\|$. Clearly $\|D_{\check{y}}\|=\|{\check{y}}\|_{\infty}\leq\alpha$. Hence, $\sigma_{\check{y}}^{2}=O(\alpha^{2})$. From the fact that $E[X_{\check{y}}^{2}]=\|{\check{y}}\|^{2}$ and using Appendix A and (3.2)-(3.3) We conclude that $\|{\check{y}}\|-O(\sigma_{{\check{y}}})\leq\mu_{\check{y}}\leq\|{\check{y}}\|+O(\sigma_{\check{y}})$. Hence, again using (3.2)-(3.3) and union bounding over the $N$ vectors in $Y$, we conclude that with probability $0.99$, uniformly for all $y\in Y$:

We now bound the cross term $Z=\frac{1}{k}b^{t}D_{{\hat{y}}}\Phi^{t}\Phi D_{{\check{y}}}b$ ($y$ is now held fixed). By disjointness of $\operatorname{supp}({\hat{y}})$ and $\operatorname{supp}({\hat{y}})$, $E[Z]=0$. Decompose $b$ into ${\check{b}}+{\hat{b}}$, where $\operatorname{supp}({\check{b}})=\operatorname{supp}({\check{y}})$ and $\operatorname{supp}({\hat{b}})=\operatorname{supp}({\hat{y}})$. For any fixed ${\hat{b}}$, the function $Z$ is linear (and hence convex) in ${\check{b}}$. Also for all possible values ${\hat{b}}^{\prime}$ of ${\hat{b}}$, $E[Z|{\hat{b}}={\hat{b}}^{\prime}]=0$. Hence, again by Talagrand,

where $\mu_{\hat{b}}^{\prime}$ is a median of $(Z|{\hat{b}}={\hat{b}}^{\prime})$, and $\sigma_{{\hat{b}}^{\prime}}=\|\frac{1}{k}({\hat{b}}^{\prime})^{t}D_{{\hat{y}}}\Phi^{t}\Phi D_{{\check{y}}}\|$. Clearly,

Again using Appendix A and $E[Z|{\hat{b}}={\hat{b}}^{\prime}]=0$ gives that $|\mu_{\hat{b}}^{\prime}|=O(\alpha)$, and again we conclude using a union bound that with probability at least $0.99$, uniformly for all $y\in Y$, $\left|\frac{1}{k}b^{t}D_{{\hat{y}}}\Phi^{t}\Phi D_{{\check{y}}}b\right|=O(\delta)$.

Tying it all together, we conclude that with probability at least $0.98$, uniformly for all $y\in Y$,

Conclusions

The obvious problems left open are those of (1) improving the dependence of $k$ in $\delta$ (from $\delta^{-4}$ to $\delta^{-2}$) and (2) removing the dependence of $k$ in $\operatorname{polylog}(n)$. Other directions of research include not only reducing the computational efficiency of random dimension reduction, but also the amount of randomness needed for the construction.

Acknowledgements

We thank Emmanuel Candes for helpful discussions.

References

Appendix A

For any real valued random variable $Z$ such that for all $t>0$

we have that $\sqrt{E(Z^{2})}-O(\sigma)\leq\mu\leq\sqrt{E(Z^{2})}+O(\sigma)$.

Define the variable $Z^{\prime}=(Z-\mu)/\sigma$.

Clearly, $E[Z^{\prime}]=O(1)$ gives $E(Z)=\mu+O(\sigma)$. In the same way we get $E[Z^{\prime 2}]=O(1)$. Thus, $E[Z^{2}]-2\mu E[Z]+\mu^{2}=O(\sigma^{2})$ and $E[Z^{2}]=(\mu\pm O(\sigma))^{2}$