Model-Based Compressive Sensing

I Introduction

We are in the midst of a digital revolution that is enabling the development and deployment of new sensors and sensing systems with ever increasing fidelity and resolution. The theoretical foundation is the Shannon/Nyquist sampling theorem, which states that a signal’s information is preserved if it is uniformly sampled at a rate at least two times faster than its Fourier bandwidth. Unfortunately, in many important and emerging applications, the resulting Nyquist rate can be so high that we end up with too many samples and must compress in order to store or transmit them. In other applications the cost of signal acquisition is prohibitive, either because of a high cost per sample, or because state-of-the-art samplers cannot achieve the high sampling rates required by Shannon/Nyquist. Examples include radar imaging and exotic imaging modalities outside visible wavelengths.

Transform compression systems reduce the effective dimensionality of an $N$-dimensional signal $x$ by re-representing it in terms of a sparse or compressible set of coefficients $\alpha$ in a basis expansion $x=\Psi\alpha$, with $\Psi$ an $N\times N$ basis matrix. By sparse we mean that only $K\ll N$ of the coefficients $\alpha$ are nonzero and need to be stored or transmitted. By compressible we mean that the coefficients $\alpha$, when sorted, decay rapidly enough to zero that $\alpha$ can be well-approximated as $K$-sparse. The sparsity and compressibility properties are pervasive in many signal classes of interest. For example, smooth signals and images are compressible in the Fourier basis, while piecewise smooth signals and images are compressible in a wavelet basis ; the JPEG and JPEG2000 standards are examples of practical transform compression systems based on these bases.

Compressive sensing (CS) provides an alternative to Shannon/Nyquist sampling when the signal under acquisition is known to be sparse or compressible . In CS, we measure not periodic signal samples but rather inner products with $M\ll N$ measurement vectors. In matrix notation, the measurements $y=\Phi x=\Phi\Psi\alpha$, where the rows of the $M\times N$ matrix $\Phi$ contain the measurement vectors. While the matrix $\Phi\Psi$ is rank deficient, and hence loses information in general, it can be shown to preserve the information in sparse and compressible signals if it satisfies the so-called restricted isometry property (RIP) . Intriguingly, a large class of random matrices have the RIP with high probability. To recover the signal from the compressive measurements $y$, we search for the sparsest coefficient vector $\alpha$ that agrees with the measurements. To date, research in CS has focused primarily on reducing both the number of measurements $M$ (as a function of $N$ and $K$) and on increasing the robustness and reducing the computational complexity of the recovery algorithm. Today’s state-of-the-art CS systems can robustly recover $K$-sparse and compressible signals from just $M=\mathcal{O}\left(K\log(N/K)\right)$ noisy measurements using polynomial-time optimization solvers or greedy algorithms.

While this represents significant progress from Nyquist-rate sampling, our contention in this paper is that it is possible to do even better by more fully leveraging concepts from state-of-the-art signal compression and processing algorithms. In many such algorithms, the key ingredient is a more realistic structured sparsity model that goes beyond simple sparsity by codifying the inter-dependency structure among the signal coefficients $\alpha$. Obviously, sparsity and compressibility correspond to simple signal models where each coefficient is treated independently; for example in a sparse model, the fact that the coefficient $\alpha_{i}$ is large has no bearing on the size of any $\alpha_{j}$, $j\neq i$. We will reserve the use of the term “model” for situations where we are enforcing structured dependencies between the values and the locations of the coefficients $\alpha_{i}$. For instance, modern wavelet image coders exploit not only the fact that most of the wavelet coefficients of a natural image are small but also the fact that the values and locations of the large coefficients have a particular structure. Coding the coefficients according to a structured sparsity model enables these algorithms to compress images close to the maximum amount possible – significantly better than a naïve coder that just processes each large coefficient independently. We have previously developed a new CS recovery algorithm that promotes structure in the sparse representation by tailoring the recovered signal according to a sparsity-promoting probabilistic model, such as an Ising graphical model . Such probabilistic models favor certain configurations for the magnitudes and indices of the significant coefficients of the signal.

In this paper, we expand on this concept by introducing a model-based CS theory that parallels the conventional theory and provides concrete guidelines on how to create structured signal recovery algorithms with provable performance guarantees. By reducing the number of degrees of freedom of a sparse/compressible signal by permitting only certain configurations of the large and zero/small coefficients, structured sparsity models provide two immediate benefits to CS. First, they enable us to reduce, in some cases significantly, the number of measurements $M$ required to stably recover a signal. Second, during signal recovery, they enable us to better differentiate true signal information from recovery artifacts, which leads to a more robust recovery.

To precisely quantify the benefits of model-based CS, we introduce and study several new theoretical concepts that could be of more general interest. We begin with structured sparsity models for $K$-sparse signals and make precise how the structure reduces the number of potential sparse signal supports in $\alpha$. Then using the model-based restricted isometry property from , we prove that such structured sparse signals can be robustly recovered from noisy compressive measurements. Moreover, we quantify the required number of measurements $M$ and show that for some structured sparsity models $M$ is independent of $N$. These results unify and generalize the limited related work to date on structured sparsity models for strictly sparse signals . We then introduce the notion of a structured compressible signal, whose coefficients $\alpha$ are no longer strictly sparse but have a structured power-law decay. To establish that structured compressible signals can be robustly recovered from compressive measurements, we generalize the standard RIP to a new restricted amplification property (RAmP). Using the RAmP, we show that the required number of measurements $M$ for recovery of structured compressible signals is independent of $N$.

To take practical advantage of this new theory, we demonstrate how to integrate structured sparsity models into two state-of-the-art CS recovery algorithms, CoSaMP and iterative hard thresholding (IHT) . The key modification is surprisingly simple: we merely replace the nonlinear sparse approximation step in these greedy algorithms with a structured sparse approximation. Thanks to our new theory, both new model-based recovery algorithms have provable robustness guarantees for both structured sparse and structured compressible signals.

To validate our theory and algorithms and demonstrate their general applicability and utility, we present two specific instances of model-based CS and conduct a range of simulation experiments. The first structured sparsity model accounts for the fact that the large wavelet coefficients of piecewise smooth signals and images tend to live on a rooted, connected tree structure . Using the fact that the number of such trees is much smaller than $N\choose K$, the number of $K$-sparse signal supports in $N$ dimensions, we prove that a tree-based CoSaMP algorithm needs only $M=\mathcal{O}\left(K\right)$ measurements to robustly recover tree-sparse and tree-compressible signals. This provides a significant reduction against the standard CS requirement $M=\mathcal{O}\left(K\log(N/K)\right)$ as the signal length $N$ increases. Figure 1 indicates the potential performance gains on a tree-compressible, piecewise smooth signal.

