Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit

Introduction

Sparse recovery problems arise in many applications ranging from medical imaging to error correction. Suppose $v$ is an unknown $d$-dimensional signal with at most $n\ll d$ nonzero components:

As discussed in , exact recovery is possible with just $N=2n$. However, recovery using only this property is not numerically feasible; the sparse recovery problem in general is known to be NP-hard. Nevertheless, massive recent work in the emerging area of Compressed Sensing demonstrated that for several natural classes of measurement matrices $\Phi$, the signal $v$ can be exactly reconstructed from its measurements $\Phi v$ with

In other words, the number of measurements $N\ll d$ should be almost linear in the sparsity $n$. Survey contains some of these results; the Compressed Sensing webpage documents progress in this area.

The two major algorithmic approaches to sparse recovery are methods based on $L_{1}$-minimization and iterative methods (Matching Pursuits). We now briefly describe these methods. Then we propose a new iterative method that has advantages of both approaches.

This approach to sparse recovery has been advocated over decades by Donoho and his collaborators (see e.g. ). The sparse recovery problem can be stated as the problem of finding the sparsest signal $v$ with the given measurements $\Phi v$:

where $\|u\|_{0}:=|{\rm supp}(u)|$. Donoho and his associates advocated the principle that for some measurement matrices $\Phi$, the highly non-convex combinatorial optimization problem $(L_{0})$ should be equivalent to its convex relaxation

The recent progress in the emerging area of Compressed Sensing pushed forward this program (see survey ). A necessary and sufficient condition of exact sparse recovery is that the map $\Phi$ be one-to-one on the set of $n$-sparse vectors. Candès and Tao proved that a stronger quantitative version of this condition guarantees the equivalence of the problems $(L_{0})$ and $(L_{1})$.

A measurement matrix $\Phi$ satisfies the Restricted Isometry Condition (RIC) with parameters $(m,\varepsilon)$ for $\varepsilon\in(0,1)$ if we have

The Restricted Isometry Condition states that every set of $m$ columns of $\Phi$ forms approximately an orthonormal system. One can interpret the Restricted Isometry Condition as an abstract version of the Uniform Uncertainty Principle in harmonic analysis (, see also discussions in and ).

Assume that the measurement matrix $\Phi$ satisfies the Restricted Isometry Condition with parameters $(3n,0.2)$. Then every $n$-sparse vector $x$ can be exactly recovered from its measurements $\Phi x$ as a unique solution to the convex optimization problem $(L_{1})$.

In a lecture on Compressive Sampling, Candès sharpened this to work for the Restricted Isometry Condition with parameters $(2n,\sqrt{2}-1)$. Measurement matrices that satisfy the Restricted Isometry Condition with number of measurements as in (1.1) include random Gaussian, Bernoulli and partial Fourier matrices. Section 2 contains more detailed information.

2. Orthogonal Matching Pursuit (OMP)

An alternative approach to sparse recovery is via iterative algorithms, which find the support of the $n$-sparse signal $v$ progressively. Once $S={\rm supp}(v)$ is found correctly, it is easy to compute the signal $v$ from its measurements $x=\Phi v$ as $v=(\Phi_{S})^{-1}x$, where $\Phi_{S}$ denotes the measurement matrix $\Phi$ restricted to columns indexed by $S$.

A basic iterative algorithm is Orthogonal Matching Pursuit (OMP), popularized and analyzed by Gilbert and Tropp in , see for a more general setting. OMP recovers the support of $v$, one index at a time, in $n$ steps. Under a hypothetical assumption that $\Phi$ is an isometry, i.e. the columns of $\Phi$ are orthonormal, the signal $v$ can be exactly recovered from its measurements $x=\Phi v$ as $v=\Phi^{*}x$.

The problem is that the $N\times d$ matrix $\Phi$ is never an isometry in the interesting range where the number of measurements $N$ is smaller than the ambient dimension $d$. Even though the matrix is not an isometry, one can still use the notion of coherence in recovery of sparse signals. In that setting, greedy algorithms are used with incoherent dictionaries to recover such signals, see , , . In our setting, for random matrices one expects the columns to be approximately orthogonal, and the observation vector $u=\Phi^{*}x$ to be a good approximation to the original signal $v$.

