Combining geometry and combinatorics: A unified approach to sparse signal recovery

Introduction

With the rise in high-speed data transmission and the exponential increase in data storage, it is imperative that we develop effective data compression techniques, techniques which accomodate both the volume and speed of data streams. A new approach to compressing $n$-dimensional vectors (or signals) begins with linear observations or measurements. For a signal $x$, its compressed representation is equal to $\bm{\Phi}x$, where $\bm{\Phi}$ is a carefully chosen $m\times n$ matrix, $m\ll n$, often chosen at random from some distribution. We call the vector $\bm{\Phi}x$ the measurement vector or a sketch of $x$. Although the dimension of $\bm{\Phi}x$ is much smaller than that of $x$, it retains many of the essential properties of $x$.

There are several reasons why linear compression or sketching is of interest. First, we can easily maintain a linear sketch $\bm{\Phi}x$ under linear updates to the signal $x$. For example, after incrementing the $i$-th coordinate $x_{i}$, we simply update the sketch as $\bm{\Phi}(x+e_{i})=\bm{\Phi}x+\bm{\Phi}e_{i}$. Similarly, we also easily obtain a sketch of a sum of two signals given the sketches for individual signals $x$ and $y$, since $\bm{\Phi}(x+y)=\bm{\Phi}x+\bm{\Phi}y$. Both properties are very useful in several computational areas, notably computing over data streams [AMS99, Mut03, Ind07], network measurement [EV03], query optimization and answering in databases [AMS99].

Another scenario where linear compression is of key importance is compressed sensing [CRT06, Don06a], a rapidly developing area in digital signal processing. In this setting, $x$ is a physical signal one wishes to sense (e.g., an image obtained from a digital camera) and the linearity of the observations stems from a physical observation process. Rather than first observing a signal in its entirety and then compressing it, it may be less costly to sense the compressed version directly via a physical process. A camera “senses” the vector by computing a dot product with a number of pre-specified measurement vectors. See [TLW+06, DDT+08] for a prototype camera built using this framework. Other applications of linear sketching include database privacy [DMT07].

There are many algorithms for recovering sparse approximations (or their variants) of signals from their sketches. The early work on this topic includes the algebraic approach of [Man92](cf. [GGI+02a]). Most of the known algorithms, however, can be roughly classified as either combinatorial or geometric.

Combinatorial approach. In the combinatorial approach, the measurement matrix $\bm{\Phi}$ is sparse and often binary. Typically, it is obtained from an adjacency matrix of a sparse bipartite random graph. The recovery algorithm proceeds by iteratively, identifying and eliminating ‘‘large’’ coefficients In the non-sketching world, such methods algorithms are often called “weak greedy algorithms”, and have been studied thoroughly by Temlyakov [Tem02] of the vector $x$. The identification uses non-adaptive binary search techniques. Examples of combinatorial sketching and recovery algorithms include [GGI+02b, CCFC02, CM04, GKMS03, DWB05, SBB06b, SBB06a, CM06, GSTV06, GSTV07, Ind08, XH07] and others.

The typical advantages of the combinatorial approach include fast recovery, often sub-linear in the signal length $n$ if $k\ll n$ and fast and incremental (under coordinate updates) computation of the sketch vector $\bm{\Phi}x$. In addition, it is possible to construct efficient (albeit suboptimal) measurement matrices explicitly, at least for simple type of signals. For example, it is known [Ind08, XH07] how to explicitly construct matrices with $k2^{(\log\log n)^{O(1)}}$ measurements, for signals $x$ that are exactly $k$-sparse. The main disadvantage of the approach is the suboptimal sketch length.

Geometric approach. This approach was first proposed in the papers [CRT06, Don06a] and has been extensively investigated since then (see [Gro06] for a bibliography). In this setting, the matrix $\bm{\Phi}$ is dense, with at least a constant fraction of non-zero entries. Typically, each row of the matrix is independently selected from a sub-exponential $n$-dimensional distribution, such as Gaussian or Bernoulli. The key property of the matrix $\bm{\Phi}$ which yields efficient recovery algorithms is the Restricted Isometry Property [CRT06], which requires that for any $k$-sparse vector $x$ we have $\|\bm{\Phi}x\|_{2}=(1\pm\delta)\|x\|_{2}$. If a matrix $\bm{\Phi}$ satisfies this property, then the recovery process can be accomplished by finding a vector $x_{*}$ using the following linear program:

The advantages of the geometric approach include a small number of measurements ($O(k\log(n/2k))$ for Gaussian matrices and $O(k\log^{O(1)}n)$ for Fourier matrices) and resiliency to measurement errors Historically, the geometric approach resulted also in the first deterministic or uniform recovery algorithms, where a fixed matrix $\bm{\Phi}$ was guaranteed to work for all signals $x$. In contrast, the early combinatorial sketching algorithms only guaranteed $1-1/n$ probability of correctness for each signal $x$. However, the papers [GSTV06, GSTV07] showed that combinatorial algorithms can achieve deterministic or uniform guarantees as well.. The main disadvantage is the running time of the recovery procedure, which involves solving a linear program with $n$ variables and $n+m$ constraints. The computation of the sketch $\bm{\Phi}x$ can be done efficiently for some matrices (e.g., Fourier); however, an efficient sketch update is not possible. In addition, the problem of finding an explicit construction of efficient matrices satisfying the RIP property is open [Tao07]; the best known explicit construction [DeV07] yields $\Omega(k^{2})$ measurements.

