New and improved Johnson-Lindenstrauss embeddings via the Restricted Isometry Property

Felix Krahmer, Rachel Ward

Introduction

The Johnson-Lindenstrauss (JL) Lemma states that any set of pp points in high dimensional Euclidean space can be embedded into O(ε−2log⁡(p))O(\varepsilon^{-2}\log(p)) dimensions, without distorting the distance between any two points by more than a factor between 1−ε1-\varepsilon and 1+ε1+\varepsilon. In its original form, the Johnson-Lindenstrauss Lemma reads as follows.

As shown in , the bound for the size of mm is tight up to an O(log⁡(1/ε))O(\log(1/\varepsilon)) factor. In the original paper of Johnson and Lindenstrauss, it was shown that a random orthogonal projection, suitably normalized, provides such an embedding with high probability . Later, this property was also verified for Gaussian random matrices, among other random matrix constructions . As a consequence, the JL Lemma has become a valuable tool for dimensionality reduction in a myriad of applications ranging from computer science , numerical linear algebra , manifold learning , and compressed sensing , , .

In most of these frameworks, the map ff under consideration is a linear map represented by an m×Nm\times N matrix Φ\Phi. In this case, one can consider the set of differences E={xi−xj}E=\{x_{i}-x_{j}\}; to prove the theorem, one then needs to show that

When Φ\Phi is a random matrix, the proof that Φ\Phi satisfies the JL lemma with high probability boils down to showing a concentration inequality of the type

In order to reduce storage space and implementation time of such embeddings, the design of structured random JL embeddings has been an active area of research in recent years ; see or for a good overview of these efforts. Of particular importance in this context is whether fast (i.e. O(Nlog⁡(N))O(N\log(N))) multiplication algorithms are available for the resulting matrices. Fast JL embeddings with optimal embedding dimension m=O(ε−2log⁡(p))m=O(\varepsilon^{-2}\log(p)) were first constructed by Ailon and Chazelle , but their embeddings are fast only for p≲eN1/3p\lesssim e^{N^{1/3}} vectors. This restriction on the number of vectors was later weakened to p≲eN1/2p\lesssim e^{N^{1/2}} . In , fast JL embeddings were constructed without any restrictions on the number of vectors, but the authors only provide sub-optimal embedding dimension m=O(ε−4log⁡(p)log⁡4(N))m=O(\varepsilon^{-4}\log(p)\log^{4}(N)). In this paper, we provide the first unrestricted fast JL construction with optimal embedding dimension up to logarithmic factors in NN. Note that in the range p≳eN1/2p\gtrsim e^{N^{1/2}} not covered by the constructions in , a logarithmic factor in NN is bounded by log⁡(log⁡(p))\log(\log(p)), and thus plays a minor role.

The restricted isometry constant δk\delta_{k} is defined as the smallest value of δ\delta for which (4) holds.

In particular, if Φ\Phi has (2k,δ2k)(2k,\delta_{2k})-RIP with δ2k≤2/(3+74)≈.4627\delta_{2k}\leq 2/(3+\sqrt{\frac{7}{4}})\approx.4627, and if y=Φxy=\Phi x admits a kk-sparse solution x#x^{\#}, then x#=arg⁡min⁡Φz=y∥z∥1x^{\#}=\arg\min_{\Phi z=y}\|z\|_{1} .

The similarity between the expressions in (2) and (4) suggests a connection between the JL lemma and the Restricted Isometry Property. A first result in this direction was established in , wherein it was shown that random matrices satisfying a concentration inequality of type (3) (and hence the JL Lemma) satisfy the RIP of optimal order. More precisely, the authors prove the following theorem.

Suppose that m,Nm,N, and 0<δ<10<\delta<1 are given. If the probability distribution generating the m×Nm\times N matrices Φ\Phi satisfies the concentration inequality (3) with ε=δ\varepsilon=\delta and absolute constant c0c_{0}, then there exist absolute constants c1,c2c_{1},c_{2} such that with probability ≥1−2e−c2δ2m,\geq 1-2e^{-c_{2}\delta^{2}m}, the RIP (4) holds for Φ\Phi with the prescribed δ\delta and any k≤c1δ2m/log⁡(N/k)k\leq c_{1}\delta^{2}m/\log(N/k).

In this sense, the JL Lemma implies the Restricted Isometry Property.

Contribution of this work.

We prove a converse result to Theorem 1.3: We show that RIP matrices, with randomized column signs, provide Johnson-Lindenstrauss embeddings that are optimal up to logarithmic factors in the ambient dimension. In particular, RIP matrices of optimal order provide Johnson-Lindenstrauss embeddings of optimal order as such, up to a logarithmic factor in NN (see Theorem 3.1). Note that without randomization, such a converse is impossible as vectors in the null space of the fixed parent matrix are always mapped to zero.