The biggest coordinate of the observation vector $u$ in magnitude should thus be a nonzero coordinate of the signal $v$. We thus find one point of the support of $v$. Then OMP can be described as follows. First, we initialize the residual $r=x$. At each iteration, we compute the observation vector $u=\Phi^{*}r$. Denoting by $I$ the coordinates selected so far, we solve a least squares problem and update the residual

to remove any contribution of the coordinates in $I$. OMP then iterates this procedure $n$ times, and outputs a set $I$ of size $n$, which should equal the support of the signal $v$.

Tropp and Gilbert analyzed the performance of OMP for Gaussian measurement matrices $\Phi$; a similar result holds for general subgaussian matrices. They proved that, for every fixed $n$-sparse $d$-dimensional signal $v$, and an $N\times d$ random Gaussian measurement matrix $\Phi$, OMP recovers (the support of) $v$ from the measurements $x=\Phi v$ correctly with high probability, provided the number of measurements is $N\sim n\log d$.

3. Advantages and challenges of both approaches

The $L_{1}$-minimization method has strongest known guarantees of sparse recovery. Once the measurement matrix $\Phi$ satisfies the Restricted Isometry Condition, this method works correctly for all sparse signals $v$. No iterative methods have been known to feature such uniform guarantees, with the exception of Chaining Pursuit and the HHS Algorithm which however only work with specifically designed structured measurement matrices.

The Restricted Isometry Condition is a natural abstract deterministic property of a matrix. Although establishing this property is often nontrivial, this task is decoupled from the analysis of the recovery algorithm.

$L_{1}$-minimization is based on linear programming, which has its advantages and disadvantages. One thinks of linear programming as a black box, and any development of fast solvers will reduce the running time of the sparse recovery method. On the other hand, it is not very clear what this running time is, as there is no strongly polynomial time algorithm in linear programming yet. All known solvers take time polynomial not only in the dimension of the program $d$, but also on certain condition numbers of the program. While for some classes of random matrices the expected running time of linear programming solvers can be bounded (see the discussion in and subsequent work in ), estimating condition numbers is hard for specific matrices. For example, there is no result yet showing that the Restricted Isometry Condition implies that the condition numbers of the corresponding linear program is polynomial in $d$.

Orthogonal Matching Pursuit is quite fast, both theoretically and experimentally. It makes $n$ iterations, where each iteration amounts to a multiplication by a $d\times N$ matrix $\Phi^{*}$ (computing the observation vector $u$), and solving a least squares problem in dimensions at most $N\times n$ (with matrix $\Phi_{I}$). This yields strongly polynomial running time. In practice, OMP is observed to perform faster and is easier to implement than $L_{1}$-minimization . For more details, see .

Orthogonal Matching Pursuit is quite transparent: at each iteration, it selects a new coordinate from the support of the signal $v$ in a very specific and natural way. In contrast, the known $L_{1}$-minimization solvers, such as the simplex method and interior point methods, compute a path toward the solution. However, the geometry of $L_{1}$ is clear, whereas the analysis of greedy algorithms can be difficult simply because they are iterative.

On the other hand, Orthogonal Matching Pursuit has weaker guarantees of exact recovery. Unlike $L_{1}$-minimization, the guarantees of OMP are non-uniform: for each fixed sparse signal $v$ and not for all signals, the algorithm performs correctly with high probability. Rauhut has shown that uniform guarantees for OMP are impossible for natural random measurement matrices .

Moreover, OMP’s condition on measurement matrices given in is more restrictive than the Restricted Isometry Condition. In particular, it is not known whether OMP succeeds in the important class of partial Fourier measurement matrices.

These open problems about OMP, first stated in and often reverberated in the Compressed Sensing community, motivated the present paper. We essentially settle them in positive by the following modification of Orthogonal Matching Pursuit.

4. Regularized OMP

