CoSaMP: Iterative signal recovery from incomplete and inaccurate samples

Introduction

Most signals of interest contain scant information relative to their ambient dimension, but the classical approach to signal acquisition ignores this fact. We usually collect a complete representation of the target signal and process this representation to sieve out the actionable information. Then we discard the rest. Contemplating this ugly inefficiency, one might ask if it is possible instead to acquire compressive samples. In other words, is there some type of measurement that automatically winnows out the information from a signal? Incredibly, the answer is sometimes yes.

Compressive sampling refers to the idea that, for certain types of signals, a small number of nonadaptive samples carries sufficient information to approximate the signal well. Research in this area has two major components:

How many samples are necessary to reconstruct signals to a specified precision? What type of samples? How can these sampling schemes be implemented in practice?

Given the compressive samples, what algorithms can efficiently construct a signal approximation?

The literature already contains a well-developed theory of sampling, which we summarize below. Although algorithmic work has been progressing, the state of knowledge is less than complete. We assert that a practical signal reconstruction algorithm should have all of the following properties.

It should accept samples from a variety of sampling schemes.

It should succeed using a minimal number of samples.

It should be robust when samples are contaminated with noise.

It should provide optimal error guarantees for every target signal.

It should offer provably efficient resource usage.

To our knowledge, no approach in the literature simultaneously accomplishes all five goals.

This paper presents and analyzes a novel signal reconstruction algorithm that achieves these desiderata. The algorithm is called CoSaMP, from the acrostic Compressive Sampling Matching Pursuit. As the name suggests, the new method is ultimately based on orthogonal matching pursuit (OMP) , but it incorporates several other ideas from the literature to accelerate the algorithm and to provide strong guarantees that OMP cannot. Before we describe the algorithm, let us deliver an introduction to the theory of compressive sampling.

We say that a signal $\bm{x}$ is $s$-sparse when $\left\|{\bm{x}}\right\|_{0}\leq s$. Sparse signals are an idealization that we do not encounter in applications, but real signals are quite often compressible, which means that their entries decay rapidly when sorted by magnitude. As a result, compressible signals are well approximated by sparse signals. We can also talk about signals that are compressible with respect to other orthonormal bases, such as a Fourier or wavelet basis. In this case, the sequence of coefficients in the orthogonal expansion decays quickly. It represents no loss of generality to focus on signals that are compressible with respect to the standard basis, and we do so without regret. For a more precise definition of compressibility, turn to Section 2.6.

The minimum number of measurements $m\geq 2s$ on account of the following simple argument. The sampling matrix must not map two different $s$-sparse signals to the same set of samples. Therefore, each collection of $2s$ columns from the sampling matrix must be nonsingular. It is easy to see that certain Vandermonde matrices satisfy this property, but these matrices are not really suitable for signal acquisition because they contain square minors that are very badly conditioned. As a result, some sparse signals are mapped to very similar sets of samples, and it is unstable to invert the sampling process numerically.

Instead, Candès and Tao proposed the stronger condition that the geometry of sparse signals should be preserved under the action of the sampling matrix . To quantify this idea, they defined the $r$th restricted isometry constant of a matrix $\bm{\Phi}$ as the least number $\delta_{r}$ for which

To acquire $s$-sparse signals, one therefore hopes to achieve a small restricted isometry constant $\delta_{2s}$ with as few samples as possible. A striking fact is that many types of random matrices have excellent restricted isometry behavior. For example, we can often obtain $\delta_{2s}\leq 0.1$ with

measurements, where $\alpha$ is a small integer. Unfortunately, no deterministic sampling matrix is known to satisfy a comparable bound. Even worse, it is computationally difficult to check the inequalities (1.1), so it may never be possible to exhibit an explicit example of a good sampling matrix.

As a result, it is important to understand how random sampling matrices behave. The two quintessential examples are Gaussian matrices and partial Fourier matrices.

If the entries of $\sqrt{m}\bm{\Phi}$ are independent and identically distributed standard normal variables then

If $\sqrt{m}\bm{\Phi}$ is a uniformly random set of $m$ rows drawn from the $N\times N$ unitary discrete Fourier transform (DFT), then

except with probability $N^{-1}$. See for the proof. Experts believe that the power on the first logarithm should be no greater than two .

The Gaussian matrix is important because it has optimal restricted isometry behavior. Indeed, for any $m\times N$ matrix,

on account of profound geometric results of Kashin and Garnaev–Gluskin . Even though partial Fourier matrices may require additional samples to achieve a small restricted isometry constant, they are more interesting for the following reasons .

There are technologies that acquire random Fourier measurements at unit cost per sample.

Other types of sampling matrices, such as the random demodulator , enjoy similar qualities. These traits are essential for the translation of compressive sampling from theory into practice.

2. Signal Recovery Algorithms

The major algorithmic challenge in compressive sampling is to approximate a signal given a vector of noisy samples. The literature describes a huge number of approaches to solving this problem. They fall into three rough categories:

These methods build up an approximation one step at a time by making locally optimal choices at each step. Examples include OMP , stagewise OMP (StOMP) , and regularized OMP (ROMP) .

These techniques solve a convex program whose minimizer is known to approximate the target signal. Many algorithms have been proposed to complete the optimization, including interior-point methods , projected gradient methods , and iterative thresholding .

These methods acquire highly structured samples of the signal that support rapid reconstruction via group testing. This class includes Fourier sampling , chaining pursuit , and HHS pursuit , as well as some algorithms of Cormode–Muthukrishnan and Iwen .

At present, each type of algorithm has its native shortcomings. Many of the combinatorial algorithms are extremely fast—sublinear in the length of the target signal—but they require a large number of somewhat unusual samples that may not be easy to acquire. At the other extreme, convex relaxation algorithms succeed with a very small number of measurements, but they tend to be computationally burdensome. Greedy pursuits—in particular, the ROMP algorithm—are intermediate in their running time and sampling efficiency.

CoSaMP, the algorithm described in this paper, is at heart a greedy pursuit. It also incorporates ideas from the combinatorial algorithms to guarantee speed and to provide rigorous error bounds . The analysis is inspired by the work on ROMP and the work of Candès–Romberg–Tao on convex relaxation methods. In particular, we establish the following result.

We can interpret the error guarantee as follows. In the absence of noise, the algorithm can recover an $s$-sparse signal to arbitrarily high precision. Performance degrades gracefully as the energy in the noise increases. Performance also degrades gracefully for compressible signals. The theorem is ultimately vacuous for signals that cannot be approximated by sparse signals, but compressive sampling is not an appropriate technique for this class.

The running time bound indicates that each matrix--vector multiplication reduces the error by a constant factor (if we amortize over the entire execution). That is, the algorithm has linear convergence Mathematicians sometimes refer to linear convergence as “exponential convergence.”. We find that the total runtime is roughly proportional to (the negation of) the reconstruction signal-to-noise ratio

3. Notation