The second structured sparsity model accounts for the fact that the large coefficients of many sparse signals cluster together . Such a so-called block sparse model is equivalent to a joint sparsity model for an ensemble of $J$, length-$N$ signals , where the supports of the signals’ large coefficients are shared across the ensemble. Using the fact that the number of clustered supports is much smaller than $JN\choose JK$, we prove that a block-based CoSaMP algorithm needs only $M=\mathcal{O}\left(JK+K\log(\frac{N}{K})\right)$ measurements to robustly recover block-sparse and block-compressible signals. In contrast, standard CS requires $M=\mathcal{O}\left(JK\log(N/K)\right)$; block sparsity reduces the dependence of $M$ on the signal length $N$, particularly for large block sizes $J$.

Our new theory and methods relate to a small body of previous work aimed at integrating structured sparsity into CS. Several groups have developed structured sparse signal recovery algorithms ; however, their approaches have either been ad hoc or focused on a single structured sparsity model. Most previous work on unions of subspaces has focused exclusively on strictly sparse signals and has considered neither compressibility nor feasible recovery algorithms. A related CS modeling framework for structured sparse and compressible signals collects the $N$ samples of a signal into $D$ groups, $D\leq N$, and allows signals where $K$ out of $D$ groups have nonzero coefficients. This framework is immediately applicable to block-sparse signals and signal ensembles with common sparse supports. While provides recovery algorithms, measurement bounds, and recovery guarantees similar to those provided in Section VI, our proposed framework has the ability to focus on arbitrary subsets of the $D\choose K$ groups that yield more elaborate structures, such as connected subtrees for wavelet coefficients. To the best of our knowledge, our general algorithmic framework for model-based recovery, the concept of a model-compressible signal, and the associated RAmP are new to the literature.

This paper is organized as follows. A review of the CS theory in Section II lays out the foundational concepts that we extend to the model-based case in subsequent sections. Section III develops the concept of structured sparse signals and introduces the concept of structured compressible signals. We also quantify how structured sparsity models improve the measurement and recovery process by exploiting the model-based RIP for structured sparse signals and by introducing the RAmP for structured compressible signals. Section IV indicates how to tune CoSaMP to incorporate structured sparsity models and establishes its robustness properties for structured sparse and structured compressible signals; the modifications to the IHT algorithm are very similar, so we defer them to an appendix to reduce redundancy. Sections V and VI then specialize our theory to the special cases of wavelet tree and block sparse signal models, respectively, and report on a series of numerical experiments that validate our theoretical claims. We conclude with a discussion in Section VII. To make the paper more readable, all proofs are relegated to a series of appendices.

II Background on Compressive Sensing

Many natural and manmade signals are not strictly sparse, but can be approximated as such; we call such signals compressible. Consider a signal $x$ whose coefficients, when sorted in order of decreasing magnitude, decay according to the power law

The approximation of compressible signals by sparse signals is the basis of transform coding as is used in algorithms like JPEG and JPEG2000 . In this framework, we acquire the full $N$-sample signal $x$; compute the complete set of transform coefficients $\alpha$ via $\alpha=\Psi^{-1}x$; locate the $K$ largest coefficients and discard the $(N-K)$ smallest coefficients; and encode the $K$ values and locations of the largest coefficients. While a widely accepted standard, this sample-then-compress framework suffers from three inherent inefficiencies. First, we must start with a potentially large number of samples $N$ even if the ultimate desired $K$ is small. Second, the encoder must compute all of the $N$ transform coefficients $\alpha$, even though it will discard all but $K$ of them. Third, the encoder faces the overhead of encoding the locations of the large coefficients.

II-B Compressive measurements and the restricted isometry property

Compressive sensing (CS) integrates the signal acquisition and compression steps into a single process . In CS we do not acquire $x$ directly but rather acquire $M<N$ linear measurements $y=\Phi x$ using an $M\times N$ measurement matrix $\Phi$. We then recover $x$ by exploiting its sparsity or compressibility. Our goal is to push $M$ as close as possible to $K$ in order to perform as much signal “compression” during acquisition as possible.

In order to recover a good estimate of $x$ (the $K$ largest $x_{i}$’s, for example) from the $M$ compressive measurements, the measurement matrix $\Phi$ should satisfy the restricted isometry property (RIP) .

An $M\times N$ matrix $\Phi$ has the $K$-restricted isometry property ($K$-RIP) with constant $\delta_{K}$ if, for all $x\in\Sigma_{K}$,

In words, the $K$-RIP ensures that all submatrices of $\Phi$ of size $M\times K$ are close to an isometry, and therefore distance (and information) preserving. Practical recovery algorithms typically require that $\Phi$ have a slightly stronger $2K$-RIP, $3K$-RIP, or higher-order RIP in order to preserve distances between $K$-sparse vectors (which are $2K$-sparse in general), three-way sums of $K$-sparse vectors (which are $3K$-sparse in general), and other higher-order structures.

II-C Recovery algorithms

Since there are infinitely many signal coefficient vectors $x^{\prime}$ that produce the same set of compressive measurements $y=\Phi x$, to recover the “right” signal we exploit our a priori knowledge of its sparsity or compressibility. For example, we could seek the sparsest $x$ that agrees with the measurements $y$:

This corresponds to a linear program that can be solved in polynomial time . Adaptations to deal with additive noise in $y$ or $x$ include basis pursuit with denoising (BPDN) , complexity-based regularization , and the Dantzig Selector .

The second approach finds the sparsest $x$ agreeing with the measurements $y$ through an iterative, greedy search. Algorithms such as matching pursuit, orthogonal matching pursuit , StOMP , iterative hard thresholding (IHT) , CoSaMP , and Subspace Pursuit (SP) all revolve around a best $L$-term approximation for the estimated signal, with $L$ varying for each algorithm; typically $L$ is $\mathcal{O}\left(K\right)$.

II-D Performance bounds on signal recovery

For a noise-free, $K$-sparse signal, these algorithms offer perfect recovery, meaning that the signal $\widehat{x}$ recovered from the compressive measurements $y=\Phi x$ is exactly $\widehat{x}=x$.

For a $K$-sparse signal $x$ whose measurements are corrupted by noise $n$ of bounded norm (that is, we measure $y=\Phi x+n$) the mean-squared error of the signal $\widehat{x}$ is

For an $s$-compressible signal $x$ whose measurements are corrupted by noise $n$ of bounded norm, the mean-squared error of the recovered signal $\widehat{x}$ is

Using (3) we can simplify this expression to

III Structured Sparsity and Compressibility

While many natural and manmade signals and images can be described to first-order as sparse or compressible, the support of their large coefficients often has an underlying inter-dependency structure. This phenomenon has received only limited attention by the CS community to date . In this section, we introduce a model-based theory of CS that captures such structure. A model reduces the degrees of freedom of a sparse/compressible signal by permitting only certain configurations of supports for the large coefficient. As we will show, this allows us to reduce, in some cases significantly, the number of compressive measurements $M$ required to stably recover a signal.

To state a formal definition of a structured sparsity model, let $x|_{\Omega}$ represent the entries of $x$ corresponding to the set of indices $\Omega\subseteq\{1,\ldots,N\}$, and let $\Omega^{C}$ denote the complement of the set $\Omega$.

A structured sparsity model $\mathcal{M}_{K}$ is defined as the union of $m_{K}$ canonical $K$-dimensional subspaces