This new algorithm for sparse recovery will perform correctly for all measurement matrices $\Phi$ satisfying the Restricted Isometry Condition, and for all sparse signals.

When we are trying to recover the signal $v$ from its measurements $x=\Phi v$, we can use the observation vector $u=\Phi^{*}x$ as a good local approximation to the signal $v$. Namely, the observation vector $u$ encodes correlations of the measurement vector $x$ with the columns of $\Phi$. Note that $\Phi$ is a dictionary, and so since the signal $v$ is sparse, $x$ has a sparse representation with respect to the dictionary. By the Restricted Isometry Condition, every $n$ columns form approximately an orthonormal system. Therefore, every $n$ coordinates of the observation vector $u$ look like correlations of the measurement vector $x$ with the orthonormal basis and therefore are close in the Euclidean norm to the corresponding $n$ coefficients of $v$. This is documented in Proposition 3.2 below.

The local approximation property suggests to make use of the $n$ biggest coordinates of the observation vector $u$, rather than one biggest coordinate as OMP did. We thus force the selected coordinates to be more regular (ie. closer to uniform) by selecting only the coordinates with comparable sizes. To this end, a new regularization step will be needed to ensure that each of these coordinates gets an even share of information. This leads to the following algorithm for sparse recovery:

Regularized Orthogonal Matching Pursuit (ROMP)

The identification and regularization steps of ROMP can be performed efficiently. In particular, the regularization step does not imply combinatorial complexity, but actually can be done in linear time. The running time of ROMP is thus comparable to that of OMP in theory, and is often better than OMP in practice. We discuss the runtime in detail in Section 4.

The main theorem of this paper states that ROMP yields exact sparse recovery provided that the measurement matrix satisfies the Restricted Isometry Condition.

Remarks. 1. Theorem 1.3 guarantees exact sparse recovery. Indeed, it is easy to compute the signal $v$ from its measurements $x=\Phi v$ and the set $I$ given by ROMP as $v=(\Phi_{I})^{-1}x$, where $\Phi_{I}$ denotes the measurement matrix $\Phi$ restricted to columns indexed by $I$.

2. Theorem 1.3 gives uniform guarantees of sparse recovery. Indeed, once the measurement matrix satisfies a deterministic condition (RIC), then our algorithm ROMP correctly recovers every sparse vector from its measurements. Uniform guarantees have been shown to be impossible for OMP , and it has been an open problem to find a version of OMP with uniform guarantees (see ). Theorem 1.3 says that ROMP essentially settles this problem.

3. The logarithmic factor in $\varepsilon$ may be an artifact of the proof. At this moment, we do not know how to remove it.

4. Measurement matrices known to satisfy the Restricted Isometry Condition include random Gaussian, Bernoulli and partial Fourier matrices, with number of measurements $N$ almost linear in the sparsity $n$, i.e. as in (1.1). Section 2 contains detailed information. It has been unknown whether OMP gives sparse recovery for partial Fourier measurements (even with non-uniform guarantees). ROMP gives sparse recovery for these measurements, and even with uniform guarantees.

The rest of the paper is organized as follows. In Section 2 we describe known classes of measurement matrices satisfying the Restricted Isometry Condition. In Section 3 we give the proof of Theorem 1.3. In Section 4 we discuss implementation, running time, and empirical performance of ROMP.

Acknowledgment

We would like to thank the referees for a thorough reading of the manuscript and making useful suggestions which greatly improved the paper.

Measurement matrices satisfying the Restricted Isometry Condition

The only known measurement matrices known to satisfy the Restricted Isometry Condition with number of measurements as in (1.1) are certain classes of random matrices. The problem of deterministic constructions is still open. The known classes include: subgaussian random matrices (in particular, Gaussian and Bernoulli), and random partial bounded orthogonal matrices (in particular, partial Fourier matrices).

Throughout the paper, $C,c,C_{1},C_{2},c_{1},c_{2},\ldots$ denote positive absolute constants unless otherwise specified.