Finally, we define the pseudoinverse of a tall, full-rank matrix $\bm{A}$ by the formula $\bm{A}^{\dagger}=(\bm{A}^{*}\bm{A})^{-1}\bm{A}^{*}$.

4. Organization

The rest of the paper has the following structure. In Section 2 we introduce the CoSaMP algorithm, we state the major theorems in more detail, and we discuss implementation and resource requirements. Section 3 describes some consequences of the restricted isometry property that pervade our analysis. The central theorem is established for sparse signals in Sections 4 and 5. We extend this result to general signals in Section 6. Finally, Section 7 places the algorithm in the context of previous work. The first appendix presents variations on the algorithm. The second appendix contains a bound on the number of iterations required when the algorithm is implemented using exact arithmetic.

The CoSaMP Algorithm

This section gives an overview of the algorithm, along with explicit pseudocode. It presents the major theorems on the performance of the algorithm. Then it covers details of implementation and bounds on resource requirements.

The most difficult part of signal reconstruction is to identify the locations of the largest components in the target signal. CoSaMP uses an approach inspired by the restricted isometry property. Suppose that the sampling matrix $\bm{\Phi}$ has restricted isometry constant $\delta_{s}\ll 1$. For an $s$-sparse signal $\bm{x}$, the vector $\bm{y}=\bm{\Phi}^{*}\bm{\Phi}\bm{x}$ can serve as a proxy for the signal because the energy in each set of $s$ components of $\bm{y}$ approximates the energy in the corresponding $s$ components of $\bm{x}$. In particular, the largest $s$ entries of the proxy $\bm{y}$ point toward the largest $s$ entries of the signal $\bm{x}$. Since the samples have the form $\bm{u}=\bm{\Phi}\bm{x}$, we can obtain the proxy just by applying the matrix $\bm{\Phi}^{*}$ to the samples.

The algorithm invokes this idea iteratively to approximate the target signal. At each iteration, the current approximation induces a residual, the part of the target signal that has not been approximated. As the algorithm progresses, the samples are updated so that they reflect the current residual. These samples are used to construct a proxy for the residual, which permits us to identify the large components in the residual. This step yields a tentative support for the next approximation. We use the samples to estimate the approximation on this support set using least squares. This process is repeated until we have found the recoverable energy in the signal.

2. Overview

As input, the CoSaMP algorithm requires four pieces of information:

Access to the sampling operator via matrix–vector multiplication.

A vector of (noisy) samples of the unknown signal.

The sparsity of the approximation to be produced.

The algorithm is initialized with a trivial signal approximation, which means that the initial residual equals the unknown target signal. During each iteration, CoSaMP performs five major steps:

Identification. The algorithm forms a proxy of the residual from the current samples and locates the largest components of the proxy.

Support Merger. The set of newly identified components is united with the set of components that appear in the current approximation.

Estimation. The algorithm solves a least-squares problem to approximate the target signal on the merged set of components.

Pruning. The algorithm produces a new approximation by retaining only the largest entries in this least-squares signal approximation.

Sample Update. Finally, the samples are updated so that they reflect the residual, the part of the signal that has not been approximated.

These steps are repeated until the halting criterion is triggered. In the body of this work, we concentrate on methods that use a fixed number of iterations. Appendix A discusses some other simple stopping rules that may also be useful in practice.

Pseudocode for CoSaMP appears as Algorithm 2.1. This code describes the version of the algorithm that we analyze in this paper. Nevertheless, there are several adjustable parameters that may improve performance: the number of components selected in the identification step and the number of components retained in the pruning step. For a brief discussion of other variations on the algorithm, turn to Appendix A.

3. Performance Guarantees

This section describes our theoretical analysis of the behavior of CoSaMP. The next section covers the resource requirements of the algorithm. Afterward, Section 2.5 combines these materials to establish Theorem A.

Our results depend on a set of hypotheses that has become common in the compressive sampling literature. Let us frame the standing assumptions:

We also define the unrecoverable energy $\nu$ in the signal. This quantity measures the baseline error in our approximation that occurs because of noise in the samples or because the signal is not sparse.

$\nu=\left\|{\bm{x}-\bm{x}_{s}}\right\|_{2}+\frac{1}{\sqrt{s}}\left\|{\bm{x}-\bm{x}_{s}}\right\|_{1}+\left\|{\bm{e}}\right\|_{2}.$ (2.1)

We postpone a more detailed discussion of the unrecoverable energy until Section 2.6.

Our key result is that CoSaMP makes significant progress during each iteration where the approximation error is large relative to the unrecoverable energy.

For each iteration $k\geq 0$, the signal approximation $\bm{a}^{k}$ is $s$-sparse and

The proof of Theorem 2.1 will occupy us for most of this paper. In Section 4, we establish an analog for sparse signals. The version for general signals appears as a corollary in Section 6.

Theorem 2.1 has some immediate consequences for the quality of reconstruction with respect to standard signal metrics. In this setting, a sensible definition of the signal-to-noise ratio (SNR) is

Both quantities are measured in decibels. Theorem 2.1 implies that, after $k$ iterations, the reconstruction SNR satisfies

In words, each iteration reduces the reconstruction SNR by about 3 decibels until the error nears the noise floor. To reduce the error to its minimal value, the number of iterations is proportional to the SNR.

Let us consider a slightly different scenario. Suppose that the signal $\bm{x}$ is $s$-sparse, so the unrecoverable energy $\nu=\left\|{\bm{e}}\right\|_{2}$. Define the dynamic range

Assume moreover that the minimum nonzero component of the signal is at least $40\nu$. Using the fact that $\left\|{\bm{x}}\right\|_{2}\leq\sqrt{s}\left\|{\bm{x}}\right\|_{\infty}$, it is easy to check that $\left\|{\bm{x}-\bm{a}}\right\|_{2}\leq\min\left|{x_{i}}\right|$ as soon as the number $k$ of iterations satisfies

It follows that the support of the approximation $\bm{a}$ must contain every entry in the support of the signal $\bm{x}$.

This discussion suggests that the number of iterations might be substantial if we require a very low reconstruction SNR or if the signal has a very wide dynamic range. This initial impression is not entirely accurate. We have established that, when the algorithm performs arithmetic to high enough precision, then a fixed number of iterations suffices to reduce the approximation error to the same order as the unrecoverable energy. Here is a result for exact computations.

Suppose that CoSaMP is implemented with exact arithmetic. After at most $6(s+1)$ iterations, CoSaMP produces an $s$-sparse approximation $\bm{a}$ that satisfies

If we solve the least-squares problems to high precision, an analogous result holds, but the approximation guarantee contains an extra term that comes from solving the least-squares problems imperfectly. In practice, it may be more efficient overall to solve the least-squares problems to low precision. The correct amount of care seems to depend on the relative costs of forming the signal proxy and solving the least-squares problem, which are the two most expensive steps in the algorithm. We discuss this point in the next section. Ultimately, the question is best settled with empirical studies.