Connections. There has been some recent progress in obtaining the advantages of both approaches by decoupling the algorithmic and combinatorial aspects of the problem. Specifically, the papers [NV07, DM08, NT08] show that one can use greedy methods for data compressed using dense matrices satisfying the RIP property. Similarly [GLR08], using the results of [KT07], show that sketches from (somewhat) sparse matrices can be recovered using linear programming.

In addition, in Figure 2 we present very recent results discovered during the course of our research. Some of the running times of the algorithms depend on the “precision parameter” $D$, which is always bounded from the above by $\|x\|_{2}$ if the coordinates of $x$ are integers.

In this paper we give a sequence of results which indicate that the combinatorial and geometric approaches are, in a rigorous sense, different manifestations of a common underlying phenomenon. This enables us to achieve a unifying perspective on both approaches, as well as obtaining several new concrete algorithmic results.

Consider any $m\times n$ matrix $\bm{\Phi}$ that is the adjacency matrix of an $(k,\epsilon)$-unbalanced expander $G=(A,B,E)$, $|A|=n$, $|B|=m$, with left degree $d$, such that $1/\epsilon,d$ are smaller than $n$. Then the scaled matrix $\bm{\Phi}/d^{1/p}$ satisfies the ${\rm RIP}_{{p},{k},{\delta}}$ property, for $1\leq p\leq 1+1/\log n$ and $\delta=C\epsilon$ for some absolute constant $C>1$.

The fact that the unbalanced expanders yield matrices with RIP-$p$ property is not an accident. In particular, we show in Section 2 that any binary matrix $\bm{\Phi}$ in which each column has $d$ ones In fact, the latter assumption can be removed without loss of generality. The reason is that, from the RIP-$1$ property alone, it follows that each column must have roughly the same number of ones. The slight unbalance in the number of ones does not affect our results much; however, it does complicate the notation somewhat. As a result, we decided to keep this assumption throughout the paper. and which satisfies RIP-$1$ property with proper parameters, must be an adjacency matrix of a good unbalanced expander. That is, an RIP-$p$ matrix and the adjacency matrix of an unbalanced expander are essentially equivalent. Therefore, RIP-$1$ provides an interesting “analytic” formulation of expansion for unbalanced graphs. Also, without significantly improved explicit constructions of unbalanced expanders with parameters that match the probabilistic bounds (a longstanding open problem), we do not expect significant improvements in the explicit constructions of RIP-$1$ matrices.

Consider any $m\times n$ binary matrix $\bm{\bm{\Phi}}$ such that each column has exactly $d$ ones. If for some scaling factor $S>0$ the matrix $S\bm{\bm{\Phi}}$ satisfies the ${\rm RIP}_{{1},{s},{\delta}}$ property, then the matrix $\bm{\bm{\Phi}}$ is an adjacency matrix of an $(s,\epsilon)$-unbalanced expander, for

In the next step in Section 3, we show that the RIP-1 property of a binary matrix (or, equivalently, the expansion property) alone suffices to guarantee that the linear program $P_{1}$ recovers a good sparse approximation. In particular, we show the following

Let $\bm{\Phi}$ be an $m\times n$ matrix of an unbalanced $(2k,\epsilon)$-expander. Let $\alpha(\epsilon)=(2\epsilon)/(1-2\epsilon)$. Consider any two vectors $x,x_{*}$, such that $\bm{\Phi}x=\bm{\Phi}x_{*}$, and $\|x_{*}\|_{1}\leq\|x\|_{1}$. Then

where $x_{k}$ is the optimal $k$-term representation for $x$.

We also provide a noise-resilient version of the theorem; see Section 3 for details.

By combining Theorem 3 with the best known probabilistic constructions of expanders (Section 2) we obtain a scheme for sparse approximation recovery with parameters as in Figure 1. The scheme achieves the best known bounds for several parameters: the scheme is deterministic (one matrix works for all vectors $x$), the number of measurements is $O(k\log(n/k))$, the update time is $O(\log(n/k))$ and the encoding time is $O(n\log(n/k))$. In particular, this provides the first known scheme which achieves the best known measurement and encoding time bounds simultaneously. In contrast, the Gaussian and Fourier matrices are known to achieve only one optimal bound at a time. The fast encoding time also speeds-up the LP decoding, given that the linear program is typically solved using the interior-point method, which repeatedly performs matrix-vector multiplications. In addition to theoretical guarantees, random sparse matrices offer an attractive empirical performance. We show in Section 4 that the empirical behavior of binary sparse matrices with LP decoding is consistent with the analytic performance of Gaussian random matrices.

In the final part of the paper, we show that adjacency matrices of unbalanced expanders can be augmented to facilitate sub-linear time combinatorial recovery. This fact has been implicit in the earlier work [GSTV07, Ind08]; here we verify that indeed the expansion property is the sufficient condition guaranteeing correctness of those algorithms. As a result, we obtain an explicit construction of matrices with $O(k2^{(\log\log n)^{O(1)}})$ rows that are amenable to a sublinear decoding algorithm for all vectors (similar to that in [GSTV07]). Previous explicit constructions for sublinear time algorithms either had $\Omega(k^{2})$ rows [CM06] or had $O(k2^{(\log\log n)^{O(1)}})$ rows [Ind08, XH07] but were restricted to $k$-sparse signals or their slight generalizations. An additional (and somewhat unexpected) consequence is that the algorithm of [Ind08] is simple, effectively mimicking the well-known “parallel bit-flip” algorithm for decoding low-density parity-check codes.

where $x_{k}$ is the optimal $k$-term representation for $x$.

Figure 3 summarizes the connections among all of our results. We show the relationship between the combinatorial and geometric approaches to sparse signal recovery