A partial bounded orthogonal matrix $\Phi$ is formed by $N$ randomly uniformly chosen rows of an orthogonal $d\times d$ matrix $\Psi$, whose entries are bounded by $C_{2}/\sqrt{d}$, for some constant $C_{2}$. An example of $\Psi$ is the discrete Fourier transform matrix. Taking measurements $\Phi v$ with a partial Fourier matrix thus amounts to observing $N$ random frequencies of the signal $v$.

The following theorem documents known results on the Restricted Isometry Condition for these classes of random matrices.

Consider an $N\times d$ measurement matrix $\Phi$, and let $n\geq 1$, $\varepsilon\in(0,1/2)$, and $\delta\in(0,1)$.

1. If $\Phi$ is a subgaussian matrix, then with probability $1-\delta$ the matrix $\frac{1}{\sqrt{N}}\Phi$ satisfies the Restricted Isometry Condition with parameters $(n,\varepsilon)$ provided that

2. If $\Phi$ is a partial bounded orthogonal matrix, then with probability $1-\delta$ the matrix $\sqrt{\frac{d}{N}}\,\Phi$ satisfies the Restricted Isometry Condition with parameters $(n,\varepsilon)$ provided that

In both cases, the constant $C$ depends only on the confidence level $\delta$ and the constants $C_{1},c_{1},C_{2}$ from the definition of the corresponding classes of matrices.

Remarks. 1. The first part of this theorem is proved in . The second part is from ; a similar estimate with somewhat worse exponents in the logarithms was proved in . See these results for the exact dependence of $C$ on the confidence level $\delta$ (although usually $\delta$ would be chosen to be some small constant itself.)

2. In Theorem 1.3, we needed to use RIC for $\varepsilon=c_{1}/\sqrt{\log n}$. An immediate consequence of Theorem 2.1 is that subgaussian matrices satisfy such RIC for the number of measurements

and partial bounded orthogonal matrices for

These numbers of measurements guarantee exact sparse recovery using ROMP.

Proof of Theorem 1.3

We shall prove a stronger version of Theorem 1.3, which states that at every iteration of ROMP, at least $50\%$ of the newly selected coordinates are from the support of the signal $v$.

Assume $\Phi$ satisfies the Restricted Isometry Condition with parameters $(2n,\varepsilon)$ for $\varepsilon=0.03/\sqrt{\log n}$. Let $v\neq 0$ be an $n$-sparse vector with measurements $x=\Phi v$. Then at any iteration of ROMP, after the regularization step, we have $J_{0}\neq\emptyset$, $J_{0}\cap I=\emptyset$ and

In other words, at least $50\%$ of the coordinates in the newly selected set $J_{0}$ belong to the support of $v$.

In particular, at every iteration ROMP finds at least one new coordinate in the support of the signal $v$. Coordinates outside the support can also be found, but (3.1) guarantees that the number of such “false” coordinates is always smaller than those in the support. This clearly implies Theorem 1.3.

Before proving Theorem 3.1 we explain how the Restricted Isometry Condition will be used in our argument. RIC is necessarily a local principle, which concerns not the measurement matrix $\Phi$ as a whole, but its submatrices of $n$ columns. All such submatrices $\Phi_{I}$, $I\subset\{1,\ldots,d\}$, $|I|\leq n$ are almost isometries. Therefore, for every $n$-sparse signal $v$, the observation vector $u=\Phi^{*}\Phi v$ approximates $v$ locally, when restricted to a set of cardinality $n$. The following proposition formalizes these local properties of $\Phi$ on which our argument is based.

Assume a measurement matrix $\Phi$ satisfies the Restricted Isometry Condition with parameters $(2n,\varepsilon)$. Then the following holds.

Since ${\rm supp}(v)\subset\Gamma$, we have

The conclusion of Part 1 follows since $I\subset\Gamma$.

Part 3. The desired inequality is equivalent to:

Let $K=I\cup J$ so that $|K|\leq 2n$. For any $x\in\range(\Phi_{I}),y\in\range(\Phi_{J})$, there are $a,b$ so that