In the hypotheses, a bound on the restricted isometry constant $\delta_{2s}$ also suffices. Indeed, Corollary 3.4 of the sequel implies that $\delta_{4s}\leq 0.1$ holds whenever $\delta_{2s}\leq 0.025$.

Choosing $\bm{y}=\bm{x}-\bm{x}_{s/2}$ and $t=s/2$, we reach

4. Implementation and Resource Requirements

CoSaMP was designed to be a practical method for signal recovery. An efficient implementation of the algorithm requires some ideas from numerical linear algebra, as well as some basic techniques from the theory of algorithms. This section discusses the key issues and develops an analysis of the running time for the two most common scenarios.

We focus on the least-squares problem in the estimation step because it is the major obstacle to a fast implementation of the algorithm. The algorithm guarantees that the matrix $\bm{\Phi}_{T}$ never has more than $3s$ columns, so our assumption $\delta_{4s}\leq 0.1$ implies that the matrix $\bm{\Phi}_{T}$ is extremely well conditioned. As a result, we can apply the pseudoinverse $\bm{\Phi}_{T}^{\dagger}=(\bm{\Phi}_{T}^{*}\bm{\Phi}_{T})^{-1}\bm{\Phi}_{T}^{*}$ very quickly using an iterative method, such as Richardson’s iteration [1, Sec. 7.2.3] or conjugate gradient [1, Sec. 7.4]. These techniques have the additional advantage that they only interact with the matrix $\bm{\Phi}_{T}$ through its action on vectors. It follows that the algorithm performs better when the sampling matrix has a fast matrix–vector multiply.

The remaining steps of the algorithm are standard. Let us estimate the operation counts.

Forming the proxy is dominated by the cost of the matrix–vector multiply $\bm{\Phi}^{*}\bm{v}$.

We use Richardson’s iteration or conjugate gradient to compute $\bm{\Phi}_{T}^{\dagger}\bm{u}$. Initializing the least-squares algorithm requires a matrix–vector multiply with $\bm{\Phi}_{T}^{*}$. Each iteration of the least-squares method requires one matrix–vector multiply each with $\bm{\Phi}_{T}$ and $\bm{\Phi}_{T}^{*}$. Since $\bm{\Phi}_{T}$ is a submatrix of $\bm{\Phi}$, the matrix–vector multiplies can also be obtained from multiplication with the full matrix. We prove in Section 5 that a constant number of least-squares iterations suffices for Theorem 2.1 to hold.

This step is dominated by the cost of the multiplication of $\bm{\Phi}$ with the $s$-sparse vector $\bm{a}^{k}$.

The following result summarizes this discussion.

5. Proof of Theorem A

In the statement of the theorem, perform the substitution $s/2\mapsto s$. Finally, we replace $\delta_{8s}$ with $\delta_{2s}$ by means of Corollary 3.4, which states that $\delta_{cr}\leq c\cdot\delta_{2r}$ for any positive integers $c$ and $r$.

6. The Unrecoverable Energy

Since the unrecoverable energy $\nu$ plays a central role in our analysis of CoSaMP, it merits some additional discussion. In particular, it is informative to examine the unrecoverable energy in a compressible signal. Let $p$ be a number in the interval $(0,1)$. We say that $\bm{x}$ is $p$-compressible with magnitude $R$ if the sorted components of the signal decay at the rate

When $p$ is small, the first term in the unrecoverable energy decays rapidly as the sparsity level $s$ increases. For the class of $p$-compressible signals, the bound (2.3) on the unrecoverable energy is sharp, modulo the exact values of the constants.

With these inequalities, we can see that CoSaMP recovers compressible signals efficiently. Let us calculate the number of iterations required to reduce the approximation error from $\left\|{\bm{x}}\right\|_{2}$ to the optimal level (2.3). For compressible signals, the energy $\left\|{\bm{x}}\right\|_{2}\leq 2R$, so

Restricted Isometry Consequences

When the sampling matrix satisfies the restricted isometry inequalities (1.1), it has several other properties that we require repeatedly in the proof that the CoSaMP algorithm is correct. Our first observation is a simple translation of (1.1) into other terms.

Suppose $\bm{\Phi}$ has restricted isometry constant $\delta_{r}$. Let $T$ be a set of $r$ indices or fewer. Then

where the last two statements contain an upper and lower bound, depending on the sign chosen.

The restricted isometry inequalities (1.1) imply that the singular values of $\bm{\Phi}_{T}$ lie between $\sqrt{1-\delta_{r}}$ and $\sqrt{1+\delta_{r}}$. The bounds follow from standard relationships between the singular values of $\bm{\Phi}_{T}$ and the singular values of basic functions of $\bm{\Phi}_{T}$. ∎

A second consequence is that disjoint sets of columns from the sampling matrix span nearly orthogonal subspaces. The following result quantifies this observation.

Suppose $\bm{\Phi}$ has restricted isometry constant $\delta_{r}$. Let $S$ and $T$ be disjoint sets of indices whose combined cardinality does not exceed $r$. Then

Abbreviate $R=S\cup T$, and observe that $\bm{\Phi}_{S}^{*}\bm{\Phi}_{T}$ is a submatrix of $\bm{\Phi}_{R}^{*}\bm{\Phi}_{R}-\mathbf{I}$. The spectral norm of a submatrix never exceeds the norm of the entire matrix. We discern that

because the eigenvalues of $\bm{\Phi}_{R}^{*}\bm{\Phi}_{R}$ lie between $1-\delta_{r}$ and $1+\delta_{r}$. ∎

This result will be applied through the following corollary.

Suppose $\bm{\Phi}$ has restricted isometry constant $\delta_{r}$. Let $T$ be a set of indices, and let $\bm{x}$ be a vector. Provided that $r\geq\left|{T\cup\operatorname{supp}(\bm{x})}\right|$,

Define $S=\operatorname{supp}(\bm{x})\setminus T$, so we have $\bm{x}|_{S}=\bm{x}|_{T^{c}}$. Thus,

As a second corollary, we show that $\delta_{2r}$ gives weak control over the higher restricted isometry constants.

Let $c$ and $r$ be positive integers. Then $\delta_{cr}\leq c\cdot\delta_{2r}$.

The result is clearly true for $c=1,2,$ so we assume $c\geq 3$. Let $S$ be an arbitrary index set of size $cr$, and let $\bm{M}=\bm{\Phi}_{S}^{*}\bm{\Phi}_{S}-\mathbf{I}$. It suffices to check that $\left\|{\bm{M}}\right\|\leq c\cdot\delta_{2r}$. To that end, we break the matrix $\bm{M}$ into $r\times r$ blocks, which we denote $\bm{M}_{ij}$. A block version of Gershgorin’s theorem states that $\left\|{\bm{M}}\right\|$ satisfies at least one of the inequalities