Signals from $\mathcal{M}_{K}$ are called $K$-structured sparse. Clearly, $\mathcal{M}_{K}\subseteq\Sigma_{K}$ and contains $m_{K}\leq{N\choose K}$ subspaces.

In Sections V and VI below we consider two concrete structured sparsity models. The first model accounts for the fact that the large wavelet coefficients of piecewise smooth signals and images tend to live on a rooted, connected tree structure . The second model accounts for the fact that the large coefficients of sparse signals often cluster together into blocks .

III-B Model-based RIP

If we know that the signal $x$ being acquired is $K$-structured sparse, then we can relax the RIP constraint on the CS measurement matrix $\Phi$ and still achieve stable recovery from the compressive measurements $y=\Phi x$ .

An $M\times N$ matrix $\Phi$ has the $\mathcal{M}_{K}$-restricted isometry property ($\mathcal{M}_{K}$-RIP) with constant $\delta_{\mathcal{M}_{K}}$ if, for all $x\in\mathcal{M}_{K}$, we have

Blumensath and Davies have quantified the number of measurements $M$ necessary for a random CS matrix to have the $\mathcal{M}_{K}$-RIP with a given probability.

where $c$ is a positive constant, an $M\times N$ i.i.d. subgaussian random matrix has the $\mathcal{M}_{K}$-RIP with constant $\delta_{\mathcal{M}_{K}}$ with probability at least $1-e^{-t}$.

This bound can be used to recover the conventional CS result by substituting $m_{K}={N\choose K}\approx(Ne/K)^{K}$. Similarly, as the number of subspaces $m_{K}$ that arise from the structure imposed can be significantly smaller than the standard ${N\choose K}$, the number of rows needed for a random matrix to have the $\mathcal{M}_{K}$-RIP can be significantly lower than the number of rows needed for the standard RIP. The $\mathcal{M}_{K}$-RIP property is sufficient for robust recovery of structured sparse signals, as we show below in Section IV-B.

III-C Structured compressible signals

The set of $s$-structured compressible signals is defined as

Define $|x|_{\mathfrak{M}_{s}}$ as the smallest value of $G$ for which this condition holds for $x$ and $s$.

We say that $x\in\mathfrak{M}_{s}$ is an $s$-structured compressible signal under the structured sparsity model $\mathcal{M}_{K}$. These approximation classes have been characterized for certain structured sparsity models; see Section V for an example. We will select the value of $s$ for which the distance between the approximation errors $\sigma_{\mathcal{M}_{K}}(x)$ and the corresponding bounds $GK^{-1/s}$ is minimal.

III-D Nested model approximations and residual subspaces

In conventional CS, the same requirement (RIP) is a sufficient condition for the stable recovery of both sparse and compressible signals. In model-based recovery, however, the class of structured compressible signals is much larger than that of structured sparse signals, since the union of subspaces defined by structured sparse signals does not contain all canonical $K$-dimensional subspaces.

To address this difference, we introduce some additional tools to develop a sufficient condition for the stable recovery of structured compressible signals. We will pay particular attention to structured sparsity models $\mathcal{M}_{K}$ that generate nested approximations, since they are more amenable to analysis and computation.

In words, a structured sparsity model generates nested approximations if the support of the best $K^{\prime}$-term structured sparse approximation contains the support of the best $K$-term structured sparse approximation for all $K<K^{\prime}$. An important example of a NAP-generating structured sparse model is the standard compressible signal model of (3).

When a structured sparsity model obeys the NAP, the support of the difference between the best $jK$-term structured sparse approximation and the best $(j+1)K$-term structured sparse approximation of a signal can be shown to lie in a small union of subspaces, thanks to the structure enforced by the model. This structure is captured by the set of subspaces that are included in each subsequent approximation, as defined below.

Under the NAP, each structured compressible signal $x$ can be partitioned into its best $K$-term structured sparse approximation $x_{T_{1}}$, the additional components present in the best $2K$-term structured sparse approximation $x_{T_{2}}$, and so on, with $x=\sum_{j=1}^{\lceil N/K\rceil}x_{T_{j}}$ and $x_{T_{j}}\in\mathcal{R}_{j,K}(\mathcal{M})$ for each $j$. Each signal partition $x_{T_{j}}$ is a $K$-sparse signal, and thus $\mathcal{R}_{j,K}(\mathcal{M})$ is a union of subspaces of dimension $K$. We will denote by $R_{j}$ the number of subspaces that compose $\mathcal{R}_{j,K}(\mathcal{M})$ and omit the dependence on $\mathcal{M}$ in the sequel for brevity.

Intuitively, the norms of the partitions $\|x_{T_{j}}\|_{2}$ decay as $j$ increases for signals that are structured compressible. As the next subsection shows, this observation is instrumental in relaxing the isometry restrictions on the measurement matrix $\Phi$ and bounding the recovery error for $s$-structured compressible signals when the model obeys the NAP.

III-E The restricted amplification property (RAmP)

For exactly $K$-structured sparse signals, we discussed in Section III-B that the number of compressive measurements $M$ required for a random matrix to have the $\mathcal{M}_{K}$-RIP is determined by the number of canonical subspaces $m_{K}$ via (1). Unfortunately, such structured sparse concepts and results do not immediately extend to structured compressible signals. Thus, we develop a generalization of the $\mathcal{M}_{K}$-RIP that we will use to quantify the stability of recovery for structured compressible signals.

One way to analyze the robustness of compressible signal recovery in conventional CS is to consider the tail of the signal outside its $K$-term approximation as contributing additional “noise” to the measurements of size $\|\Phi(x-x_{K})\|_{2}$ . Consequently, the conventional $K$-sparse recovery performance result can be applied with the augmented noise $n+\Phi(x-x_{K})$.

This technique can also be used to quantify the robustness of structured compressible signal recovery. The key quantity we must control is the amplification of the structured sparse approximation residual through $\Phi$. The following property is a new generalization of the RIP and model-based RIP.

A matrix $\Phi$ has the $(\epsilon_{K},r)$-restricted amplification property (RAmP) for the residual subspaces $\mathcal{R}_{j,K}$ of model $\mathcal{M}$ if

for any $u\in\mathcal{R}_{j,K}$ for each $1\leq j\leq\lceil N/K\rceil$.

The regularity parameter $r>0$ caps the growth rate of the amplification of $u\in\mathcal{R}_{j,K}$ as a function of $j$. Its value can be chosen so that the growth in amplification with $j$ balances the decay of the norm in each residual subspace $\mathcal{R}_{j,K}$ with $j$.

We can quantify the number of compressive measurements $M$ required for a random measurement matrix $\Phi$ to have the RAmP with high probability; we prove the following in Appendix A.

Let $\Phi$ be an $M\times N$ matrix with i.i.d. subgaussian entries and let the set of residual subspaces $\mathcal{R}_{j,K}$ of the structured sparsity model $\mathcal{M}$ contain $R_{j}$ subspaces of dimension $K$ for each $1\leq j\leq\lceil N/K\rceil$. If

then the matrix $\Phi$ has the $(\epsilon_{K},r)$-RAmP with probability $1-e^{-t}$.