By the proof of Part 2 above and since $\langle ab\rangle=0$, we have

The proof is by induction on the iteration of ROMP. The induction claim is that for all previous iterations, the set of newly chosen indices $J_{0}$ is nonempty, disjoint from the set of previously chosen indices $I$, and (3.1) holds.

Let $I$ be the set of previously chosen indices at the start of a given iteration. The induction claim easily implies that

Let $J_{0}$, $J$, be the sets found by ROMP in the current iteration. By the definition of the set $J_{0}$, it is nonempty.

Let $r\neq 0$ be the residual at the start of this iteration. We shall approximate $r$ by a vector in $\range(\Phi_{{\rm supp}(v)\setminus I})$. That is, we want to approximately realize the residual $r$ as measurements of some signal which lives on the still unfound coordinates of the the support of $v$. To that end, we consider the subspace

The Restricted Isometry Condition in the form of Part 3 of Proposition 3.2 ensures that $F$ and $E_{0}$ are almost orthogonal. Thus $E_{0}$ is close to the orthogonal complement of $F$ in $H$,

We will also consider the signal we seek to identify at the current iteration, its measurements, and its observation vector:

Lemma 3.5 will show that $\|(u-u_{0})|_{T}\|_{2}$ for any small enough subset $T$ is small, and Lemma 3.8 will show that $\|u|_{J_{0}}\|_{2}$ is not too small. First, we show that the residual $r$ has a simple description:

By definition of the residual in the algorithm, $r=P_{F^{\perp}}x$. Since $x\in H$, we conclude from the orthogonal decomposition $H=F+E$ that $x=P_{F}x+P_{E}x$. Thus $r=x-P_{F}x=P_{E}x$. ∎

To guarantee a correct identification of $v_{0}$, we first state two approximation lemmas that reflect in two different ways the fact that subspaces $E_{0}$ and $E$ are close to each other. This will allow us to carry over information from $E_{0}$ to $E$.

By definition of $F$, we have $x-x_{0}=\Phi(v-v_{0})\in F$. Therefore, by Lemma 3.3, $r=P_{E}x=P_{E}x_{0}$, and so

Now we use Part 3 of Proposition 3.2 for the sets $I$ and ${\rm supp}(v)\setminus I$ whose union has cardinality at most $2n$ by (3.2). It follows that $\|P_{F}P_{E_{0}}x_{0}\|_{2}\leq 2.2\varepsilon\|x_{0}\|_{2}$ as desired. ∎

Consider the observation vectors $u_{0}=\Phi^{*}x_{0}$ and $u=\Phi^{*}r$. Then for any set $T\subset\{1,\ldots,d\}$ with $|T|\leq 2n$, we have

Since $x_{0}=\Phi v_{0}$, we have by Lemma 3.4 and the Restricted Isometry Condition that

To complete the proof, it remains to apply Part 2 of Proposition 3.2, which yields $\|(u_{0}-u)|_{T}\|_{2}\leq(1+\varepsilon)\|x_{0}-r\|_{2}$. ∎

We next show that the energy (norm) of $u$ when restricted to $J$, and furthermore to $J_{0}$, is not too small. By the approximation lemmas, this will yield that ROMP selects at least a fixed percentage of energy of the still unidentified part of the signal. By the regularization step of ROMP, since all selected coefficients have comparable magnitudes, we will conclude that not only a portion of energy but also of the support is selected correctly. This will be the desired conclusion.

We have $\|u|_{J}\|_{2}\geq 0.8\|v_{0}\|_{2}$.

Let $S$ = ${\rm supp}(v)\setminus I$. Since $|S|\leq n$, the maximality property of $J$ in the algorithm implies that

Furthermore, since $v_{0}|_{S}=v_{0}$, by Part 1 of Proposition 3.2 we have

Putting these two inequalities together and using Lemma 3.5, we conclude that

We next bound the norm of $u$ restricted to the smaller set $J_{0}$. We do this by first noticing a general property of regularization:

We will construct at most $O(\log m)$ subsets $A_{k}$ with comparable coordinates as in (3.4), and such that at least one of these sets will have large energy as in (3.5).