The derivation is entirely analogous with the usual proof of Gershgorin’s theorem, so we omit the details. For each diagonal block, we have $\left\|{\bm{M}_{ii}}\right\|\leq\delta_{r}$ because of the restricted isometry inequalities (1.1). For each off-diagonal block, we have $\left\|{\bm{M}_{ij}}\right\|\leq\delta_{2r}$ because of Proposition 3.2. Substitute these bounds into the block Gershgorin theorem and rearrange to complete the proof. ∎

Finally, we present a result that measures how much the sampling matrix inflates nonsparse vectors. This bound permits us to establish the major results for sparse signals and then transfer the conclusions to the general case.

Suppose that $\bm{\Phi}$ verifies the upper inequality of (1.1), viz.

We repeat the geometric argument of Rudelson that is presented in .

and notice that, by hypothesis, the operator norm

The content of the proposition is the claim that

To establish this point, it suffices to check that $K\subset S$.

Instead, we prove the reverse inclusion for the polars: $S^{\circ}\subset K^{\circ}$. The norm with unit ball $S^{\circ}$ is calculated as

Consider a vector $\bm{u}$ in the unit ball $S^{\circ}$, and let $I$ be a set of $r$ coordinates where $\bm{u}$ is largest in magnitude. We must have

or else $\left|{u_{i}}\right|>\frac{1}{\sqrt{r}}$ for each $i\in I$. But then $\left\|{\bm{u}}\right\|_{S^{\circ}}\geq\left\|{\bm{u}|_{I}}\right\|_{2}>1$, a contradiction. Therefore, we may write

In summary, $S^{\circ}\subset K^{\circ}$. ∎

The Iteration Invariant: Sparse Case

We now commence the proof of Theorem 2.1. For the moment, let us assume that the signal is actually sparse. Section 6 removes this assumption.

The result states that each iteration of the algorithm reduces the approximation error by a constant factor, while adding a small multiple of the noise. As a consequence, when the approximation error is large in comparison with the noise, the algorithm makes substantial progress in identifying the unknown signal.

Assume that $\bm{x}$ is $s$-sparse. For each $k\geq 0$, the signal approximation $\bm{a}^{k}$ is $s$-sparse, and

The argument proceeds in a sequence of short lemmas, each corresponding to one step in the algorithm. Throughout this section, we retain the assumption that $\bm{x}$ is $s$-sparse.

Fix an iteration $k\geq 1$. We write $\bm{a}=\bm{a}^{k-1}$ for the signal approximation at the beginning of the iteration. Define the residual $\bm{r}=\bm{x}-\bm{a}$, which we interpret as the part of the signal we have not yet recovered. Since the approximation $\bm{a}$ is always $s$-sparse, the residual $\bm{r}$ must be $2s$-sparse. Notice that the vector $\bm{v}$ of updated samples can be viewed as noisy samples of the residual:

2. Identification

The identification phase produces a set of components where the residual signal still has a lot of energy.

The set $\Omega=\operatorname{supp}(\bm{y}_{2s})$ contains at most $2s$ indices, and

The identification phase forms a proxy $\bm{y}=\bm{\Phi}^{*}\bm{v}$ for the residual signal. The algorithm then selects a set $\Omega$ of $2s$ components from $\bm{y}$ that have the largest magnitudes. The goal of the proof is to show that the energy in the residual on the set $\Omega^{c}$ is small in comparison with the total energy in the residual.

Define the set $R=\operatorname{supp}(\bm{r})$. Since $R$ contains at most $2s$ elements, our choice of $\Omega$ ensures that $\left\|{\bm{y}|_{R}}\right\|_{2}\leq\left\|{\bm{y}|_{\Omega}}\right\|_{2}$. By squaring this inequality and canceling the terms in $R\cap\Omega$, we discover that

Since the coordinate subsets here contain few elements, we can use the restricted isometry constants to provide bounds on both sides.

First, observe that the set $\Omega\setminus R$ contains at most $2s$ elements. Therefore, we may apply Proposition 3.1 and Corollary 3.3 to obtain

Likewise, the set $R\setminus\Omega$ contains $2s$ elements or fewer, so Proposition 3.1 and Corollary 3.3 yield

Since the residual is supported on $R$, we can rewrite $\bm{r}|_{R\setminus\Omega}=\bm{r}|_{\Omega^{c}}$. Finally, combine the last three inequalities and rearrange to obtain

Invoke the numerical hypothesis that $\delta_{2s}\leq\delta_{4s}\leq 0.1$ to complete the argument. ∎

3. Support Merger

The next step of the algorithm merges the support of the current signal approximation $\bm{a}$ with the newly identified set of components. The following result shows that components of the signal $\bm{x}$ outside this set have very little energy.

Let $\Omega$ be a set of at most $2s$ indices. The set $T=\Omega\cup\operatorname{supp}(\bm{a})$ contains at most $3s$ indices, and

Since $\operatorname{supp}(\bm{a})\subset T$, we find that

where the inequality follows from the containment $T^{c}\subset\Omega^{c}$. ∎

4. Estimation

The estimation step of the algorithm solves a least-squares problem to obtain values for the coefficients in the set $T$. We need a bound on the error of this approximation.

Let $T$ be a set of at most $3s$ indices, and define the least-squares signal estimate $\bm{b}$ by the formulae

This result assumes that we solve the least-squares problem in infinite precision. In practice, the right-hand side of the bound contains an extra term owing to the error from the iterative least-squares solver. In Section 5, we study how many iterations of the least-squares solver are required to make the least-squares error negligible in the present argument.

Using the expression $\bm{u}=\bm{\Phi}\bm{x}+\bm{e}$ and the fact $\bm{\Phi}_{T}^{\dagger}\bm{\Phi}_{T}=\mathbf{I}_{T}$, we calculate that

The cardinality of $T$ is at most $3s$, and $\bm{x}$ is $s$-sparse, so Proposition 3.1 and Corollary 3.3 imply that

Finally, invoke the hypothesis that $\delta_{3s}\leq\delta_{4s}\leq 0.1$. ∎

5. Pruning

The final step of each iteration is to prune the intermediate approximation to its largest $s$ terms. The following lemma provides a bound on the error in the pruned approximation.

The pruned approximation $\bm{b}_{s}$ satisfies

The intuition is that $\bm{b}_{s}$ is close to $\bm{b}$, which is close to $\bm{x}$. Rigorously,

The second inequality holds because $\bm{b}_{s}$ is the best $s$-sparse approximation to $\bm{b}$. In particular, the $s$-sparse vector $\bm{x}$ is a worse approximation. ∎

6. Proof of Theorem 4.1

We now complete the proof of the iteration invariant for sparse signals, Theorem 4.1. At the end of an iteration, the algorithm forms a new approximation $\bm{a}^{k}=\bm{b}_{s}$, which is evidently $s$-sparse. Applying the lemmas we have established, we easily bound the error:

To obtain the second bound in Theorem 4.1, simply solve the error recursion and note that

Analysis of Iterative Least-squares