The order of the bound of Theorem 2 is lower than $\mathcal{O}\left(K\log(N/K)\right)$ as long as the number of subspaces $R_{j}$ grows slower than $N^{K}$.

Armed with the RaMP, we can state the following result, which will provide robustness for the recovery of structured compressible signals; see Appendix B for the proof.

Let $x\in\mathfrak{M}_{s}$ be an $s$-structured compressible signal under a structured sparsity model $\mathcal{M}$ that obeys the NAP. If $\Phi$ has the $(\epsilon_{K},r)$-RAmP and $r=s-1$, then we have

where $C_{s}$ is a constant that depends only on $s$.

IV Model-Based Signal Recovery Algorithms

To take practical advantage of our new theory for model-based CS, we demonstrate how to integrate structured sparsity models into two state-of-the-art CS recovery algorithms, CoSaMP (in this section) and iterative hard thresholding (IHT) (in Appendix C to avoid repetition). The key modification is simple: we merely replace the best $K$-term sparse approximation step in these greedy algorithms with a best $K$-term structured sparse approximation. Since at each iteration we need only search over the $m_{K}$ subspaces of $\mathcal{M}_{K}$ rather than ${N\choose K}$ subspaces of $\Sigma_{K}$, fewer measurements will be required for the same degree of robust signal recovery. Or, alternatively, using the same number of measurements, more accurate recovery can be achieved.

After presenting the modified CoSaMP algorithm, we prove robustness guarantees for both structured sparse and structured compressible signals. To this end, we must define an enlarged union of subspaces that includes sums of elements in the structured sparsity model.

The $B$-order sum for the set $\mathcal{M}_{K}$, with $B>1$ an integer, is defined as

IV-B Performance of structured sparse signal recovery

We now study the performance of model-based CoSaMP signal recovery on structured sparse and structured compressible signals. A robustness guarantee for noisy measurements of structured sparse signals can be obtained using the model-based RIP (10). Our performance guarantee for structured sparse signal recovery will require that the measurement matrix $\Phi$ be a near-isometry for all subspaces in $\mathcal{M}^{B}_{K}$ for some $B>1$. This requirement is a direct generalization of the $2K$-RIP, $3K$-RIP, and higher-order RIPs from the conventional CS theory. The following theorem is proven in Appendix D.

Let $x\in\mathcal{M}_{K}$ and let $y=\Phi x+n$ be a set of noisy CS measurements. If $\Phi$ has an $\mathcal{M}^{4}_{K}$-RIP constant of $\delta_{\mathcal{M}^{4}_{K}}\leq 0.1$, then the signal estimate $\widehat{x}_{i}$ obtained from iteration $i$ of the model-based CoSaMP algorithm satisfies

This guarantee matches that of the CoSaMP algorithm [11, Theorem 4.1]; however, our guarantee is only for structured sparse signals rather than for all sparse signals.

IV-C Performance of structured compressible signal recovery

Using the new tools introduced in Section III, we can provide a robustness guarantee for noisy measurements of structured compressible signals, using the RAmP as a condition on the measurement matrix $\Phi$.

Let $x\in\mathfrak{M}_{s}$ be an $s$-structured compressible signal from a structured sparsity model $\mathcal{M}$ that obeys the NAP, and let $y=\Phi x+n$ be a set of noisy CS measurements. If $\Phi$ has the $\mathcal{M}^{4}_{K}$-RIP with $\delta_{\mathcal{M}^{4}_{K}}\leq 0.1$ and the $(\epsilon_{K},r)$-RAmP with $\epsilon_{K}\leq 0.1$ and $r=s-1$, then the signal estimate $\widehat{x}_{i}$ obtained from iteration $i$ of the model-based CoSaMP algorithm satisfies

Proof sketch. To prove the theorem, we first bound the optimal structured sparse recovery error for an $s$-structured compressible signal $x\in\mathfrak{M}_{s}$ when the matrix $\Phi$ has the $(\epsilon_{K},r)$-RAmP with $r\leq s-1$ (see Theorem 3). Then, using Theorem 13, we can easily prove the result by following the analogous proof in . ∎

The standard CoSaMP algorithm also features a similar guarantee for structured compressible signals, with the constant changing from 35 to 20.

IV-D Robustness to model mismatch

We now analyze the robustness of model-based CS recovery to model mismatch, which occurs when the signal being recovered from compressive measurements does not conform exactly to the structured sparsity model used in the recovery algorithm.

We begin with optimistic results for signals that are “close” to matching the recovery structured sparsity model. First consider a signal $x$ that is not $K$-structured sparse as the recovery algorithm assumes but rather $(K+\kappa)$-structured sparse for some small integer $\kappa$. This signal can be decomposed into $x_{K}$, the signal’s $K$-term structured sparse approximation, and $x-x_{K}$, the error of this approximation. For $\kappa\leq K$, we have that $x-x_{K}\in\mathcal{R}_{2,K}$. If the matrix $\Phi$ has the $(\epsilon_{K},r)$-RAmP, then it follows than

Using equations (13) and (15), we obtain the following guarantee for the $i^{th}$ iteration of model-based CoSaMP:

By noting that $\|x-x_{K}\|_{2}$ is small, we obtain a guarantee that is close to (13).

Second, consider a signal $x$ that is not $s$-structured compressible as the recovery algorithm assumes but rather $(s-\epsilon)$-structured compressible. The following bound can be obtained under the conditions of Theorem 5 by modifying the argument in Appendix B:

As $\epsilon$ becomes smaller, the factor $\frac{\lceil N/K\rceil^{\epsilon}-1}{\epsilon}$ approaches $\log\lceil N/K\rceil$, matching (14). In summary, as long as the deviations from the structured sparse and structured compressible classes are small, our model-based recovery guarantees still apply within a small bounded constant factor.

IV-E Computational complexity of model-based recovery

The computational complexity of a structured signal recovery algorithm differs from that of a standard algorithm by two factors. The first factor is the reduction in the number of measurements $M$ necessary for recovery: since most current recovery algorithms have a computational complexity that is linear in the number of measurements, any reduction in $M$ reduces the total complexity. The second factor is the cost of the structured sparse approximation. The $K$-term approximation used in most current recovery algorithms can be implemented with a simple sorting operation ($\mathcal{O}\left(N\log N\right)$ complexity, in general). Ideally, the structured sparsity model should support a similarly efficient approximation algorithm.

To validate our theory and algorithms and demonstrate their general applicability and utility, we now present two specific instances of model-based CS and conduct a range of simulation experiments.

V Example: Wavelet Tree Model

Wavelet decompositions have found wide application in the analysis, processing, and compression of smooth and piecewise smooth signals because these signals are $K$-sparse and compressible, respectively . Moreover, the wavelet coefficients can be naturally organized into a tree structure, and for many kinds of natural and manmade signals the largest coefficients cluster along the branches of this tree. This motivates a connected tree model for the wavelet coefficients .

While CS recovery for wavelet-sparse signals has been considered previously , the resulting algorithms integrated the tree constraint in an ad-hoc fashion. Furthermore, the algorithms provide no recovery guarantees or bounds on the necessary number of compressive measurements.