Let $y=(y_{1},\ldots,y_{m})$, and consider a partition of $\{1,\ldots,m\}$ using sets with comparable coordinates:

Let $k_{0}=\left\lceil\log m\right\rceil+1$, so that $|y_{i}|\leq\frac{1}{m}\|y\|_{2}$ for all $i\in A_{k}$, $k>k_{0}$. Then the set $U=\bigcup_{k\leq k_{0}}A_{k}$ contains most of the energy of $y$:

Therefore there exists $k\leq k_{0}$ such that

In our context, Lemma 3.7 applied to the vector $u|_{J}$ along with Lemma 3.6 directly implies:

To show the first claim, that $J_{0}$ is nonempty, we note that $v_{0}\neq 0$. Indeed, otherwise by (3.3) we have $I\subset{\rm supp}(v)$, so by the definition of the residual in the algorithm, we would have $r=0$ at the start of the current iteration, which is a contradiction. Then $J_{0}\neq\emptyset$ by Lemma 3.8.

The second claim, that $J_{0}\cap I=\emptyset$, is also simple. Indeed, recall that by the definition of the algorithm, $r=P_{F^{\perp}}\in F^{\perp}=(\range(\Phi_{I}))^{\perp}$. It follows that the observation vector $u=\Phi^{*}r$ satisfies $u|_{I}=0$. Since by its definition the set $J$ contains only nonzero coordinates of $u$ we have $J\cap I=\emptyset$. Since $J_{0}\subset J$, the second claim $J_{0}\cap I=\emptyset$ follows.

The nontrivial part of the theorem is its last claim, inequality (3.1). Suppose it fails. Namely, suppose that $|J_{0}\cap{\rm supp}(v)|<\frac{1}{2}|J_{0}|$, and thus

Set $\Lambda=J_{0}\backslash{\rm supp}(v)$. By the comparability property of the coordinates in $J_{0}$ and since $|\Lambda|>\frac{1}{2}|J_{0}|$, there is a fraction of energy in $\Lambda$:

where the last inequality holds by Lemma 3.8.

On the other hand, we can approximate $u$ by $u_{0}$ as

Since $\Lambda\subset J$ and using Lemma 3.5, we have

Furthermore, by definition (3.3) of $v_{0}$, we have $v_{0}|_{\Lambda}=0$. So, by Part 1 of Proposition 3.2,

Using the last two inequalities in (3.7), we conclude that

This is a contradiction to (3.6) so long as $\varepsilon\leq 0.03/\sqrt{\log n}$. This proves Theorem 3.1. ∎

Implementation and empirical performance of ROMP

The Identification step of ROMP, i.e. selection of the subset $J$, can be done by sorting the coordinates of $u$ in the nonincreasing order and selecting $n$ biggest. Many sorting algorithms such as Mergesort or Heapsort provide running times of $O(d\log d)$.

The Regularization step of ROMP, i.e. selecting $J_{0}\subset J$, can be done fast by observing that $J_{0}$ is an interval in the decreasing rearrangement of coefficients. Moreover, the analysis of the algorithm shows that instead of searching over all intervals $J_{0}$, it suffices to look for $J_{0}$ among $O(\log n)$ consecutive intervals with endpoints where the magnitude of coefficients decreases by a factor of $2$. (these are the sets $A_{k}$ in the proof of Lemma 3.7). Therefore, the Regularization step can be done in time $O(n)$.

In addition to these costs, the $k$-th iteration step of ROMP involves multiplication of the $d\times N$ matrix $\Phi^{*}$ by a vector, and solving the least squares problem with the $N\times|I|$ matrix $\Phi_{I}$, where $|I|\leq 2k\leq 2n$. For unstructured matrices, these tasks can be done in time $dN$ and $O(n^{2}N)$ respectively. Since the submatrix of $\Phi$ when restricted to the index set $I$ is near an isometry, using an iterative method such as the Conjugate Gradient Method allows us to solve the least squares method in a constant number of iterations (up to a specific accuracy.) Using such a method then reduces the time of solving the least squares problem to just $O(nN)$. Thus in the cases where ROMP terminates after a fixed number of iterations, the total time to solve all required least squares problems would be just $O(nN)$. For structured matrices, such as partial Fourier, these times can be improved even more using fast multiply techniques.