This observation has several consequences in the area of compressed sensing, and also allows us to obtain improved JL embedding results for several matrix constructions with existing RIP bounds . Of particular interest is the random partial Fourier or the random partial Hadamard matrix, which is formed by choosing a random subset of mm rows from the N×NN\times N discrete Fourier or Hadamard matrix respectively, and with high probability has (k,δ)(k,\delta)-RIP if the embedding dimension m≳δ−2klog⁡4(N)m\gtrsim\delta^{-2}k\log^{4}(N). For these matrices with randomized column signs, the running time for matrix-vector multiplication is O(Nlog⁡(N))O(N\log(N)) as opposed to the running time of O(Nm)O(Nm) for purely random matrices. For such constructions, the previous best-known embedding dimension to ensure that (2) holds with probability 1−η1-\eta, given by Ailon and Liberty , is m≍ε−4log⁡(p/η)log⁡4(N)m\asymp\varepsilon^{-4}\log(p/\eta)\log^{4}(N). We can improve their result to have optimal dependence on the distortion, ε\varepsilon, showing that m≍ε−2log⁡(p/η)log⁡4(N)m\asymp\varepsilon^{-2}\log(p/\eta)\log^{4}(N) rows suffice for the embedding.

This paper is structured as follows: Section 2 introduces necessary notation. In Section 3, we state our main results, and Section 4 gives concrete examples of how these results improve on the best-known JL bounds for several matrix constructions as well as applications of our findings in compressed sensing. In Section 5 we give the relevant concentration inequalities and explicit RIP-based matrix inequalities that are needed for the proofs, which are then carried out in Section 6.

Notation

The main results

Along the way, our method provides a direct converse to Theorem 1.3:

Concrete examples and applications

Using Theorem 3.1, we can improve on the best Johnson-Lindenstrauss bounds for several matrix constructions that are known to have the Restricted Isometry Property:

For measures with discrete support, such constructions are equivalent to choosing mm rows at random from an N×NN\times N matrix with orthonormal rows and uniformly bounded entries. Examples include the random partial Fourier matrix or random partial Hadamard matrix, formed from the discrete Fourier matrix or discrete Hadamard matrix respectively. (In the Fourier case, we distribute the resulting real and complex parts in different coordinates, inducing an additional factor of 22.) Note that the structure of these matrices allows for fast matrix vector multiplication. Recently, Ailon and Liberty verified the JL Lemma for such constructions, with column signs randomized, when m≳ε−4log⁡(p)log⁡4(N)m\gtrsim\varepsilon^{-4}\log(p)\log^{4}(N). Our result improves the factor of ε−4\varepsilon^{-4} in their result to the optimal dependence ε−2\varepsilon^{-2}. We note that while their proof also uses the RIP, it also requires arguments from that are specific to discrete bounded orthonormal systems.

Examples of bounded orthonormal systems connected to continuous measures include the trigonometric polynomials and Chebyshev polynomials, which are orthogonal with respect to the uniform and Chebyshev measures, respectively. The Legendre system, while not uniformly bounded, can still be transformed via preconditioning to a bounded orthonormal system with respect to the Chebyshev measure . Note that all of these constructions have an associated fast transform.

Partial circulant matrices.

Other classes of structured random matrices known to have the RIP include partial circulant matrices . In one such set-up, the first row of the N×NN\times N matrix is a Gaussian or Rademacher random vector, and each subsequent row is created by rotating one element to the right relative to the preceding row vector. Again, mm rows of this matrix are sampled, but in contrast to partial Fourier or Hadamard matrices, the selection need not be random. Using that convolution corresponds to multiplication in the Fourier domain, these matrices have associated fast matrix-vector multiplication routines. In , such matrices were shown to have the RIP with high probability for m≳max⁡(δ−1k32log⁡32(N),δ−2klog⁡4(N))m\gtrsim\operatorname{max}\left(\delta^{-1}k^{\frac{3}{2}}\log^{\frac{3}{2}}(N),\delta^{-2}k\log^{4}(N)\right).

On the other hand, such a matrix composed with a diagonal matrix of random signs was shown to be a JL embedding with high probability as long as m≳ε−2log⁡2(p)m\gtrsim\varepsilon^{-2}\log^{2}(p) . Through Theorem 3.1, the same results also obtain if m\gtrsim\operatorname{max}\Big{(}\varepsilon^{-1}\log^{3/2}\left(\frac{4p}{\eta}\right)\log^{\frac{3}{2}}(N),\varepsilon^{-2}\log\left(\frac{4p}{\eta}\right)\log^{4}(N)\Big{)}. For large pp, this is an improvement compared to .

Deterministic constructions.