We first describe tree sparsity in the context of sparse wavelet decompositions. We focus on one-dimensional signals and binary wavelet trees, but all of our results extend directly to $d$-dimensional signals and $2^{d}$-ary wavelet trees.

Consider a signal $x$ of length $N=2^{I}$, for an integer value of $I$. The wavelet representation of $x$ is given by

where $\nu$ is the scaling function and $\psi_{i,j}$ is the wavelet function at scale $i$ and offset $j$. The wavelet transform consists of the scaling coefficient $v_{0}$ and wavelet coefficients $w_{i,j}$ at scale $i$, $0\leq i\leq I-1$, and position $j$, $0\leq j\leq 2^{i}-1$. In terms of our earlier matrix notation, $x$ has the representation $x=\Psi\alpha$, where $\Psi$ is a matrix containing the scaling and wavelet functions as columns, and $\alpha=[v_{0}~w_{0,0}~w_{1,0}~w_{1,1}~w_{2,0}\ldots]^{T}$ is the vector of scaling and wavelet coefficients. We are, of course, interested in sparse and compressible $\alpha$.

The nested supports of the wavelets at different scales create a parent/child relationship between wavelet coefficients at different scales. We say that $w_{i-1,\lfloor j/2\rfloor}$ is the parent of $w_{i,j}$ and that $w_{i+1,2j}$ and $w_{i+1,2j+1}$ are the children of $w_{i,j}$. These relationships can be expressed graphically by the wavelet coefficient tree in Figure 2.

Wavelet functions act as local discontinuity detectors, and using the nested support property of wavelets at different scales, it is straightforward to see that a signal discontinuity will give rise to a chain of large wavelet coefficients along a branch of the wavelet tree from a leaf to the root. Moreover, smooth signal regions will give rise to regions of small wavelet coefficients. This “connected tree” property has been well-exploited in a number of wavelet-based processing and compression algorithms. In this section, we will specialize the theory developed in Sections III and IV to a connected tree model $\mathcal{T}$.

A set of wavelet coefficients $\Omega$ forms a connected subtree if, whenever a coefficient $w_{i,j}\in\Omega$, then its parent $w_{i-1,\lceil j/2\rceil}\in\Omega$ as well. Each such set $\Omega$ defines a subspace of signals whose support is contained in $\Omega$; that is, all wavelet coefficients outside $\Omega$ are zero. In this way, we define the structured sparsity model $\mathcal{T}_{K}$ as the union of all $K$-dimensional subspaces corresponding to supports $\Omega$ that form connected subtrees.

Define the set of $K$-tree sparse signals as

To quantify the number of subspaces in $\mathcal{T}_{K}$, it suffices to count the number of distinct connected subtrees of size $K$ in a binary tree of size $N$. We prove the following result in Appendix E.

The number of subspaces in $\mathcal{T}_{K}$ obeys $T_{K}\leq\frac{4^{K+4}}{Ke^{2}}$ for $K\geq\log_{2}N$ and $T_{K}\leq\frac{(2e)^{K}}{K+1}$ for $K<\log_{2}N$.

To simplify the presentation in the sequel, we will simply use the weaker bound $T_{K}\leq\frac{(2e)^{K}}{K+1}$ for all values of $K$ and $N$.

V-B Tree-based approximation

Fortuitously, an efficient solver exists, called the condensing sort and select algorithm (CSSA) . Recall that subtree approximation coincides with standard $K$-term approximation (and hence can be solved by simply sorting the wavelet coefficients) when the wavelet coefficients are monotonically nonincreasing along the tree branches out from the root. The CSSA solves (16) in the case of general wavelet coefficient values by condensing the nonmonotonic segments of the tree branches using an iterative sort-and-average routine during a greedy search through the nodes. For each node in the tree, the algorithm calculates the average wavelet coefficient magnitude for each subtree rooted at that node, and records the largest average among all the subtrees as the energy for that node. The CSSA then searches for the unselected node with the largest energy and adds the subtree corresponding to the node’s energy to the estimated support as a supernode: a single node that provides a condensed representation of the corresponding subtree . Condensing a large coefficient far down the tree accounts for the potentially large cost (in terms of the total budget of tree nodes $K$) of growing the tree to that point.

Since the first step of the CSSA involves sorting all of the wavelet coefficients, overall it requires $\mathcal{O}\left(N\log N\right)$ computations. However, once the CSSA grows the optimal tree of size $K$, it is trivial to determine the optimal trees of size $<K$ and computationally efficient to grow the optimal trees of size $>K$ .

The constrained optimization (16) can be rewritten as an unconstrained problem by introducing the Lagrange multiplier $\lambda$ :

where $\overline{\mathcal{T}}=\cup_{n=1}^{N}\mathcal{T}_{n}$ and $\bar{\alpha}$ are the wavelet coefficients of $\bar{x}$. Except for the inconsequential $\lambda K$ term, this optimization coincides with Donoho’s complexity penalized sum of squares , which can be solved in only $\mathcal{O}\left(N\right)$ computations using coarse-to-fine dynamic programming on the tree. Its primary shortcoming is the nonobvious relationship between the tuning parameter $\lambda$ and and the resulting size $K$ of the optimal connected subtree.

V-C Tree-compressible signals

Specializing Definition 2 from Section III-C to $\mathcal{T}$, we make the following definition.

Define the set of $s$-tree compressible signals as

Furthermore, define $|x|_{\mathfrak{T}_{s}}$ as the smallest value of $G$ for which this condition holds for $x$ and $s$.

Tree approximation classes contain signals whose wavelet coefficients have a loose (and possibly interrupted) decay from coarse to fine scales. These classes have been well-characterized for wavelet-sparse signals and are intrinsically linked with the Besov spaces $B^{s}_{q}(L_{p}())$. Besov spaces contain functions of one or more continuous variables that have (roughly speaking) $s$ derivatives in $L_{p}()$; the parameter $q$ provides finer distinctions of smoothness. When a Besov space signal $x_{a}\in B^{s}_{p}(L_{p}())$ with $s>1/p-1/2$ is sampled uniformly and converted to a length-$N$ vector $x$, its wavelet coefficients belong to the tree approximation space $\mathfrak{T}_{s}$, with

where “$\asymp$” denotes an equivalent norm. The same result holds if $s=1/p-1/2$ and $q\leq p$.

V-D Stable tree-based recovery from compressive measurements

For tree-sparse signals, by applying Theorem 1 and Proposition 1, we find that a subgaussian random matrix has the $\mathcal{T}_{K}$-RIP property with constant $\delta_{\mathcal{T}_{K}}$ and probability $1-e^{-t}$ if the number of measurements obeys

Thus, the number of measurements necessary for stable recovery of tree-sparse signals is linear in $K$, without the dependence on $N$ present in conventional non-model-based CS recovery.

For tree-compressible signals, we must quantify the number of subspaces $R_{j}$ in each residual set $\mathcal{R}_{j,K}$ for the approximation class. We can then apply the theory of Section IV-C with Proposition 1 to calculate the smallest allowable $M$ via Theorem 5.