Unbalanced expanders and RIP matrices

In this section we show that RIP-$p$ matrices for $p\approx 1$ can be constructed using high-quality expanders. The formal definition of the latter is as follows.

A $(k,\epsilon)$-unbalanced expander is a bipartite simple graph $G=(A,B,E)$ with left degree $d$ such that for any $X\subset A$ with $|X|\leq k$, the set of neighbors $N(X)$ of $X$ has size $|N(X)|\geq(1-\epsilon)d|X|$.

In constructing such graphs, our goal is to make $|B|$, $d$, and $\epsilon$ as small as possible, while making $k$ as close to $|B|$ as possible.

The following well-known proposition can be shown using the probabilistic method.

For any $n/2\geq k\geq 1$, $\epsilon>0$, there exists a $(k,\epsilon)$-unbalanced expander with left degree $d=O(\log(n/k)/\epsilon$ and right set size $O(kd/\epsilon)=O(k\log(n/k)/\epsilon^{2}$).

For any $n\geq k\geq 1$ and $\epsilon>0$, one can explicitly construct a $(k,\epsilon)$-unbalanced expander with left degree $d=2^{O(\log(\log(n)/\epsilon)))^{3}}$, left set size $n$ and right set size $m=kd/\epsilon^{O(1)}$.

The construction is given in [CRVW02], Theorem 7.3. Note that the theorem refers to notion of lossless conductors, which is equivalent to unbalanced expanders, modulo representing all relevant parameters (set sizes, degree, etc.) in the log-scale. After an additional $O(nd)$-time postprocessing, we can ensure that the graph is simple; i.e., it contains no duplicate edges. ∎

2. Construction of RIP matrices

An $m\times n$ matrix $\bm{\Phi}$ is said to satisfy ${\rm RIP}_{{p},{k},{\delta}}$ if, for any $k$-sparse vector $x$, we have

Observe that the definitions of ${\rm RIP}_{{1},{k},{\delta}}$ and ${\rm RIP}_{{2},{k},{\delta}}$ matrices are incomparable. In what follows below, we present sparse binary matrices with $O(k\log(n/k))$ rows that are RIP1,k,δ; it has been shown recently [Cha08] that sparse binary matrices cannot be RIP2,k,δ unless the number of rows is $\Omega(k^{2})$. In the other direction, consider an appropriately scaled random Gaussian matrix $G$ of $R\approx k\log(n)$ rows. Such a matrix is known to be RIP2,k,δ. To see that this matrix is not RIP1,k,δ, consider the signal $x$ consisting of all zeros except a single 1 and the signal $y$ consisting of all zeros except $k$ terms with coefficient $1/k$. Then $\|x\|_{1}=\|y\|_{1}$ but $\|Gx\|_{1}\approx\sqrt{k}\|Gy\|_{1}$.

Theorem 1 Consider any $m\times n$ matrix $\bm{\Phi}$ that is the adjacency matrix of an $(k,\epsilon)$-unbalanced expander $G=(A,B,E)$ with left degree $d$, such that $1/\epsilon,d$ are smaller than $n$. Then the scaled matrix $\bm{\Phi}/d^{1/p}$ satisfies the ${\rm RIP}_{{p},{k},{\delta}}$ property, for $1\leq p\leq 1+1/\log n$ and $\delta=C\epsilon$ for some absolute constant $C>1$.

The proof proceeds in two stages. In the first part, we show that the theorem holds for the case of $p=1$. In the second part, we extend the theorem to the case where $p$ is slightly larger than $1$.

The case of $p=1$. We order the edges $e_{t}=(i_{t},j_{t})$, $t=1\ldots dn$ of $G$ in a lexicographic manner. It is helpful to imagine that the edges $e_{1},e_{2}\ldots$ of $E$ are being added to the (initially empty) graph. An edge $e_{t}=(i_{t},j_{t})$ causes a collision if there exists an earlier edge $e_{s}=(i_{s},j_{s}),s<t$, such that $j_{t}=j_{s}$. We define $E^{\prime}$ to be the set of edges which do not cause collisions, and $E^{\prime\prime}=E-E^{\prime}$.

To upper bound the latter quantity, observe that the vectors satisfy the following constraints:

The coordinates of $z$ are monotonically non-increasing.

For each prefix set $P_{i}=\{1\ldots di\}$, $i\leq k$, we have $\|r_{|P_{i}}\|_{1}\leq\epsilon di$ - this follows from the expansion properties of the graph $G$.

$r_{|P_{1}}=0$, since the graph is simple.

It is now immediate that for any $r,z$ satisfying the above constraints, we have $r\cdot z\leq\|z\|_{1}\epsilon$. Since $\|z\|_{1}=d\|x\|_{1}$, the lemma follows. ∎

Lemma 9 immediately implies that $\left\|{\bm{\Phi}x}\right\|_{1}\geq d\left\|{x}\right\|_{1}(1-2\epsilon)$. Since for any $x$ we have $\left\|{\bm{\Phi}x}\right\|_{1}\leq d\left\|{x}\right\|_{1}$, it follows that $\bm{\Phi}/d$ satisfies the ${\rm RIP}_{{1},{k},{2\epsilon}}$ property.

The case of $p\leq 1+1/\log n$. Let $u=\bm{\Phi}x$. We will show that if $\epsilon$ is small enough, then the value of $\left\|{u}\right\|_{p}^{p}$ is close to $d\left\|{x}\right\|_{p}^{p}$.

We start from a few useful technical claims.