In other cases, however, ROMP may need more than a constant number of iterations before terminating, say the full $O(n)$ iterations. In this case, it may be more efficient to maintain the QR factorization of $\Phi_{I}$ and use the Modified Gram-Schmidt algorithm. With this method, solving all the least squares problems takes total time just $O(n^{2}N).$ However, storing the QR factorization is quite costly, so in situations where storage is limited it may be best to use the iterative methods mentioned above.

ROMP terminates in at most $2n$ iterations. Therefore, for unstructured matrices using the methods mentioned above and in the interesting regime $N\geq\log d$, the total running time of ROMP is O(dNn). This is the same bound as for OMP .

2. Non-sparse signals

In many applications, one needs to recover a signal $v$ which is not sparse but close to being sparse in some way. Such are, for example, compressible signals, whose coefficients decay at a certain rate (see , ). To make ROMP work for such signals, one can replace the stopping criterion of exact recovery $r=0$ by “repeat $n$ times or until $r=0$, whichever occurs first”. Note that we could amend the algorithm for sparse signals in this way as well, allowing for a specific level of accuracy to be attained before terminating.

We recently proved that ROMP is stable and guarantees approximate recovery of non-sparse signals with noisy measurements; this will be discussed in a forthcoming paper.

3. Experiments

First we describe the setup of our experiments. For many values of the ambient dimension $d$, the number of measurements $N$, and the sparsity $n$, we reconstruct random signals using ROMP. For each set of values, we generate an $N\times d$ Gaussian measurement matrix $\Phi$ and then perform $500$ independent trials. The results we obtained using Bernoulli measurement matrices were very similar. In a given trial, we generate an $n$-sparse signal $v$ in one of two ways. In either case, we first select the support of the signal by choosing $n$ components uniformly at random (independent from the measurement matrix $\Phi$). In the cases where we wish to generate flat signals, we then set these components to one. Our work as well as the analysis of Gilbert and Tropp show that this is a challenging case for ROMP (and OMP). In the cases where we wish to generate sparse compressible signals, we set the $i^{th}$ component of the support to plus or minus $i^{-1/p}$ for a specified value of $0<p<1$. We then execute ROMP with the measurement vector $x=\Phi v$.

Figure 1 depicts the percentage (from the $500$ trials) of sparse flat signals that were reconstructed exactly. This plot was generated with $d=256$ for various levels of sparsity $n$. The horizontal axis represents the number of measurements $N$, and the vertical axis represents the exact recovery percentage. We also performed this same test for sparse compressible signals and found the results very similar to those in Figure 1. Our results show that performance of ROMP is very similar to that of OMP which can be found in .

Figure 2 depicts a plot of the values for $N$ and $n$ at which $99\%$ of sparse flat signals are recovered exactly. This plot was generated with $d=256$. The horizontal axis represents the number of measurements $N$, and the vertical axis the sparsity level $n$.

Theorem 1.3 guarantees that ROMP runs with at most $O(n)$ iterations. Figure 3 depicts the number of iterations executed by ROMP for $d=10,000$ and $N=200$. ROMP was executed under the same setting as described above for sparse flat signals as well as sparse compressible signals for various values of $p$, and the number of iterations in each scenario was averaged over the $500$ trials. These averages were plotted against the sparsity of the signal. As the plot illustrates, only $2$ iterations were needed for flat signals even for sparsity $n$ as high as $40$. The plot also demonstrates that the number of iterations needed for sparse compressible is higher than the number needed for sparse flat signals, as one would expect. The plot suggests that for smaller values of $p$ (meaning signals that decay more rapidly) ROMP needs more iterations. However it shows that even in the case of $p=0.5$, only $6$ iterations are needed even for sparsity $n$ as high as $20$.

References