To develop an efficient implementation of CoSaMP, it is critical to use an iterative method when we solve the least-squares problem in the estimation step. The two natural choices are Richardson’s iteration and conjugate gradient. The efficacy of these methods rests on the assumption that the sampling operator has small restricted isometry constants. Indeed, since the set $T$ constructed in the support merger step contains at most $3s$ components, the hypothesis $\delta_{4s}\leq 0.1$ ensures that the condition number

This condition number is closely connected with the performance of Richardson’s iteration and conjugate gradient. In this section, we show that Theorem 4.1 holds if we perform a constant number of iterations of either least-squares algorithm.

For completeness, let us explain how Richardson’s iteration can be applied to solve the least-squares problems that arise in CoSaMP. Suppose we wish to compute $\bm{A}^{\dagger}\bm{u}$ where $\bm{A}$ is a tall, full-rank matrix. Recalling the definition of the pseudoinverse, we realize that this amounts to solving a linear system of the form

This problem can be approached by splitting the Gram matrix:

where $\bm{M}=\bm{A}^{*}\bm{A}-\mathbf{I}$. Given an initial iterate $\bm{z}^{0}$, Richardon’s method produces the subsequent iterates via the formula

Evidently, this iteration requires only matrix–vector multiplies with $\bm{A}$ and $\bm{A}^{*}$. It is worth noting that Richardson’s method can be accelerated [1, Sec. 7.2.5], but we omit the details.

It is quite easy to analyze Richardson’s iteration [1, Sec. 7.2.1]. Observe that

In words, the iteration converges linearly.

In our setting, $\bm{A}=\bm{\Phi}_{T}$ where $T$ is a set of at most $3s$ indices. Therefore, the restricted isometry inequalities (1.1) imply that

We have assumed that $\delta_{3s}\leq\delta_{4s}\leq 0.1$, which means that the iteration converges quite fast. Once again, the restricted isometry behavior of the sampling matrix plays an essential role in the performance of the CoSaMP algorithm.

Conjugate gradient provides even better guarantees for solving the least-squares problem, but it is somewhat more complicated to describe and rather more difficult to analyze. We refer the reader to [1, Sec. 7.4] for more information. The following lemma summarizes the behavior of both Richardson’s iteration and conjugate gradient in our setting.

Conjugate gradient produces a sequence of iterates that satisfy

2. Initialization

Iterative least-squares algorithms must be seeded with an initial iterate, and their performance depends heavily on a wise selection thereof. CoSaMP offers a natural choice for the initializer: the current signal approximation. As the algorithm progresses, the current signal approximation provides an increasingly good starting point for solving the least-squares problem.

Let $\bm{x}$ be an $s$-sparse signal with noisy samples $\bm{u}=\bm{\Phi}\bm{x}+\bm{e}$. Let $\bm{a}^{k-1}$ be the signal approximation at the end of the $(k-1)$th iteration, and let $T$ be the set of components identified by the support merger. Then

By construction of $T$, the approximation $\bm{a}^{k-1}$ is supported inside $T$, so

Using Lemma 4.4, we may calculate how far $\bm{a}^{k-1}$ lies from the solution to the least-squares problem.

Roughly, the error in the initial iterate is controlled by the current approximation error. ∎

3. Iteration Count

We need to determine how many iterations of the least-squares algorithm are required to ensure that the approximation produced is sufficiently good to support the performance of CoSaMP.

Suppose that we initialize the LS algorithm with $\bm{z}^{0}=\bm{a}^{k-1}$. After at most three iterations, both Richardson’s iteration and conjugate gradient produce a signal estimate $\bm{b}$ that satisfies

Combine Lemma 5.1 and Lemma 5.2 to see that three iterations of Richardson’s method yield

The bound for conjugate gradient is slightly better. Let $\bm{b}|_{T}=\bm{z}^{3}$. According to the estimation result, Lemma 4.4, we have

An application of the triangle inequality completes the argument. ∎

4. CoSaMP with Iterative least-squares

Finally, we need to check that the sparse iteration invariant, Theorem 4.1 still holds when we use an iterative least-squares algorithm.

Suppose that we use Richardson’s iteration or conjugate gradient for the estimation step, initializing the LS algorithm with the current approximation $\bm{a}^{k-1}$ and performing three LS iterations. Then Theorem 4.1 still holds.

We repeat the calculation in Section 4.6 using Corollary 5.3 instead of the simple estimation lemma. To that end, recall that the residual $\bm{r}=\bm{x}-\bm{a}^{k-1}$. Then

This bound is precisely what is required for the theorem to hold. ∎

Extension to General Signals

In this section, we finally complete the proof of the main result for CoSaMP, Theorem 2.1. The remaining challenge is to remove the hypothesis that the target signal is sparse, which we framed in Theorems 4.1 and 5.4. Although this difficulty might seem large, the solution is simple and elegant. It turns out that we can view the noisy samples of a general signal as samples of a sparse signal contaminated with a different noise vector that implicitly reflects the tail of the original signal.

Decompose $\bm{x}=\bm{x}_{s}+(\bm{x}-\bm{x}_{s})$ to obtain $\bm{u}=\bm{\Phi}\bm{x}_{s}+\widetilde{\bm{e}}$ where $\widetilde{\bm{e}}=\bm{\Phi}(\bm{x}-\bm{x}_{s})+\bm{e}$. To compute the size of the error term, we simply apply the triangle inequality and Proposition 3.5:

Finally, invoke the fact that $\delta_{s}\leq\delta_{4s}\leq 0.1$ to obtain $\sqrt{1+\delta_{s}}\leq 1.05$. ∎

This lemma is just the tool we require to complete the remaining argument.

Let $\bm{x}$ be a general signal, and use Lemma 6.1 to write the noisy vector of samples $\bm{u}=\bm{\Phi}\bm{x}_{s}+\widetilde{\bm{e}}$. Apply the sparse iteration invariant, Theorem 4.1, or the analog for iterative least-squares, Theorem 5.4. We obtain

Invoke the lower and upper triangle inequalities to obtain

Finally, recall the estimate for $\left\|{\widetilde{\bm{e}}}\right\|_{2}$ from Lemma 6.1, and simplify to reach

where $\nu$ is the unrecoverable energy (2.1). ∎

Discussion and Related Work

CoSaMP draws on both algorithmic ideas and analytic techniques that have appeared before. This section describes the other major signal recovery algorithms, and it compares them with CoSaMP. It also attempts to trace the key ideas in the algorithm back to their sources.

We begin with a short discussion of the major algorithmic approaches to signal recovery from compressive samples. We focus on provably correct methods, although we acknowledge that some ad hoc techniques provide excellent empirical results.

The initial discovery works on compressive sampling proposed to perform signal recovery by solving a convex optimization problem . Given a sampling matrix $\bm{\Phi}$ and a noisy vector of samples $\bm{u}=\bm{\Phi}\bm{x}+\bm{e}$, consider the mathematical program

provided that the sampling matrix $\bm{\Phi}$ has restricted isometry constant $\delta_{4s}\leq 0.2$. In , the hypothesis on the restricted isometry constant is sharpened to $\delta_{2s}\leq\sqrt{2}-1$. The error bound for CoSaMP is equivalent, modulo the exact value of the constants.