For any $j\in B$, define $u^{\prime}_{j}=x_{i}$ if $(i,j)\in E^{\prime}$, and $u^{\prime}_{j}=0$ otherwise. Also, define $u_{j}^{\prime\prime}=u_{j}-u^{\prime}_{j}$.

We have $\sum_{j}|u_{j}^{\prime\prime}|^{p}=\Theta(\epsilon)d\left\|{x}\right\|_{p}^{p}$.

By Lemma 9 we know that $\sum_{j}|u^{\prime}_{j}|=\sum_{(i,j)\in E^{\prime}}|x_{i}|\geq(1-\epsilon)d\left\|{x}\right\|_{1}$. Therefore

We have $\left\|{u}\right\|_{p}^{p}\geq(1-\Theta(\epsilon))d\left\|{x}\right\|_{p}^{p}$.

Define the set $S=\{j:\epsilon|u_{j}^{\prime}|\leq|u_{j}^{\prime\prime}|\}$. We have

It suffices to bound $\left\|{u^{\prime}}\right\|_{p}^{p}$ from below.

where we used Lemma 9 in the third line. Altogether

We have $\left\|{u}\right\|_{p}^{p}\leq(1+\Theta(\epsilon))d\left\|{x}\right\|_{p}^{p}$.

Define the set $T=\{j:|u_{j}^{\prime\prime}|\leq|u_{j}^{\prime}|\}$. Decompose $\left\|{u}\right\|_{p}^{p}$ into a sum

which, as in the proof of Lemma 13, can be bounded from above by

Combining Lemma 14 and Lemma 13 yields the proof of the theorem. ∎

The assumption that each column has exactly $d$ ones is not crucial, since the RIP-$1$ property itself implies that the number of ones in each column can differ by at most factor of $1+\delta$. All proofs in this paper are resilient to this slight unbalance. However, we decided to keep this assumption for the ease of notation.

Theorem 2 Consider any $m\times n$ binary matrix $\bm{\Phi}$ such that each column has exactly $d$ ones. If for some scaling factor $S>0$ the matrix $S\bm{\Phi}$ satisfies the ${\rm RIP}_{{1},{s},{\delta}}$ property, then the matrix $\bm{\Phi}$ is an adjacency matrix of an $(s,\epsilon)$-unbalanced expander, for

Note that, for small values of $\delta>0$, we have $\Big(1-\frac{1}{1+\delta}\Big)/(2-\sqrt{2})\approx\delta/(2-\sqrt{2})$.

Let $G=(A,B,E)$ be the graph with adjacency matrix $\bm{\Phi}$. Assume that there exists $X\subset A$, $|X|=k^{\prime}\leq k$ such that $|N(X)|<dk^{\prime}(1-\epsilon)$. We will construct two $n$-dimensional vectors $y,z$ such that $\left\|{y}\right\|_{1}=\left\|{z}\right\|_{1}=k^{\prime}$, but $\left\|{\bm{\Phi}z}\right\|_{1}/\left\|{\bm{\Phi}y}\right\|_{1}\leq 1-\epsilon(2-\sqrt{2})$, which is a contradiction.

The vector $y$ is simply the characteristic vector of the set $X$. Clearly, we have $\left\|{y}\right\|_{1}=k^{\prime}$ and $\left\|{\bm{\Phi}y}\right\|_{1}=dk$.

The vector $z$ is defined via a random process. For $i\in X$, define $r_{i}$ to be i.i.d. random variables uniformly distributed over $\{-1,1\}$. We define $z_{i}=r_{i}$ if $i\in X$, and $z_{i}=0$ otherwise. Note that $\left\|{z}\right\|_{1}=\left\|{y}\right\|_{1}=k^{\prime}$.

Let $C\subset N(X)$ be the “collision set”, i.e., the set of all $j\in N(X)$ such that the number $u_{j}$ of the edges from $j$ to $X$ is at least $2$. Let $|C|=l$. By the definition of the set $C$ we have $\sum_{j}u_{j}\geq 2l$. Moreover, from the assumption it follows that $\sum_{j}u_{j}\geq 2\epsilon dk^{\prime}$.

Let $v=\bm{\Phi}z$. We split $v$ into $v_{C}$ and $v_{C^{c}}$. Clearly, $\left\|{v_{C^{c}}}\right\|_{1}=k^{\prime}d-\sum_{j}u_{j}$. It suffices to show that $\left\|{v_{C}}\right\|_{1}$ is significantly smaller than $\sum_{j}u_{j}$ for some $z$.

The expected value of $\left\|{v_{C}}\right\|_{2}^{2}$ is equal to $\sum_{j}u_{j}$.

For each $j\in C$, the coordinate $v_{j}$ is a sum of $u_{j}$ independent random variables uniformly distributed over $\{-1,1\}$. The claim follows by elementary analysis. ∎

By Claim 15 we know that there exists $z$ such that $\left\|{v_{C}}\right\|_{2}\leq\sqrt{\sum_{j}u_{j}}\leq\frac{\sum_{j}u_{j}}{\sqrt{2l}}$. This implies that $\left\|{v_{C}}\right\|_{1}\leq\sqrt{l}\left\|{v_{C}}\right\|_{2}\leq\frac{\sum_{j}u_{j}}{\sqrt{2}}$. Therefore

LP decoding

In this section we show that if $A$ is an adjacency matrix of an expander graph, then the LP decoding procedure can be used for recovering sparse approximations.

Let $\Phi$ be an $m\times n$ adjacency matrix of an unbalanced $(2k,\epsilon)$-expander $G$ with left degree $d$. Let $\alpha(\epsilon)=(2\epsilon)/(1-2\epsilon)$. We also define $E(X:Y)=E\cap(X\times Y)$ to be the set of edges between the sets $X$ and $Y$.