Several deterministic constructions of RIP matrices are known, including a recent result in that requires only m≳k2−μm\gtrsim k^{2-\mu}. We refer the reader to the exposition in for a good overview in this direction; we highlight two such deterministic constructions here. Using finite fields, DeVore provides deterministic constructs of cyclic -11-valued matrices with (k,δ)(k,\delta)-RIP with m≳δ−2k2log⁡2(N)m\gtrsim\delta^{-2}k^{2}\log^{2}(N). Iwen provides deterministic constructions of -11-valued matrices whose number theoretic properties allow their products with Discrete Fourier Transform (DFT) matrices to be well approximated using a few highly sparse matrix multiplications. Both the binary-valued matrices and their products with the DFT yield (k,δ)(k,\delta)-RIP matrices with m≳δ−2k2log⁡2(N)m\gtrsim\delta^{-2}k^{2}\log^{2}(N). By Theorem 3.1, the class of matrices that results by randomizing the column signs of either of these deterministic constructions satisfies the JL Lemma with m≳ε−2log⁡2(p)log⁡2(N)m\gtrsim\varepsilon^{-2}\log^{2}(p)\log^{2}(N).

Note that the amount of randomness needed to construct such embeddings is still comparable to the first two examples, requiring NN random bits. Under the model assumption that the entries of each vector x∈Ex\in E to be embedded has random signs, however, the required randomness in the matrix is removed completely.

For each of the aforementioned examples, we summarize the number of dimensions mm that are known to be sufficient (k,δ)(k,\delta)-RIP to hold. We also list the previously best known bound for JL embedding dimension (if there is one) along with the JL bounds obtained from Theorem 3.1. Where Theorem 3.1 yields a better bound than previously known, at least for some range of parameters, we highlight the result in bold face. In each of the bounds, we list only the dependence on δ,k\delta,k, and NN, or ε,k,\varepsilon,k, and NN, omitting absolute constants.

Compressed sensing in redundant dictionaries.

As shown recently in , concentration inequalities of type (3) allow for the extension of the compressed sensing methodology to redundant dictionaries – in particular, tight frames – as opposed to orthonormal bases only. Since signals with sparse representations in redundant dictionaries comprise a much more realistic model of nature, this extension of compressed sensing is fundamental. Our results show that basically all random matrix constructions arising in the standard theory of compressed sensing (i.e., based on RIP estimates) also yield compressed sensing matrices for the redundant framework.

Compressed sensing with cross validation.

Compressed sensing algorithms are designed to recover approximately sparse signals; if this assumption is violated, they may yield solutions far from the input signal. In , a method of cross validation is introduced to detect such situations, and to obtain tight bounds on the error incurred by compressed sensing reconstruction algorithms in general. There, a subset y1=Φ1xy_{1}=\Phi_{1}x of the mm measurements y=Φxy=\Phi x are held out from the reconstruction algorithm and only the remaining measurements y2=Φ2xy_{2}=\Phi_{2}x are used to produce a candidate approximation x^\widehat{x} to the unknown xx. If the hold-out matrix Φ1\Phi_{1} satisfies the Johnson-Lindenstrauss Lemma, then the observable quantity ∥Φ1(x−x^)∥2\|\Phi_{1}(x-\widehat{x})\|_{2} can be used as a reliable proxy for the unknown error ∥x−x^∥2\|x-\widehat{x}\|_{2}. Our work shows that any RIP matrix as in the standard compressed sensing framework can be used for cross validation up to a randomization of its column signs.

Optimal asymptotics in δ𝛿\delta for RIP to hold.

As mentioned above, it can be shown using a Gelfand width argument that m≍klog⁡(Nk)m\asymp k\log(\frac{N}{k}) is the optimal asymptotics (in NN and kk) of the embedding dimension for a matrix with the restricted isometry property (4). Our results – combined with the known optimality of the asymptotics m=ε−2log⁡(p)m=\varepsilon^{-2}\log(p) for the embedding dimension in the Johnson-Lindenstrauss Lemma (1.1) – imply that up to a factor of log⁡(1δ)\log\left(\frac{1}{\delta}\right), m≍δ−2m\asymp\delta^{-2} is the optimal asymptotics in the restricted isometry constant δ\delta for fixed NN and kk as δ→0\delta\rightarrow 0. Recall that this rate is realized by many of the above examples, such as Gaussian random matrices.

Proof Ingredients

The proof of Theorem 3.1 relies on concentration inequalities for Rademacher sequences and explicit RIP-based norm estimates. The first concentration result is a classical inequality by Hoeffding .

The second concentration of measure result is a deviation bound for Rademacher chaos. There are many such bounds in the literature; the following inequality dates back to , but appeared with explicit constants and with a much simplified proof as Theorem 1717 in .