Tropp and Gilbert proposed the use of a greedy iterative algorithm called orthogonal matching pursuit (OMP) for signal recovery . The algorithm initializes the current sample vector $\bm{v}=\bm{u}$. In each iteration, it forms the signal proxy $\bm{y}=\bm{\Phi}^{*}\bm{v}$ and identifies a component of the proxy with largest magnitude. It adds the new component to the set $T$ of previously identified components. Then OMP forms a new signal approximation by solving a least-squares problem: $\bm{a}=\bm{\Phi}_{T}^{\dagger}\bm{u}$. Finally, it updates the samples $\bm{v}=\bm{u}-\bm{\Phi}\bm{a}$. These steps are repeated until a halting criterion is satisfied.

Donoho et al. invented another greedy iterative method called stagewise OMP, or StOMP . This algorithm uses the signal proxy to select multiple components at each step, using a rule inspired by ideas from wireless communications. The algorithm is faster than OMP because of the selection rule, and it sometimes provides good performance, although parameter tuning can be difficult. There are no rigorous results available for StOMP.

Very recently, Needell and Vershynin developed and analyzed another greedy approach, called regularized OMP, or ROMP . This algorithm is similar to OMP but uses a more sophisticated selection rule. Among the $s$ largest entries of the signal proxy, it identifies the largest subset whose entries differ in magnitude by at most a factor of two. The work on ROMP represents an advance because the authors establish under restricted isometry hypotheses that their algorithm can approximately recover any compressible signal from noisy samples. More precisely, suppose that the sampling matrix $\bm{\Phi}$ has restricted isometry constant $\delta_{8s}\leq 0.01/\sqrt{\log s}$. Given noisy samples $\bm{u}=\bm{\Phi}\bm{x}+\bm{e}$, ROMP produces a $2s$-sparse signal approximation $\bm{a}$ that satisfies

This result is comparable with the result for convex relaxation, aside from the extra logarithmic factor in the restricted isometry hypothesis and the error bound. The results for CoSaMP show that it does not suffer these parasitic factors, so its performance is essentially optimal.

After we initially presented this work, Dai and Milenkovic developed an algorithm called Subspace Pursuit that is very similar to CoSaMP. They established that their algorithm offers performance guarantees analogous with those for CoSaMP. See for details.

Finally, we note that there is a class of sublinear algorithms for signal reconstruction from compressive samples. A sublinear algorithm uses time and space resources that are asymptotically smaller than the length of the signal. One of the earliest such techniques is the Fourier sampling algorithm of Gilbert et al. . This algorithm uses random (but structured) time samples to recover signals that are compressible with respect to the discrete Fourier basis. Given $s\polylog(N)$ samples The term $\polylog$ indicates a function that is dominated by a polynomial in the logarithm of its argument., Fourier sampling produces a signal approximation $\bm{a}$ that satisfies

except with probability $N^{-1}$. The result for Fourier sampling holds for each signal (rather than for all). Later, Gilbert et al. developed two other sublinear algorithms, chaining pursuit and HHS pursuit , that offer uniform guarantees for all signals. Chaining pursuit has an error bound

which is somewhat worse than (7.2). HHS pursuit achieves the error bound (7.2). These methods all require more measurements than the linear and superlinear algorithms (by logarithmic factors), and these measurements must be highly structured. As a result, the sublinear algorithms may not be useful in practice.

The sublinear algorithms are all combinatorial in nature. They use ideas from group testing to identify the support of the signal quickly. There are several other combinatorial signal recovery methods due to Cormode–Muthukrishnan and Iwen . These algorithms have drawbacks similar to the sublinear approaches.

2. Relative Performance

Table 2 summarizes the relative behavior of these algorithms in terms of the following criteria.

Does the algorithm work for a variety of sampling schemes? Or does it require structured samples? The designation “RIP” means that a bound on a restricted isometry constant suffices. “Subgauss.” means that the algorithm succeeds for the class of subgaussian sampling matrices.

Does the algorithm recover all signals given a fixed sampling matrix? Or do the results require a sampling matrix to be drawn at random for each signal?

Does the algorithm succeed when (a) the signal is compressible but not sparse and (b) when the samples are contaminated with noise? In most cases, stable algorithms have error bounds similar to (7.2). See the discussion above for details.

Of the linear and superlinear algorithms, CoSaMP achieves the best performance on all these metrics. Although CoSaMP is slower than the sublinear algorithms, it makes up for this shortcoming by allowing more general sampling matrices and requiring fewer samples.

3. Key Ideas

We conclude with a historical overview of the ideas that inform the CoSaMP algorithm and its analysis.

The overall greedy iterative structure of CoSaMP has a long history. The idea of approaching sparse approximation problems in this manner dates to the earliest algorithms. In particular, methods for variable selection in regression, such as forward selection and its relatives, all take this form . Temlyakov’s survey describes the historical role of greedy algorithms in nonlinear approximation. Mallat and Zhang introduced greedy algorithms into the signal processing literature and proposed the name matching pursuit . Gilbert, Strauss, and their collaborators showed how to incorporate greedy iterative strategies into fast algorithms for sparse approximation problems, and they established the first rigorous guarantees for greedy methods . Tropp provided a new theoretical analysis of OMP in his work . Subsequently, Tropp and Gilbert proved that OMP was effective for compressive sampling .

Unlike the simplest greedy algorithms, CoSaMP identifies many components during each iteration, which allows the algorithm to run faster for many types of signals. It is not entirely clear where this idea first appeared. Several early algorithms of Gilbert et al. incorporate this approach , and it is an essential feature of the Fourier sampling algorithm . More recent compressive sampling recovery algorithms also select multiple indices, including chaining pursuit , HHS pursuit , StOMP , and ROMP .

CoSaMP uses the restricted isometry properties of the sampling matrix to ensure that the identification step is successful. Candès and Tao isolated the restricted isometry conditions in their work on convex relaxation methods for compressive sampling . The observation that restricted isometries can also be used to ensure the success of greedy methods is relatively new. This idea plays a role in HHS pursuit , but it is expressed more completely in the analysis of ROMP .

The pruning step of CoSaMP is essential to maintain the sparsity of the approximation, which is what permits us to use restricted isometries in the analysis of the algorithm. It also has significant ramifications for the running time because it impacts the speed of the iterative least-squares algorithms. This technique originally appeared in HHS pursuit .

The iteration invariant, Theorem 2.1, states that if the error is large then CoSaMP makes substantial progress. This approach to the overall analysis echoes the analysis of other greedy iterative algorithms, including the Fourier sampling method and HHS Pursuit .

Finally, mixed-norm error bounds, such as that in Theorem A, have become an important feature of the compressive sampling literature. This idea appears in the work of Candès–Romberg–Tao on convex relaxation ; it is used in the analysis of HHS pursuit ; it also plays a role in the theoretical treatment of Cohen–Dahmen–DeVore .