Without loss of generality, we can assume that $S$ consists of the largest (in magnitude) coefficients of $y$. We partition coordinates into sets $S_{0},S_{1},S_{2},\ldots S_{t}$, such that (i) the coordinates in the set $S_{l}$ are not larger (in magnitude) than the coordinates in the set $S_{l-1}$, $l\geq 1$, and (ii) all sets but $S_{t}$ have size $k$. Therefore, $S_{0}=S$. Let $\bm{\Phi}^{\prime}$ be a submatrix of $\bm{\Phi}$ containing rows from $N(S)$.

From the equivalence of expansion and RIP-1 property we know that $\|\bm{\Phi}^{\prime}y_{S}\|_{1}=\|\bm{\Phi}y_{S}\|_{1}\geq d(1-2\epsilon)\|y_{S}\|_{1}$. At the same time, we know that $\|\bm{\Phi}^{\prime}y\|_{1}=0$. Therefore

From the expansion properties of $G$ it follows that, for $l\geq 1$, we have $|N(S\cup S_{l})|\geq d(1-\epsilon)|S\cup S_{l}|$. It follows that at most $d\epsilon 2k$ edges can cross from $S_{l}$ to $N(S)$, and therefore

It follows that $d(1-2\epsilon)\|y_{S}\|_{1}\leq 2d\epsilon\|y\|_{1}$, and thus $\|y_{S}\|_{1}\leq(2\epsilon)/(1-2\epsilon)\|y\|_{1}$. ∎

2. LP recovery

The following theorem provides recovery guarantees for the program $P_{1}$, by setting $u=x$ and $v=x_{*}$.

Theorem 3 Consider any two vectors $u,v$, such that for $y=v-u$ we have $\bm{\Phi}y=0$, and $\|v\|_{1}\leq\|u\|_{1}$. Let $S$ be the set of $k$ largest (in magnitude) coefficients of $u$, then

where we used Lemma 16 in the last line. It follows that

Consider any two vectors $u,v$, such that for $y=v-u$ we have $\|\bm{\Phi}y\|_{1}=\beta\geq 0$, and $\|v\|_{1}\leq\|u\|_{1}$. Let $S$ be the set of $k$ largest (in magnitude) coefficients of $u$. Then

Experimental Results

Our theoretical analysis shows that, up to constant factors, our scheme achieves the best known bounds for sparse approximate recovery. In order to determine the exact values of those constant factors, we show, in Figure 4, the empirical probability of correct recovery of a random $k$-sparse signal of dimension $n=200$ as a function of $k=\rho m$ and $m=\delta n$. As one can verify, the empirical $O(\cdot)$ constants involved are quite low. The thick curve shows the analytic computation of the phase transition between the survival of typical $l$-faces of the cross-polytope (left) and the polytope (right) under projection by $G$ a Gaussian random matrix. This line is equivalent to a phase transition in the ability of LP decoding to find the sparsest solution to $Gx_{*}=Gx$, and, in effect, is representative of the performance of Gaussian matrices in this framework (see [Don06b] and [DT06] for more details). Gaussian measurement matrices with $m=\delta n$ rows and $n$ columns can recover signals with sparsity $k=\rho m$ below the thick curve and cannot recover signals with sparsity $k$ above the curve. This figure thus shows that the empirical behavior of binary sparse matrices with LP decoding is consistent with the analytic performance of Gaussian random matrices. Furthermore, the empirical values of the asymptotic constants seem to agree. See [BI08] for further experimental data.

Sublinear-time decoding of sparse vectors

In this section we focus on the recovery problem for the case where the signal vector $x$ is $k$-sparse. The algorithm given here is essentially subsumed by the HHS algorithm provided in the Appendix (modulo the slightly higher reconstruction time and the number of measurements of the latter). The advantage of the algorithm from this section, however, is that it is very simple to present as well as analyze. This will enable us to illustrate the basic concepts without getting into technical details needed for the case of general signals $x$.

The algorithm we present here is essentially a simplification of the algorithm given in [Ind08]. The main difference is the assumption about the measurement matrix $\bm{\Phi}$. In [Ind08] the matrix $\bm{\Phi}$ was assumed to be an “augmented” version of the adjacency matrix of an extractor. Here, we assume that is an “augmented” version of the adjacency matrix of an unbalanced expander (or alternatively, by Theorem 2, any binary matrix that satisfies an RIP-1 property). The latter assumption turns out to be the “right” one, resulting in further simplification of the algorithm.

To define the matrices, we start with a straightforward definition. If $q$ and $r$ are 0–1 vectors, we can view them as masks that determine which entries of a signal appear and which are zeroed out. For example, the signal $q\circ x$ consists of the entries of $x$ restricted to the nonzero components of $q$. The notation $\circ$ indicates the componentwise or Hadamard product of two vectors. Given 0–1 matrices $\bm{Q}$ and $\bm{R}$, we can form a matrix that encodes sequential restrictions by all pairs of their rows. We express this matrix as the row tensor product $\bm{Q}\otimes_{\rm r}\bm{R}$.