Let XX be the N×NN\times N matrix with entries xi,jx_{i,j} and assume that xi,i=0x_{i,i}=0 for all i∈[N]i\in[N]. Let ξ=(ξj)j=1N\xi=(\xi_{j})_{j=1}^{N} be a Rademacher sequence. Then, for any t>0t>0,

We also need the following basic estimate for RIP matrices (see for instance Proposition 2.52.5 in ).

The proof of our norm estimate for RIP-matrices uses Proposition 5.3, and relies on the observation commonly used in the theory of compressed sensing (see for example ) that for zz in decreasing arrangement and ∥z∥2=1\|z\|_{2}=1, for J≥2J\geq 2 one has ∥z(J)∥∞≤1s∥z(J−1)∥2\|z_{(J)}\|_{\infty}\leq\frac{1}{\sqrt{s}}\|z_{(J-1)}\|_{2} and thus ∥z(♭)∥∞≤1/s\|z_{(\flat)}\|_{\infty}\leq 1/\sqrt{s}.

The following bounds hold: ∥C∥≤δs,∥C∥F≤δs,\quad\|C\|\leq\frac{\delta}{s},\quad\|C\|_{\cal{F}}\leq\frac{\delta}{\sqrt{s}},\quad and ∥v∥2≤δs\quad\|v\|_{2}\leq\frac{\delta}{\sqrt{s}}.

To obtain (10), we use the inequality of arithmetic and geometric means; to obtain (9), we use Proposition 5.3.

Proof of the main results

We begin by proving Theorem 3.1. Without loss of generality, we assume that all x∈Ex\in E are normalized so that ∥x∥2=1\|x\|_{2}=1. Furthermore, assume that k=2sk=2s is even.

We first consider a fixed x∈Ex\in E, eventually taking a union bound over all xx. We further assume that xx is in decreasing arrangement. To achieve this, we reorder the entries of xx, and permute the columns of Φ\Phi accordingly. This has no impact on the following estimates, as the Restricted Isometry Property of the matrix Φ\Phi is invariant under permutations of its columns. We need to estimate

As Φ\Phi has the Restricted Isometry Property of order k≥sk\geq s and level δ\delta, it also has the RIP of order ss and level δ\delta, and each Φ(J)\Phi_{(J)} is almost an isometry. Hence, noting that ∥Dx(J)ξ(J)∥2=∥Dξ(J)x(J)∥2=∥x(J)∥2\|D_{x_{(J)}}\xi_{(J)}\|_{2}=\|D_{\xi_{(J)}}x_{(J)}\|_{2}=\|x_{(J)}\|_{2}, the first term can be estimated as follows.

Thus, using that δ≤ε/4,\delta\leq\varepsilon/4,

To estimate the second term, fix ξ(1)=:b\xi_{(1)}=:b and consider the random variable

with vv as in Proposition 5.4. By Hoeffding’s inequality (Proposition 5.1) combined with Proposition 5.4,

In order for this probability to be less than η/2\eta/2, we need:

In order for this probability to be less than η/2\eta/2, we need:

By assumption, δ≤ε4\delta\leq\frac{\varepsilon}{4}, so conditions (13) and (15) are satisfied by setting τ=.55,γ=.1\tau=.55,\gamma=.1, and s≥20log⁡(4p/η)s\geq 20\log{(4p/\eta)} (that is, k=2s≥40log⁡(4p/η)k=2s\geq 40\log{(4p/\eta)}). Then the second term is bounded by .2δ.2\delta in absolute value, and the last term is bounded by .55δ.55\delta. Together with the deterministic RIP-based estimate for the first term, this implies the Theorem. ∎

Remarks:

As shown in , a random matrix Φ\Phi whose entries follow a subgaussian distribution is known to have with high probability the Restricted Isometry Property of best possible order, that is, one can choose m≍δ−2klog⁡(Nk).m\asymp\delta^{-2}k\log\left(\frac{N}{k}\right). When k≥40log⁡(4pη)k\geq 40\log\left(\frac{4p}{\eta}\right), Φ\Phi is a JL embedding by Theorem 3.1, and our resulting bound for mm is optimal up to a single logarithmic factor in NN. This shows that Theorem 3.1 must also be optimal up to a single logarithmic factor in NN.

Acknowledgments

The authors would like to thank Holger Rauhut, Deanna Needell, Jan Vybíral, Mark Tygert, Mark Iwen, Justin Romberg, Mark Davenport, and Arie Israel for valuable discussions on this topic. Rachel Ward gratefully acknowledges the partial support of National Science Foundation Postdoctoral Research Fellowship. Felix Krahmer gratefully acknowledges the partial support of the Hausdorff Center for Mathematics. Finally, both authors are grateful for the support of the Institute of Advanced Study through the Park City Math Institute where this project was initiated.

References