The number of $K$-dimensional subspaces that comprise $\mathcal{R}_{j,K}$ obeys

Using Proposition 17 and Theorem 5, we obtain the following condition for the matrix $\Phi$ to have the RAmP, which is proved in Appendix F.

Let $\Phi$ be an $M\times N$ matrix with i.i.d. subgaussian entries. If

then the matrix $\Phi$ has the $(\epsilon_{K},s)$-RAmP for the structured sparsity model $\mathcal{T}$ and all $s>0.5$ with probability $1-e^{-t}$.

Both cases give a simplified bound on the number of measurements required as $M=\mathcal{O}\left(K\right)$, which is a substantial improvement over the $M=\mathcal{O}\left(K\log(N/K)\right)$ required by conventional CS recovery methods. Thus, when $\Phi$ satisfies Proposition 3, we have the guarantee (14) for sampled Besov space signals from $B^{s}_{q}(L_{p}())$.

V-E Experiments

We now present the results of a number of numerical experiments that illustrate the effectiveness of a tree-based recovery algorithm. Our consistent observation is that explicit incorporation of the structured sparsity model in the recovery process significantly improves the quality of recovery for a given number of measurements. In addition, model-based recovery remains stable when the inputs are no longer tree-sparse, but rather are tree-compressible and/or corrupted with differing levels of noise. We employ the Daubechies-6 wavelet basis for sparsity, and recover the signal using model-based CoSaMP (Algorithm 1) with a CSSA-based structured sparse approximation step in all experiments.

We first study one-dimensional signals that match the connected wavelet-tree model described above. Among such signals is the class of piecewise smooth functions, which are commonly encountered in analysis and practice.

Figure 3(a) illustrates the results of a Monte Carlo simulation study on the impact of the number of measurements $M$ on the performance of model-based and conventional recovery for a class of tree-sparse piecewise polynomial signals. Each data point was obtained by measuring the normalized recovery error of 500 sample trials. Each sample trial was conducted by generating a new piecewise polynomial signal of length $N=1024$ with five polynomial pieces of cubic degree and randomly placed discontinuities, computing its best $K$-term tree-approximation using the CSSA, and then measuring the resulting signal using a matrix with i.i.d. Gaussian entries. Model-based recovery attains near-perfect recovery at $M=3K$ measurements, while CoSaMP only matches this performance at $M=5K$.

For the same class of signals, we empirically compared the recovery times of our proposed algorithm with those of the standard approach (CoSaMP). Experiments were conducted on a Sun workstation with a 1.8GHz AMD Opteron dual-core processor and 2GB memory running UNIX, using non-optimized Matlab code and a function-handle based implementation of the random projection operator $\Phi$. As is evident from Figure 3(b), wavelet tree-based recovery is in general slower than CoSaMP. This is due to the fact that the CSSA step in the iterative procedure is more computationally demanding than simple $K-$term approximation. Nevertheless, the highest benefits of model-based CS recovery are obtained around $M=3K$; in this regime, the runtimes of the two approaches are comparable, with tree-based recovery requiring fewer iterations and yielding much smaller recovery error than standard recovery.

Figure 4 shows the growth of the overmeasuring factor $M/K$ with the signal length $N$ for conventional CS and model-based recovery. We generated 50 sample piecewise cubic signals and numerically computed the minimum number of measurements $M$ required for the recovery error $\|x-\widehat{x}\|_{2}\leq 2.5\sigma_{\mathcal{T}_{K}}(x)$, the best tree-approximation error, for every sample signal. The figure shows that while doubling the signal length increases the number of measurements required by standard recovery by $K$, the number of measurements required by model-based recovery is constant for all $N$. These experimental results verify the theoretical performance described in Proposition 3.

Finally, we turn to two-dimensional images and a wavelet quadtree model. The connected wavelet-tree model has proven useful for compressing natural images ; thus, our algorithm provides a simple and provably efficient method for recovering a wide variety of natural images from compressive measurements. An example of recovery performance is given in Figure 6. The test image (Peppers) is of size $N=128\times 128=16384$ pixels, and we computed $M=5000$ random Gaussian measurements. Model-based recovery again offers higher performance than standard signal recovery algorithms like CoSaMP, both in terms of recovery mean-squared error and visual quality.

VI Example: Block-Sparse Signals and Signal Ensembles

In a block-sparse signal, the locations of the significant coefficients cluster in blocks under a specific sorting order. Block-sparse signals have been previously studied in CS applications, including DNA microarrays and magnetoencephalography . An equivalent problem arises in CS for signal ensembles, such as sensor networks and MIMO communication . In this case, several signals share a common coefficient support set. For example, when a frequency-sparse acoustic signal is recorded by an array of microphones, then all of the recorded signals contain the same Fourier frequencies but with different amplitudes and delays. Such a signal ensemble can be re-shaped as a single vector by concatenation, and then the coefficients can be rearranged so that the concatenated vector exhibits block sparsity.

It has been shown that the block-sparse structure enables signal recovery from a reduced number of CS measurements, both for the single signal case and the signal ensemble case , through the use of specially tailored recovery algorithms. However, the robustness guarantees for the algorithms either are restricted to exactly sparse signals and noiseless measurements, do not have explicit bounds on the number of necessary measurements, or are asymptotic in nature. An optimization-based algorithm introduced in provides similar recovery guarantees to those obtained by the algorithm we present in this chapter; thus, our method can be interpreted as a greedy-based counterpart to that provided in .

In this section, we formulate the block sparsity model as a union of subspaces and pose an approximation algorithm on this union of subspaces. The approximation algorithm is used to implement block-based signal recovery. We also define the corresponding class of block-compressible signals and quantify the number of measurements necessary for robust recovery.

Define the set of $K$-block sparse signals as

VI-B Block-based approximation

To pose the block-based approximation algorithm, we need to define the mixed norm of a matrix.

The $(p,q)$ mixed norm of the matrix $X=[x_{1}~x_{2}~\ldots~x_{N}]$ is defined as

When $q=0$, $\|X\|_{(p,0)}$ simply counts the number of nonzero columns in $X$.

VI-C Block-compressible signals

We define the set of $s$-block compressible signals as

where $\mathcal{I}$ indexes the sorted column norms.

We say that $X$ is an $s$-block compressible signal if $X\in\mathfrak{S}_{s}$. For such signals, we have $\|X-X_{K}\|_{(2,2)}=\sigma_{\mathcal{S}_{K}}(x)\leq G_{1}K^{-s}$, and $\|X-X_{K}\|_{(2,1)}\leq G_{2}K^{1/2-s}$. Note that the block-compressible model does not impart a structure to the decay of the signal coefficients, so that the sets $\mathcal{R}_{j,K}$ are equal for all values of $j$; due to this property, the $(\delta_{\mathcal{S}_{K}},s)$-RAmP is implied by the $\mathcal{S}_{K}$-RIP. Taking this into account, we can derive the following result from , which is proven similarly to Theorem 13.