The construction of the measurement matrix is as follows. Let $G=(A,B,E)$ be a $(k,\epsilon)$-unbalanced expander, $\epsilon=1/8$, with left degree $d$, $|A|=n$, $|B|=m$. Without loss of generality we can assume that $m$ is a power of $2$. Let $\bm{\Psi}$ be the adjacency matrix of $G$. The measurement matrix $\bm{\Phi}=\bm{\Psi}\otimes_{\rm r}\bm{B}$ has $m\log n$ rows. The matrix $\bm{B}$ is the bit-test matrix; its $t$th column contains the binary expansion of $t$. That is, For $t=0\ldots\log n-1$ and $j=0\ldots m-1$, let $l=t+j\log n$. The $(l+1)$-th row $\bm{\Phi}_{l+1}$ of $\bm{\Phi}$ is such that, for any $i=0\ldots m-1$, $(\bm{\Phi}_{l+1})_{i+1}=(\bm{\Psi})_{j+1}\cdot\mbox{bin}_{t}(i)$, where $\mbox{bin}_{t}(i)$ is the $t$-th least significant bit in the binary representation of $i$.

The purpose of the above augmentation is as follows. Consider a specific $k$-sparse vector $x$. Let $\operatorname{supp}(x)$ denote the the support of $x$, i.e., the set of indices of the non-zero coordinates of $x$. Assume that the $(j+1)$-th row $r$ of the matrix $\bm{\Psi}$ is such that $|\operatorname{supp}(x)\cap\operatorname{supp}(r)|=1$, and let $i+1\in\operatorname{supp}(x)\cap\operatorname{supp}(r)$. Then, from the values of $\bm{\Phi}_{j\log n+1}x,\ldots,\bm{\Phi}_{j\log n+\log n}x$, we can recover both $i$ and $x_{i+1}$. This is because for each of those rows, the value of $\bm{\Phi}_{j\log n+t+1}x$ is equal to $0$ if $\mbox{bin}_{t}(i)=0$, or $x_{i+1}$ otherwise; moreover, at least one of the values is non-zero. Of course, if $|\operatorname{supp}(x)\cap\operatorname{supp}(r)|\neq 1$, then the algorithm might “recover” an incorrect pair $(i,val)$. This can be detected if $val=0$; otherwise, the error might be undetected.

Let $\bm{\Phi}$ be a matrix described above. Then for any $k$-sparse $x$, given $\bm{\Phi}x$, we can recover $x$ in time $O(m\log^{2}n)$.

The main component of the recovery algorithm is the procedure ${\sc Reduce}(\bm{\Phi}x)$, which returns a vector $y$ such that $\|x-y\|_{0}\leq\|x\|_{0}/2$. The procedure maintains a vector $votes[\cdot]$ of multisets. Initially, all entries are set to $\emptyset$.

Reduce($\bm{\Phi}x$): for $j=1\ldots m$ compute $(l,val)$ such that if $\operatorname{supp}(\bm{\Phi})_{j}\cap\operatorname{supp}(x)=\{i\}$ then $x_{i}=val$ if $val\neq 0$ then $votes[l]=votes[l]\cup\{val\}$ $y=0$ for $i=1\ldots n$ such that $votes[i]\neq\emptyset$ if $votes[i]$ contains $\geq d/2$ copies of $val$ then $y_{i}=val$ return $y$

For any $k$-sparse vector $x$, the procedure Reduce returns a vector $y$ such that $x-y$ is $k/2$-sparse.

From the expanding properties of $G$, it follows that at most $\epsilon d\|x\|_{0}$ rows of $\bm{\Psi}$ return an incorrect pair $(i,val)$. Each set of $d/2$ incorrect pairs can change the outcome for at most $2$ positions of $y$. Thus, at most $4\epsilon\|x\|_{0}=\|x\|_{0}/2$ positions of $y$ can be incorrect. ∎

The final algorithm invokes Recover$(\bm{\Phi}x)$.

Recover($\bm{\Phi}x$) if $\bm{\Phi}x=0$ then return $0$ $y=$ Reduce($\bm{\Phi}x$) $\{$ We now need to recover $x-y\}$ $z=$Recover$(\bm{\Phi}(x-y))$ return $y+z$

Conclusion

We show in this paper that the geometric and the combinatorial approaches to sparse signal recovery are different manifestations of a common underyling phenomenon. Thus, we are able to show a unified perspective on both approaches—the key unifiying elements are the adjacency matrices of unbalanced expanders.

In most of the recent applications of compressed sensing, a physical device instantiates the measurement of $x$ and, as such, these applications need measurement matrices which are conducive to physical measurement processes. This paper shows that there is another, quite different, large, natural class of measurement matrices, combined with the same (or similar) recovery algorithms for sparse signal approximation. These measurement matrices may or may not be conducive to physical measurement processes but they are quite amenable to computational or digital signal measurement. Our work suggests a number of applications in digital or computational “sensing” such as efficient numerical linear algebra and network coding.

The preliminary experimental analysis exhibits interesting high-dimensional geometric phenomena as well. Our results suggest that the projection of polytopes under Gaussian random matrices is similar to that of projection by sparse random matrices, despite the fact that Gaussian random matrices are quite different from sparse ones.

References

Appendix A HHS(p)(p) decoding

In this section we focus on the general recovery problem for an arbitrary vector $x$. We show that, as in the previous algorithm in Section 5, we can use the adjacency matrices of unbalanced expanders, suitably augmented and concatenated, to obtain an explicit construction of matrices that are designed for the sub-linear time combinatorial recovery algorithm HHS in [GSTV07]. We call this modified algorithm the HHS$(p)$ algorithm. For the sake of exposition, we highlight the necessary changes in the construction of the matrices and the algorithm and summarize the remaining, unchanged portions of the algorithm. Similarly, for the analysis of HHS$(p)$, we present only those portions which are significantly different from the original analysis and summarize those which are not. We establish the following result.