Acknowledgment

We would like to thank Martin Strauss for many inspiring discussions. He is ultimately responsible for many of the ideas in the algorithm and analysis. We would also like to thank Roman Vershynin for suggestions that drastically simplified the proofs.

Appendix A Algorithmic Variations

This appendix describes other possible halting criteria and their consequences. It also proposes some other variations on the algorithm.

There are three natural approaches to halting the algorithm. The first, which we have discussed in the body of the paper, is to stop after a fixed number of iterations. Another possibility is to use the norm $\left\|{\bm{v}}\right\|_{2}$ of the current samples as evidence about the norm $\left\|{\bm{r}}\right\|_{2}$ of the residual. A third possibility is to use the magnitude $\left\|{\bm{y}}\right\|_{\infty}$ of the entries of the proxy to bound the magnitude $\left\|{\bm{r}}\right\|_{\infty}$ of the entries of the residual.

It suffices to discuss halting criteria for sparse signals because Lemma 6.1 shows that the general case can be viewed in terms of sampling a sparse signal. Let $\bm{x}$ be an $s$-sparse signal, and let $\bm{a}$ be an $s$-sparse approximation. The residual $\bm{r}=\bm{x}-\bm{a}$. We write $\bm{v}=\bm{\Phi}\bm{r}+\bm{e}$ for the induced noisy samples of the residual and $\bm{y}=\bm{\Phi}^{*}\bm{v}$ for the signal proxy.

The discussion proceeds in two steps. First, we argue that an a priori halting criterion will result in a guarantee about the quality of the final signal approximation.

The halting criterion $\left\|{\bm{v}}\right\|_{2}\leq\varepsilon$ ensures that

The halting criterion $\left\|{\bm{y}}\right\|_{\infty}\leq\eta/\sqrt{2s}$ ensures that

Since $\bm{r}$ is $2s$-sparse, Proposition 3.1 ensures that

If $\left\|{\bm{v}}\right\|_{2}\leq\varepsilon$, it is immediate that

The definition $\bm{r}=\bm{x}-\bm{a}$ and the numerical bound $\delta_{2s}\leq\delta_{4s}\leq 0.1$ dispatch the first claim.

Let $R=\operatorname{supp}(\bm{r})$, and note that $\left|{R}\right|\leq 2s$. Proposition 3.1 results in

we find that the requirement $\left\|{\bm{y}}\right\|_{\infty}\leq\eta/\sqrt{2s}$ ensures that

The numerical bound $\delta_{2s}\leq 0.1$ completes the proof.

Second, we check that each halting criterion is triggered when the residual has the desired property.

The halting criterion $\left\|{\bm{v}}\right\|_{2}\leq\varepsilon$ is triggered as soon as

The halting criterion $\left\|{\bm{y}}\right\|_{\infty}\leq\eta/\sqrt{2s}$ is triggered as soon as

ensures that $\left\|{\bm{v}}\right\|_{2}\leq\varepsilon$. Note that $\delta_{2s}\leq 0.1$ to complete the first part of the argument.

Now let $R$ be the singleton containing the index of a largest-magnitude coefficient of $\bm{y}$. Proposition 3.1 implies that

By the first part of this theorem, the halting criterion $\left\|{\bm{y}}\right\|_{\infty}\leq\eta/\sqrt{2s}$ is triggered as soon as

Since $\bm{x}-\bm{a}$ is $2s$-sparse, we have the bound $\left\|{\bm{x}-\bm{a}}\right\|_{2}\leq\sqrt{2s}\left\|{\bm{x-a}}\right\|_{\infty}$. To wrap up, recall that $\delta_{1}\leq\delta_{2s}\leq 0.1$. ∎

A.2. Other Variations

This section briefly describes several natural variations on CoSaMP that may improve its performance.

Here is a version of the algorithm that is, perhaps, simpler than Algorithm 2.1. At each iteration, we approximate the current residual rather than the entire signal. This approach is similar to HHS Pursuit . The inner loop changes in the following manner.

As before, select $\Omega=\operatorname{supp}(\bm{y}_{2s})$.

Solve a least-squares problem with the current samples instead of the original samples to obtain an approximation of the residual signal. Formally, $\bm{b}=\bm{\Phi}_{\Omega}^{\dagger}\bm{v}$. In this case, one initializes the iterative least-squares algorithm with the zero vector to take advantage of the fact that the residual is becoming small.

Add this approximation of the residual to the previous approximation of the signal to obtain a new approximation of the signal: $\bm{c}=\bm{a}^{k-1}+\bm{b}$.

Construct the $s$-sparse signal approximation: $\bm{a}^{k}=\bm{c}_{s}$.

Update the samples as before: $\bm{v}=\bm{u}-\bm{\Phi}\bm{a}^{k}$.

We can show that this algorithm satisfies a result similar to Theorem 2.1 by adapting the current argument. We were unable to verify an analog of Theorem 2.2. Nevertheless, we believe that this version of the algorithm is also promising.

After the inner loop of the algorithm is complete, we can solve another least-squares problem in an effort to improve the final result. If $\bm{a}$ is the approximation at the end of the loop, we set $T=\operatorname{supp}(\bm{a})$. Then solve $\bm{b}=\bm{\Phi}_{T}^{\dagger}\bm{u}$ and output the $s$-sparse signal approximation $\bm{b}$. Note that the output approximation is not guaranteed to be better than $\bm{a}$ because of the noise vector $\bm{e}$, but it should never be much worse.

Another variation is to prune the merged support $T$ down to $s$ entries before solving the least-squares problem. One may use the values of the proxy $\bm{y}$ as surrogates for the unknown values of the new approximation on the set $\Omega$. Since the least-squares problem is solved at the end of the iteration, the columns of $\bm{\Phi}$ that are used in the least-squares approximation are orthogonal to the current samples $\bm{v}$. As a result, the identification step always selects new components in each iteration. We have not attempted an analysis of this algorithm.

Appendix B Iteration Count

In this appendix, we obtain an estimate on the number of iterations of the CoSaMP algorithm necessary to identify the recoverable energy in a sparse signal, assuming exact arithmetic. Except where stated explicitly, we assume that $\bm{x}$ is $s$-sparse. It turns out that the number of iterations depends heavily on the signal structure. Let us explain the intuition behind this fact.

When the entries in the signal decay rapidly, the algorithm must identify and remove the largest remaining entry from the residual before it can make further progress on the smaller entries. Indeed, the large component in the residual contaminates each component of the signal proxy. In this case, the algorithm may require an iteration or more to find each component in the signal.

To quantify these intuitions, we want to collect the components of the signal into groups that are comparable with each other. To that end, define the component bands of a signal $\bm{x}$ by the formulae

The profile of the signal is the number of bands that are nonempty:

In words, the profile counts how many orders of magnitude at which the signal has coefficients. It is clear that the profile of an $s$-sparse signal is at most $s$. See Figure 1 for images of stylized signals with different profiles.