Let $x$ be a signal from the structured sparsity model $\mathcal{S}$, and let $y=\Phi x+n$ be a set of noisy CS measurements. If $\Phi$ has the $\mathcal{S}_{K}^{4}$-RIP with $\delta_{\mathcal{S}_{K}^{4}}\leq 0.1$, then the estimate obtained from iteration $i$ of block-based CoSaMP, using the approximation algorithm (18), satisfies

Thus, the algorithm provides a recovered signal of similar quality to approximations of $X$ by a small number of nonzero columns. When the signal $x$ is $K$-block sparse, we have that $||X-X_{K}^{\mathcal{S}}\|_{(2,2)}=||X-X_{K}^{\mathcal{S}}\|_{(2,1)}=0$, obtaining the same result as Theorem 13, save for a constant factor.

VI-D Stable block-based recovery from compressive measurements

Since Theorem 6 poses the same requirement on the measurement matrix $\Phi$ for sparse and compressible signals, the same number of measurements $M$ is required to provide performance guarantees for block-sparse and block-compressible signals. The class $\mathcal{S}_{K}$ contains $S={N\choose K}$ subspaces of dimension $JK$. Thus, a subgaussian random matrix has the $\mathcal{S}_{K}$-RIP property with constant $\delta_{\mathcal{S}_{K}}$ and probability $1-e^{-t}$ if the number of measurements obeys

To compare with the standard CS measurement bound, the number of measurements required for robust recovery scales as $M=\mathcal{O}\left(JK+K\log(N/K)\right)$, which is a substantial improvement over the $M=\mathcal{O}\left(JK\log(N/K)\right)$ that would be required by conventional CS recovery methods. When the size of the block $J$ is larger than $\log(N/K)$, then this term becomes $\mathcal{O}\left(KJ\right)$; that is, it is linear on the total sparsity of the block-sparse signal.

We note in passing that the bound on the number of measurements (19) assumes a dense subgaussian measurement matrix, while the measurement matrices used in have a block-diagonal structure. To obtain measurements from an $M\times JN$ dense matrix in a distributed setting, it suffices to partition the matrix into $J$ pieces of size $M\times N$ and calculate the CS measurements at each sensor with the corresponding matrix; these individual measurements are then summed to obtain the complete measurement vector. For large $J$, (19) implies that the total number of measurements required for recovery of the signal ensemble is lower than the bound for the case where each signal recovery is performed independently for each signal ($M=\mathcal{O}\left(JK\log(N/K)\right)$).

VI-E Experiments

Figure 7 illustrates an $N=4096$ signal that exhibits block sparsity, and its recovered version from $M=960$ measurements using CoSaMP and model-based recovery. The block size $J=64$ and there were $K=6$ active blocks in the signal. We observe the clear advantage of using the block-sparsity model in signal recovery.

Figure 9(a) indicates the decay in recovery error as a function of the numbers of measurements for CoSaMP and model-based recovery. We generated sample block-sparse signals as follows: we randomly selected a set of $K$ blocks, each of size $J$, and endow them with coefficients that follow an i.i.d. Gaussian distribution. Each sample point in the curves is generated by performing 200 trials of the corresponding algorithm. As in the connected wavelet-tree case, we observe clear gains using model-based recovery, particularly for low-measurement regimes; CoSaMP matches model-based recovery only for $M\geq 5K.$

Figure 9(b) compares the recovery times of the two approaches. For this particular model, we observe that our proposed approach is in general much faster than CoSaMP. This is because of two reasons: a) the block-based approximation step involves sorting fewer coefficients, and thus is faster than $K-$term approximation; b) block-based recovery requires fewer iterations to converge to the true solution.

VII Conclusions

In this paper, we have aimed to demonstrate that there are significant performance gains to be made by exploiting more realistic and richer signal models beyond the simplistic sparse and compressible models that dominate the CS literature. Building on the unions of subspaces results of and the proof machinery of , we have taken some first steps towards what promises to be a general theory for model-based CS by introducing the notion of a structured compressible signal and the associated restricted amplification property (RAmP) condition it imposes on the measurement matrix $\Phi$. Our analysis poses the nested approximation property (NAP) as a sufficient condition that is satisfied by many structured sparsity models.

For the volumes of natural and manmade signals and images that are wavelet-sparse or compressible, our tree-based CoSaMP and IHT algorithms offer performance that significantly exceeds today’s state-of-the-art while requiring only $M=\mathcal{O}\left(K\right)$ rather than $M=\mathcal{O}\left(K\log(N/K)\right)$ random measurements. For block-sparse signals and signal ensembles with common sparse support, our block-based CoSaMP and IHT algorithms offer not only excellent performance but also require just $M=\mathcal{O}\left(JK\right)$ measurements, where $JK$ is the signal sparsity. Furthermore, block-based recovery can recover signal ensembles using fewer measurements than the number required when each signal is recovered independently; we have shown such advantages using real-world data from environmental sensor networks . Additional structured sparsity models have been developed using our general framework in and ; we have also released a Matlab toolbox containing the corresponding model-based CS recovery algorithms, available at http://dsp.rice.edu/software.

Appendix A Proof of Theorem 2

To prove this theorem, we will study the distribution of the maximum singular value of a submatrix $\Phi_{T}$ of a matrix with i.i.d. Gaussian entries $\Phi$ corresponding to the columns indexed by $T$. From this we obtain the probability that RAmP does not hold for a fixed support $T$. We will then evaluate the same probability for all supports $T$ of elements of $\mathcal{R}_{j,K}$, where the desired bound on the amplification is dependent on the value of $j$. This gives us the probability that the RAmP does not hold for a given residual subspace set $\mathcal{R}_{j,K}$. We fix the probability of failure on each of these sets; we then obtain probability that the matrix $\Phi$ does not have the RAmP using a union bound. We end by obtaining conditions on the number of rows $M$ of $\Phi$ to obtain a desired probability of failure.

We begin from the following concentration of measure for the largest singular value of a $M\times K$ submatrix $\Phi_{T}$, $|T|=K$, of an $M\times N$ matrix $\Phi$ with i.i.d. subgaussian entries that are properly normalized :

For large enough $M$, $\beta\ll 1$; thus we ignore this small constant in the sequel. By letting $\tau=j^{r}\sqrt{1+\epsilon_{K}}-1-\sqrt{\frac{K}{M}}$ (with the appropriate value of $j$ for $T$), we obtain

We use a union bound over all possible $R_{j}$ supports for $u\in\mathcal{R}_{j,K}$ to obtain the probability that $\Phi$ amplifies the norm of some $u$ by more than $j^{r}\sqrt{1+\epsilon_{K}}$:

Bound the right hand side by a constant $\mu$; this requires

for each $j$. We use another union bound among the residual subspaces $\mathcal{R}_{j.K}$ to measure the probability that the RAmP does not hold:

To bound this probability by $e^{-t}$, we need $\mu=\frac{K}{N}e^{-t}$; plugging this into (20), we obtain

for each $j$. Simplifying, we obtain that for $\Phi$ to posess the RAmP with probability $1-e^{-t}$, the following must hold for all $j$:

Since $(\sqrt{a}+\sqrt{b})^{2}\leq 2a+2b$ for $a,b>0$, then the hypothesis (12) implies (21), proving the theorem. ∎

Appendix B Proof of Theorem 3

is the difference between the best $jK$ structured sparse approximation and the best $(j-1)K$ structured sparse approximation. Additionally, each piece $x_{T_{j}}\in\mathcal{R}_{j,K}$. Therefore, since $\Phi$ satisfies the $(\epsilon_{K},s-1)$-RAmP, we obtain

Since $x\in\mathfrak{M}_{s}$, the norm of each piece can be bounded as

It is easy to show, using Euler-Maclaurin summations, that $\sum_{j=2}^{\lceil N/K\rceil}j^{-1}\leq\ln\lceil N/K\rceil$; we then obtain

Appendix C Model-based Iterative Hard Thresholding

Our proposed model-based iterative hard thresholding (IHT) is given in Algorithm 2. For this algorithm, Theorems 13, 5, and 6 can be proven with only a few modifications: $\Phi$ must have the $\mathcal{M}_{K}^{3}$-RIP with $\delta_{\mathcal{M}_{K}^{3}}\leq 0.1$, and the constant factor in the bound changes from 15 to 4 in Theorem 13, from 35 to 10 in Theorem 5, and from 20 to 5 in Theorem 6.

To illustrate the performance of the algorithm, we repeat the HeaviSine experiment from Figure 1. Recall that $N=1024$, and $M=80$ for this example. The advantages of using our tree-structured sparse approximation step (instead of mere hard thresholding) are evident from Figure 10. In practice, we have observed that our model-based algorithm converges in fewer steps than IHT and yields much more accurate results in terms of recovery error.

Appendix D Proof of Theorem 13

The proof of this theorem is identical to that of the CoSaMP algorithm in [11, Section 4.6], and requires a set of six lemmas. The sequence of Lemmas 1–6 below are modifications of the lemmas in that are restricted to the structured sparsity model. Lemma 4 does not need any changes from , so we state it without proof. The proof of Lemmas 3–6 use the properties in Lemmas 1 and 2, which are simple to prove.

Suppose $\Phi$ has $\mathcal{M}$-RIP with constant $\delta_{\mathcal{M}}$. Let $\Omega$ be a support corresponding to a subspace in $\mathcal{M}$. Then we have the following handy bounds.

Suppose $\Phi$ has $\mathcal{M}^{2}_{K}$-RIP with constant $\delta_{\mathcal{M}^{2}_{K}}$. Let $\Omega$ be a support corresponding to a subspace in $\mathcal{M}_{K}$, and let $x\in\mathcal{M}_{K}$. Then $\|\Phi_{\Omega}^{T}\Phi x|_{\Omega^{C}}\|_{2}\leq\delta_{\mathcal{M}^{2}_{K}}\|x|_{\Omega^{C}}\|_{2}.$

We begin the proof of Theorem 13 by fixing an iteration $i\geq 1$ of model-based CoSaMP. We write $\widehat{x}=\widehat{x}_{i-1}$ for the signal estimate at the beginning of the $i^{th}$ iteration. Define the signal residual $s=x-\widehat{x}$, which implies that $s\in\mathcal{M}_{K}^{2}$. We note that we can write $r=y-\Phi\widehat{x}=\Phi(x-\widehat{x})+n=\Phi s+n$.

where $\Omega\setminus\Pi$ denotes the set difference of $\Omega$ and $\Pi$. These signals are in $\mathcal{M}_{K}^{4}$ (since they arise as the difference of two elements from $\mathcal{M}_{K}^{2}$); therefore, we can apply the $\mathcal{M}^{4}_{K}$-RIP constants and Lemmas 1 and 2 to provide the following bounds on both sides (see for details):

The argument is completed by noting that $\delta_{\mathcal{M}_{K}^{2}}\leq\delta_{\mathcal{M}_{K}^{4}}\leq 0.1$. ∎

(Estimation) Let $\Lambda$ be a support corresponding to a subspace in $\mathcal{M}_{K}^{3}$, and define the least squares signal estimate $b$ by $b|_{T}=\Phi_{T}^{\dagger}y$, $b|_{T^{C}}=0$. Then

Since $\Lambda$ is a support corresponding to a subspace in $\mathcal{M}_{K}^{3}$ and $x\in\mathcal{M}_{K}$, we use Lemmas 1 and 2 to obtain

Finally, note that $\delta_{\mathcal{M}_{K}^{3}}\leq\delta_{\mathcal{M}_{K}^{4}}\leq 0.1$. ∎

Proof of Lemma 6: Since $\widehat{x}_{i}$ is the best approximation in $\mathcal{M}_{K}$ to $b$, and $x\in\mathcal{M}_{K}$, we obtain

We use these lemmas in reverse sequence for the inequalities below:

From the recursion on $\widehat{x}_{i}$, we obtain $\|x-\widehat{x}_{i}\|_{2}\leq 2^{-i}\|x\|_{2}+15\|n\|_{2}$. This completes the proof of Theorem 13.∎

Appendix E Proof of Proposition 1

When $K<\log_{2}N$, the number of subtrees of size $K$ of a binary tree of size $N$ is the Catalan number

using Stirling’s approximation. When $K>\log_{2}N$, we partition this count of subtrees into the numbers of subtrees $t_{K,h}$ of size $K$ and height $h$, to obtain

We obtain the following asymptotic identity from [49, page 51]:

We now simplify the formula slightly: we seek a bound for the sum term (which we denote by $\beta_{h}$ for brevity):

Let $m_{\max}=\frac{h}{\pi\sqrt{2K}}$, the value of $m$ for which the term inside the sum (25) is maximum; this is not necessarily an integer. Then,

where the second inequality comes from the fact that the series in the sum is strictly increasing for $m\leq\left\lfloor m_{\max}\right\rfloor$ and strictly decreasing for $m>\left\lceil m_{\max}\right\rceil$. One of the terms in the sum can be added to one of the integrals. If we have that

When the opposite of (26) is true, we have that

Since the term in the sum reaches its maximum for $m_{\max}$, we will have in all three cases that

We perform a change of variables $u=2\pi x$ and define $\sigma=h/\sqrt{2K}$ to obtain

Using the formula for the fourth central moment of a Gaussian distribution:

It is easy to show, using Euler-Maclaurin summations, that

Appendix F Proof of Proposition 3

We wish to find the value of the bound (12) for the subspace count given in (17). We obtain $M\geq\max_{1\leq j\leq\lceil N/K\rceil}M_{j}$, where

We separate the terms that are linear on $K$ and $j$, and obtain

The sequence $\{M_{j}\}_{j=1}^{\left\lceil\frac{N}{K}\right\rceil}$ is a decreasing sequence, since the denominators are decreasing sequences whenever $s>0.5$. We then have

This completes the proof of Proposition 3. ∎

Acknowledgements

We thank Petros Boufounos, Mark Davenport, Yonina Eldar, Moshe Mishali, and Robert Nowak for helpful discussions.

References