Let $R$ denote the size of the measurements for either an explicit or random construction. Then, the HHS$(p)$ algorithm runs in time $\operatorname{poly}(R)$.

For explicit constructions $R=O(k2^{\log\log n^{O(1)}})$ and for random constructions $R=O(k\operatorname{polylog}n)$.

If we truncate the output $\widehat{x}$ of the algorithm to the $k$ largest terms, then the modified output $\widehat{x}_{k}$ satisfies

This subsection describes how to construct the measurement matrix $\bm{\Psi}$. We note that this construction has the same super-structure as the random construction in [GSTV07] and, as such, we highlight the necessary changes while summarizing the similar pieces. For ease of exposition, we set $\epsilon=1$. The matrix $\bm{\Psi}$ consists of two pieces: an identification matrix $\bm{\Omega}$ and an estimation matrix $\bm{\Phi}$. We use the first part of the sketch to identify indices of significant components of the signal quickly and then we use the second part to estimate the coefficients of those terms.

We use concatenated copies of the (normalized) adjacency matrices $\bm{\Gamma}$ of $(k^{\prime},\epsilon^{\prime})$-unbalanced expanders with left degree $d$. In what follows, we let the sparsity parameter $k^{\prime}$ vary but fix $p=1+1/\log(n)$. Normalize by the factor $d^{-1/p}$. By Theorem 1, such matrices satisfy the ${\rm RIP}_{{p},{k^{\prime}},{\delta}}$ property for $p=1+1/\log n$ and $\delta\leq C\epsilon^{\prime}$. Fix $\epsilon^{\prime}$ so that $4\epsilon^{\prime}+1/(2d^{1/p})<1/2$. Finally, we note that the number of rows in such a matrix is $k^{\prime}2^{(\log\log n)^{O(1)}}$ for known explicit constructions of such matrices and is $k^{\prime}\operatorname{polylog}(n)$ for random constructions. To simplify our accounting below, we denote the number of rows by $R$ for either the explicit or the random constructions.

The identification matrix $\bm{\Omega}$ is a 0–1 matrix with dimensions ${\rm O}(R\operatorname{polylog}(n))\times n$. It consists of a combination of a three structured matrices. Formally, $\bm{\Omega}$ is the row tensor product $\bm{\Omega}=\bm{B}\otimes_{\rm r}\bm{A}$. The bit-test matrix $\bm{B}$ has dimensions ${\rm O}(\log n)\times n$, and the isolation matrix $\bm{A}$ has dimensions ${\rm O}(R\operatorname{polylog}n)\times n$.

The isolation matrix $\bm{A}$ is a 0–1 matrix with a hierarchical structure. It consists of $\log(k)$ blocks $\bm{A}^{(j)}$ labeled by $j=1,2,4,8,\dots,J$, where $J=k$. Each block, in turn, has further substructure as a row tensor product of a collection of 0–1 matrices: $\bm{A}^{(j)}_{r,s}=\bm{R}^{(j)}_{r}\otimes_{\rm r}\bm{S}^{(j)}_{s}$ for $s=1,2,4,\ldots,k$ and $r=2s,4s,8s,\ldots,n$. The first matrix $\bm{S}^{(j)}_{s}$ is called the sifting matrix and the second matrix $\bm{R}^{(j)}_{r}$ is called the noise reduction matrix. Each sifting matrix $\bm{S}^{(j)}_{s}$ is the normalized adjacency matrix of an $(s,\epsilon^{\prime})$-unbalanced expander (or, equivalently, a normalized 0–1 ${\rm RIP}_{{p},{s},{\delta}}$ matrix).

The purpose of the noise reduction matrix $\bm{R}^{(j)}_{r}$ is to attenuate the noise in a signal that has a single large component. It is also a 0–1 valued matrix constructed as follows. For each pair of indices $(r,s)$, let $\beta\approx r/s$ be a prime power corresponding to a finite field of size $\beta$. Then we set each noise reduction matrix to be the Nisan-Wigderson generator over the field of size $\beta$. That is, we index the rows of the matrix by by pairs $(a,b)$ of field elements and index columns by polynomials $q$ of degree at most $\operatorname{polylog}(n)$. Position $((a,b),q)$ is 1, if $q(a)=b$, and 0, otherwise. Such a construction produces a matrix with $O(\beta^{2})$ rows. See [DeV07] or [CM06] for similar constructions of similar sizes.

We note that the dimensions of the product of matrices $\bm{A}^{(j)}_{r,s}=\bm{R}^{(j)}_{r}\otimes_{\rm r}\bm{S}^{(j)}_{s}$ are the critical ones (not the dimensions of the individual matrices). Each block is of dimension ${\rm O}(R\operatorname{polylog}(n))\times n$.

A.1.2. The estimation matrix

A.2. The HHS(p)(p) algorithm

A.3. Proof sketches

The goal of the algorithm is to identify a small set of signal components that carry most of the $p$-energy in the signal and to estimate the coefficients of those components well. We argue that, when our signal estimate is poor, the algorithm makes substantial progress toward this goal. We focus on the analysis of the algorithm in the case when our signal estimate is poor as this is the critical case. More precisely, assume that the current approximation ${a}$ satisfies

First, we obtain a fundamental relationship between the large and the small entries in a signal in Lemmas 21, 22.

Next, we show in Lemma 23 that the sifting matrix isolates significant coefficients in a moderate amount of noise.

Finally, we use a small RIPp,L,δ matrix to estimate the values of the coefficients in Lemma 28.