First, we prove a result on the number of iterations needed to acquire an $s$-sparse signal. At the end of the section, we extend this result to general signals, which yields Theorem 2.2.

Let $\bm{x}$ be an $s$-sparse signal, and define $p={\rm profile}(\bm{x})$. After at most

iterations, CoSaMP produces an approximation $\bm{a}$ that satisfies

For a fixed $s$, the bound on the number of iterations achieves its maximum value at $p=s$. Since $\log_{4/3}5.6<6$, the number of iterations never exceeds $6(s+1)$.

Let us instate some notation that will be valuable in the proof of the theorem. We write $p={\rm profile}(\bm{x})$. For each $k=0,1,2,\dots$, the signal $\bm{a}^{k}$ is the approximation after the $k$th iteration. We abbreviate $S_{k}=\operatorname{supp}(\bm{a}^{k})$, and we define the residual signal $\bm{r}^{k}=\bm{x}-\bm{a}^{k}$. The norm of the residual can be viewed as the approximation error.

For a nonnegative integer $j$, we may define an auxiliary signal

In other words, $\bm{y}^{j}$ is the part of $\bm{x}$ contained in the bands $B_{j}$, $B_{j+1}$, $B_{j+2}$, …. For each $j\in J$, we have the estimate

by definition of the bands. These auxiliary signals play a key role in the analysis.

B.2. Proof of Theorem B.1

The proof of the theorem involves a sequence of lemmas. The first object is to establish an alternative that holds in each iteration. One possibility is that the approximation error is small, which means that the algorithm is effectively finished. Otherwise, the approximation error is dominated by the energy in the unidentified part of the signal, and the subsequent approximation error is a constant factor smaller.

For each iteration $k=0,1,2,\dots$, at least one of the following alternatives holds. Either

Define $T_{k}$ as the merged support that occurs during iteration $k$. The pruning step ensures that the support $S_{k}$ of the approximation at the end of the iteration is a subset of the merged support, so

At the end of the $k$th iteration, the pruned vector $\bm{b}_{s}$ becomes the next approximation $\bm{a}^{k}$, so the estimation and pruning results, Lemmas 4.4 and 4.5, together imply that

Note that the same relation holds trivially for iteration $k=0$ because $\bm{r}^{0}=\bm{x}$ and $S_{0}=\emptyset$.

Suppose that there is an iteration $k\geq 0$ where

We can introduce this bound directly into the inequality (B.6) to obtain the first conclusion (B.3).

Suppose on the contrary that in iteration $k$ we have

Introducing this relation into the inequality (B.6) leads quickly to the conclusion (B.4). We also have the chain of relations

Therefore, the sparse iteration invariant, Theorem 4.1 ensures that (B.5) holds. ∎

The next lemma contains the critical part of the argument. Under the second alternative in the previous lemma, we show that the algorithm completely identifies the support of the signal, and we bound the number of iterations required to do so.

Fix $K=\lfloor p\log_{4/3}(1+4.6\sqrt{s/p})\rfloor$. Assume that (B.4) and (B.5) are in force for each iteration $k=0,1,2,\dots,K$. Then $\operatorname{supp}(\bm{a}^{K})=\operatorname{supp}(\bm{x})$.

First, we check that, once the norm of the residual is smaller than each element of a band, the components in that band persist in the support of each subsequent approximation. Define $J$ to be the set of nonempty bands, and fix a band $j\in J$. Suppose that, for some iteration $k$, the norm of the residual satisfies

Then it must be the case that $B_{j}\subset\operatorname{supp}(\bm{a}^{k})$. If not, then some component $i\in B_{j}$ appears in the residual: $r^{k}_{i}=x_{i}$. This supposition implies that

an evident contradiction. Since (B.5) guarantees that the norm of the residual declines in each iteration, (B.7) ensures that the support of each subsequent approximation contains $B_{j}$.

Next, we bound the number of iterations required to find the next nonempty band $B_{j}$, given that we have already identified the bands $B_{i}$ where $i<j$. Formally, assume that the support $S_{k}$ of the current approximation contains $B_{i}$ for each $i<j$. In particular, the set of missing components $S_{k}^{c}\subset\operatorname{supp}(\bm{y}^{j})$. It follows from relation (B.4) that

We discover that the total number of iterations required to identify all the (nonempty) bands is at most

For each iteration $k\geq\lfloor k_{\star}\rfloor$, it follows that $\operatorname{supp}(\bm{a}^{k})=\operatorname{supp}(\bm{x})$.

It remains to bound $k_{\star}$ in terms of the profile $p$ of the signal. For convenience, we focus on a slightly different quantity. First, observe that $p=\left|{J}\right|$. Using the geometric mean–arithmetic mean inequality, we discover that

To bound the remaining sum, we recall the relation (B.2). Then we invoke Jensen’s inequality and simplify the result.

The final equality holds because the total number of elements in all the bands equals the signal sparsity $s$. Combining these bounds, we reach

Take logarithms, multiply by $p$, and divide through by $\log\beta$. We conclude that the required number of iterations $k_{\star}$ is bounded as

Finally, we check that the algorithm produces a small approximation error within a reasonable number of iterations.

Abbreviate $K=\lfloor p\log(1+4.6\sqrt{s/p})\rfloor$. Suppose that (B.3) never holds during the first $K$ iterations of the algorithm. Under this circumstance, Lemma B.2 implies that both (B.4) and (B.5) hold during each of these $K$ iterations. It follows from Lemma B.3 that the support $S_{K}$ of the $K$th approximation equals the support of $\bm{x}$. Since $S_{K}$ is contained in the merged support $T_{K}$, we see that the vector $\bm{x}|_{T_{K}^{c}}=\bm{0}$. Therefore, the estimation and pruning results, Lemmas 4.4 and 4.5, show that

It follows that there is an iteration $k\leq K$ where (B.3) is in force. Repeated applications of the iteration invariant, Theorem 4.1, allow us to conclude that

B.3. Proof of Theorem 2.2

Finally, we extend the sparse iteration count result to the general case.

Let $\bm{x}$ be an arbitrary signal, and define $p={\rm profile}(\bm{x}_{s})$. After at most

iterations, CoSaMP produces an approximation $\bm{a}$ that satisfies

Let $\bm{x}$ be a general signal, and let $p={\rm profile}(\bm{x}_{s})$. Lemma 6.1 allows us to write the noisy vector of samples $\bm{u}=\bm{\Phi}\bm{x}_{s}+\widetilde{\bm{e}}$. The sparse iteration count result, Theorem B.1, states that after at most

iterations, the algorithm produces an approximation $\bm{a}$ that satisfies

Apply the lower triangle inequality to the left-hand side. Then recall the estimate for the noise in Lemma 6.1, and simplify to reach

where $\nu$ is the unrecoverable energy. ∎

Invoke Theorem B.4. Recall that the estimate for the number of iterations is maximized with $p=s$, which gives an upper bound of $6(s+1)$ iterations, independent of the signal. ∎

References