We note that for $p=1+\iota$, the constant $C_{p}=\Big(\frac{1}{p-1}\Big)^{1/p}$ is bounded below by $1/\iota$ and from above by $1/\iota^{2}$ so if $p=1+1/\log(n)$, then $\log(n)\leq C_{p}\leq\log^{2}(n)$.

If we combine this relation with our previous calculation, we find

To prove the second inequality, we use the Cauchy-Schwarz inequality to see that

Finally, we add this inequality to (A.4) to obtain

Now, we turn to the identification of significant signal entries. We pull these off in bands of decreasing magnitude. For exposition, we assume that $\|{g}\|_{1}=1$. We can think of the action of one submatrix $\bm{S}^{(j)}_{t}$ as

mapping each input signal to a collection of output signals.

To prove this result, we must show a sequence of shorter results. Our goal is to isolate significant spikes from one another and to reduce the contribution of the rest of the signal to these isolations.

An $(s,r)$ band is a band of spikes of magnitude between $1/r$ and $1/(2r)$. There is some $s$ (a power of 2) such that the number of spikes in the band is between $s$ and $2s$. If the $p$-contribution, $sr^{-p}$, of this band to the current residual error $\|{x}-{a}\|_{p}^{p}$ is greater than the average band’s contribution, $\frac{1}{\log n}\Bigg(k^{1/p-1}\|{x}-{a}\|_{p}\Bigg)^{p}$; i.e.,

The next lemma tells us how many spikes lie above a significant band.

If an $(s,r)$ band is significant, then there are at most $s\operatorname{polylog}(n)$ terms of magnitude greater than $1/r$ in the current residual.

If the $s$ spikes of magnitude $1/r$ have $p$-energy at least as big as $s^{\prime}$ spikes of magnitude $1/r^{\prime}>1/r$, then

and so the number of larger spikes must be less than $s$, $s^{\prime}\leq s$, as $p>1$. Because there are only $\log n$ possible $s^{\prime}$, if we take a union over all $s^{\prime}$, we achieve the desired result. ∎

We say that a spike with index $i$ is isolated from other spikes with indices in a set $I$ by a matrix $\bm{\Gamma}$ if there is a row in $\bm{\Gamma}$ that has a one in column $i$ and no other one in any column in $I$.

We are ready to begin the proof of Lemma 23. Our goal is to show that $\bm{\Gamma}=\bm{S}^{(j)}_{s}$ isolates from each other the vast majority of spikes of magnitude $\geq 1/r$, and, hence, isolates the majority of those with magnitude equal to $1/r$ from the larger spikes.

Let ${h}={g}_{i}{e}_{i}+\nu$ with $\|\nu\|_{1}\leq 1/s$. We apply the noise reduction operator $\bm{R}=\bm{R}^{(j)}_{r}$ to ${h}$. In at least half of the output signals, we have signals of the form

Consider the submatrix $\bm{R^{\prime}}$ of $\bm{R}$ obtained by extracting the $\beta$ rows of $\bm{R}$ that have a 1 in column $i$ and then removing column $i$ itself. By construction, $\bm{R}^{\prime}$ has at most $\operatorname{polylog}(n)$ ones in any column, so its 1-to-1 operator norm is at most $\operatorname{polylog}(n)$. Thus the average over rows of $\bm{R}^{\prime}$ of $\left\|{\bm{R}^{\prime}\nu}\right\|_{1}$ is at most $\frac{\operatorname{polylog}(n)}{\beta s}\approx 1/r\operatorname{polylog}(n)$, or at most $1/2r$ if we incorporate extra log factors into $\beta$.

In addition, the number of rows in the matrix $\bm{A}^{(j)}_{r,s}$ is $s\cdot\max(1,(r/s)^{2})$, where $s\leq k$ and $sr^{-p}\geq k^{p-1}$. If $s/r\geq 1$, the number of rows is $s\leq k$. If $s/r\leq 1$, the number of rows is $r^{2}/s$. Observe that $sr^{-p}\geq k^{1-p}$ and $r/s\geq 1$ imply $1/r\geq 1/k$. Since $sr^{-p}\geq k^{1-p}$, it follows that

The estimation portion of the algorithm is exactly the same as the HHS algorithm. Let $L$ be the set of identified candidates, $|L|\leq K$. By the RIP-$p$ property, the submatrix of $\Phi$ given by columns indexed by $L$ is invertible and the inverse $\Phi_{L}^{+}$ is bounded in the appropriate operator norm. we can apply $\Phi_{L}^{+}$ to $\Phi x$ by Jacobi iteration in time $O(K^{2}\operatorname{polylog}(n))$ to estimate the coefficients in $L$. In the exact case of $x_{\overline{L}}=0$, we recover $x$ exactly. Otherwise, the RIP-$p$ property assures that the $p$-norm of the vector of estimates is bounded in terms of the $p$-norm of $x-x_{K}$; i.e., this estimation process introduces errors that are small compared with the error associated with missing the $n-K$ other terms in $x$. Its proof is a small modification of the proof in [GSTV07], which itself was originally proven by Rudelson.

Then $\left\|{\widetilde{x}_{k}-x}\right\|_{p}\leq\left\|{x_{k}-x}\right\|_{p}+2\epsilon k^{1/p-1}\left\|{x_{k}-x}\right\|_{1}$ and $\left\|{\widetilde{x}_{k}-x}\right\|_{1}\leq(1+3\epsilon)\left\|{x_{k}-x}\right\|_{1}$.

Let $H$ denote $\operatorname{supp}(x_{k})\cup\operatorname{supp}(\widetilde{x}_{k})$, so that $|H|\leq 2k$. Let $H^{{\rm c}}$ denote the complement of $H$.