Subspace Pursuit for Compressive Sensing Signal Reconstruction

I Introduction

Compressive sensing (CS) is a sampling method closely connected to transform coding which has been widely used in modern communication systems involving large scale data samples. A transform code converts input signals, embedded in a high dimensional space, into signals that lie in a space of significantly smaller dimensions. Examples of transform coders include the well known wavelet transforms and the ubiquitous Fourier transform.

Compressive sensing techniques perform transform coding successfully whenever applied to so-called compressible and/or $K$-sparse signals, i.e., signals that can be represented by $K\ll N$ significant coefficients over an $N$-dimensional basis. Encoding of a $K$-sparse, discrete-time signal $\mathbf{x}$ of dimension $N$ is accomplished by computing a measurement vector $\mathbf{y}$ that consists of $m\ll N$ linear projections of the vector $\mathbf{x}$. This can be compactly described via

Here, $\mathbf{\Phi}$ represents an $m\times N$ matrix, usually over the field of real numbers. Within this framework, the projection basis is assumed to be incoherent with the basis in which the signal has a sparse representation .

The work by Donoho and Candès et. al. demonstrated that CS reconstruction is, indeed, a polynomial time problem – albeit under the constraint that more than $2K$ measurements are used. The key observation behind these findings is that it is not necessary to resort to $l_{0}$ optimization to recover x from the under-determined inverse problem; a much easier $l_{1}$ optimization, based on Linear Programming (LP) techniques, yields an equivalent solution, as long as the sampling matrix $\mathbf{\Phi}$ satisfies the so called restricted isometry property (RIP) with a constant parameter.

While LP techniques play an important role in designing computationally tractable CS decoders, their complexity is still highly impractical for many applications. In such cases, the need for faster decoding algorithms - preferably operating in linear time - is of critical importance, even if one has to increase the number of measurements. Several classes of low-complexity reconstruction techniques were recently put forward as alternatives to linear programming (LP) based recovery, which include group testing methods , and algorithms based on belief propagation .

Recently, a family of iterative greedy algorithms received significant attention due to their low complexity and simple geometric interpretation. They include the Orthogonal Matching Pursuit (OMP), the Regularized OMP (ROMP) and the Stagewise OMP (StOMP) algorithms. The basic idea behind these methods is to find the support of the unknown signal sequentially. At each iteration of the algorithms, one or several coordinates of the vector x are selected for testing based on the correlation values between the columns of $\mathbf{\Phi}$ and the regularized measurement vector. If deemed sufficiently reliable, the candidate column indices are subsequently added to the current estimate of the support set of x. The pursuit algorithms iterate this procedure until all the coordinates in the correct support set are included in the estimated support set. The computational complexity of OMP strategies depends on the number of iterations needed for exact reconstruction: standard OMP always runs through $K$ iterations, and therefore its reconstruction complexity is roughly $O\left(KmN\right)$ (see Section IV-C for details). This complexity is significantly smaller than that of LP methods, especially when the signal sparsity level $K$ is small. However, the pursuit algorithms do not have provable reconstruction quality at the level of LP methods. For OMP techniques to operate successfully, one requires that the correlation between all pairs of columns of $\mathbf{\Phi}$ is at most $1/2K$ , which by the Gershgorin Circle Theorem represents a more restrictive constraint than the RIP. The ROMP algorithm can reconstruct all $K$-sparse signals provided that the RIP holds with parameter $\delta_{2K}\leq 0.06/\sqrt{\log K}$, which strengthens the RIP requirements for $l_{1}$-linear programming by a factor of $\sqrt{\log K}$.

The main contribution of this paper is a new algorithm, termed the subspace pursuit (SP) algorithm. It has provable reconstruction capability comparable to that of LP methods, and exhibits the low reconstruction complexity of matching pursuit techniques for very sparse signals. The algorithm can operate both in the noiseless and noisy regime, allowing for exact and approximate signal recovery, respectively. For any sampling matrix $\mathbf{\Phi}$ satisfying the RIP with a constant parameter independent of $K$, the SP algorithm can recover arbitrary $K$-sparse signals exactly from its noiseless measurements. When the measurements are inaccurate and/or the signal is not exactly sparse, the reconstruction distortion is upper bounded by a constant multiple of the measurement and/or signal perturbation energy. For very sparse signals with $K\leq\text{const}\cdot\sqrt{N}$, which, for example, arise in certain communication scenarios, the computational complexity of the SP algorithm is upper bounded by $O\left(mNK\right)$, but can be further reduced to $O\left(mN\log K\right)$ when the nonzero entries of the sparse signal decay slowly.

The basic idea behind the SP algorithm is borrowed from coding theory, more precisely, the $A^{*}$ order-statistic algorithm for additive white Gaussian noise channels. In this decoding framework, one starts by selecting the set of $K$ most reliable information symbols. This highest reliability information set is subsequently hard-decision decoded, and the metric of the parity checks corresponding to the given information set is evaluated. Based on the value of this metric, some of the low-reliability symbols in the most reliable information set are changed in a sequential manner. The algorithm can therefore be seen as operating on an adaptively modified coding tree. If the notion of “most reliable symbol” is replaced by “column of sensing matrix exhibiting highest correlation with the vector ${\bf{y}}$”, the notion of “parity-check metric” by “residual metric”, then the above method can be easily changed for use in CS reconstruction. Consequently, one can perform CS reconstruction by selecting a set of $K$ columns of the sensing matrix with highest correlation that span a candidate subspace for the sensed vector. If the distance of the received vector to this space is deemed large, the algorithm incrementally removes and adds new basis vectors according to their reliability values, until a sufficiently close candidate word is identified. SP employs a search strategy in which a constant number of vectors is expurgated from the candidate list. This feature is mainly introduced for simplicity of analysis: one can easily extend the algorithm to include adaptive expurgation strategies that do not necessarily operate on fixed-sized lists.

In compressive sensing, the major challenge associated with sparse signal reconstruction is to identify in which subspace, generated by not more than $K$ columns of the matrix $\mathbf{\Phi}$, the measured signal $\mathbf{y}$ lies. Once the correct subspace is determined, the non-zero signal coefficients are calculated by applying the pseudoinversion process. The defining character of the SP algorithm is the method used for finding the $K$ columns that span the correct subspace: SP tests subsets of $K$ columns in a group, for the purpose of refining at each stage an initially chosen estimate for the subspace. More specifically, the algorithm maintains a list of $K$ columns of $\mathbf{\Phi}$, performs a simple test in the spanned space, and then refines the list. If $\mathbf{y}$ does not lie in the current estimate for the correct spanning space, one refines the estimate by retaining reliable candidates, discarding the unreliable ones while adding the same number of new candidates. The “reliability property” is captured in terms of the order statistics of the inner products of the received signal with the columns of $\Phi$, and the subspace projection coefficients.

As a consequence, the main difference between ROMP and the SP reconstruction strategy is that the former algorithm generates a list of candidates sequentially, without back-tracing: it starts with an empty list, identifies one or several reliable candidates during each iteration, and adds them to the already existing list. Once a coordinate is deemed to be reliable and is added to the list, it is not removed from it until the algorithm terminates. This search strategy is overly restrictive, since candidates have to be selected with extreme caution. In contrast, the SP algorithm incorporates a simple method for re-evaluating the reliability of all candidates at each iteration of the process.

At the time of writing this manuscript, the authors became aware of the related work by J. Tropp, D. Needell and R. Vershynin , describing a similar reconstruction algorithm. The main difference between the SP algorithm and the CoSAMP algorithm of is in the manner in which new candidates are added to the list. In each iteration, in the SP algorithm, only $K$ new candidates are added, while the CoSAMP algorithm adds $2K$ vectors. This makes the SP algorithm computationally more efficient, but the underlying analysis more complex. In addition, the restricted isometry constant for which the SP algorithm is guaranteed to converge is larger than the one presented in . Finally, this paper also contains an analysis of the number of iterations needed for reconstruction of a sparse signal (see Theorem 6 for details), for which there is no counterpart in the CoSAMP study.

The remainder of the paper is organized as follows. Section II introduces relevant concepts and terminology for describing the proposed CS reconstruction technique. Section III contains the algorithmic description of the SP algorithm, along with a simulation-based study of its performance when compared with OMP, ROMP, and LP methods. Section IV contains the main result of the paper pertaining to the noiseless setting: a formal proof for the guaranteed reconstruction performance and the reconstruction complexity of the SP algorithm. Section V contains the main result of the paper pertaining to the noisy setting. Concluding remarks are given in Section VI, while proofs of most of the theorems are presented in the Appendix of the paper.

II Preliminaries

We are concerned with the problem of low-complexity recovery of the unknown signal $\mathbf{x}$ from the measurement $\mathbf{y}$. A natural formulation of the recovery problem is within an $l_{0}$ norm minimization framework, which seeks a solution to the problem

Unfortunately, the above $l_{0}$ minimization problem is NP-hard, and hence cannot be used for practical applications .

One way to avoid using this computationally intractable formulation is to consider a $l_{1}$-regularized optimization problem,

denotes the $l_{1}$ norm of the vector x.

The main advantage of the $l_{1}$ minimization approach is that it is a convex optimization problem that can be solved efficiently by linear programming (LP) techniques. This method is therefore frequently referred to as $l_{1}$-LP reconstruction , and its reconstruction complexity equals $O\left(m^{2}N^{3/2}\right)$ when interior point methods are employed . See for other methods to further reduce the complexity of $l_{1}$-LP.

The reconstruction accuracy of the $l_{1}$-LP method is described in terms of the restricted isometry property (RIP), formally defined below.

We define $\delta_{K}$, the RIP constant, as the infimum of all parameters $\delta$ for which the RIP holds, i.e.

where $\lambda_{\min}\left(\mathbf{\Phi}_{I}^{*}\mathbf{\Phi}_{I}\right)$ and $\lambda_{\max}\left(\mathbf{\Phi}_{I}^{*}\mathbf{\Phi}_{I}\right)$ denote the minimal and maximal eigenvalues of $\mathbf{\Phi}_{I}^{*}\mathbf{\Phi}_{I}$, respectively.

Most known families of matrices satisfying the RIP property with optimal or near-optimal performance guarantees are random. Examples include:

Random matrices with i.i.d. entries that follow either the Gaussian distribution, Bernoulli distribution with zero mean and variance $1/n$, or any other distribution that satisfies certain tail decay laws. It was shown in that the RIP for a randomly chosen matrix from such ensembles holds with overwhelming probability whenever

where $C$ is a function of the RIP constant.

Random matrices from the Fourier ensemble. Here, one selects $m$ rows from the $N\times N$ discrete Fourier transform matrix uniformly at random. Upon selection, the columns of the matrix are scaled to unit norm. The resulting matrix satisfies the RIP with overwhelming probability, provided that

where $C$ depends only on the RIP constant.

There exists an intimate connection between the LP reconstruction accuracy and the RIP property, first described by Candés and Tao in . If the sampling matrix $\mathbf{\Phi}$ satisfies the RIP with constants $\delta_{K}$, $\delta_{2K}$, and $\delta_{3K}$, such that

then the $l_{1}$-LP algorithm will reconstruct all $K$-sparse signals exactly. This sufficient condition (1) can be improved to

For subsequent derivations, we need two results summarized in the lemmas below. The first part of the claim, as well as a related modification of the second claim also appeared in . For completeness, we include the proof of the lemma in Appendix -A.

(Monotonicity of $\delta_{K}$) For any two integers $K\leq K^{\prime}$,

The lemma implies that $\delta_{K}\leq\delta_{2K}\leq\delta_{3K}$, which consequently simplifies (1) to $\delta_{3K}<1/3$. Both (1) and (2) represent sufficient conditions for exact reconstruction.

In order to describe the main steps of the SP algorithm, we introduce next the notion of the projection of a vector and its residue.

denotes the pseudo-inverse of the matrix $\mathbf{\Phi}_{I}$, and ∗ stands for matrix transposition.

The residue vector of the projection equals

We find the following properties of projections and residues of vectors useful for our subsequent derivations.

The proof of Lemma 2 can be found in Appendix -B.

III The SP Algorithm

We performed extensive computer simulations in order to compare the accuracy of different reconstruction algorithms empirically. In the compressive sensing framework, all sparse signals are expected to be exactly reconstructed as long as the level of the sparsity is below a certain threshold. However, the computational complexity to test this uniform reconstruction ability is $O\left(N^{K}\right)$, which grows exponentially with $K$. Instead, for empirical testing, we adopt the simulation strategy described in which calculates the empirical frequency of exact reconstruction for the Gaussian random matrix ensemble. The steps of the testing strategy are listed below.

For given values of the parameters $m$ and $N$, choose a signal sparsity level $K$ such that $K\leq m/2$;

Randomly generate a $m\times N$ sampling matrix $\mathbf{\Phi}$ from the standard i.i.d. Gaussian ensemble;

Select a support set $T$ of size $\left|T\right|=K$ uniformly at random, and generate the sparse signal vector $\mathbf{x}$ by either one of the following two methods:

Draw the elements of the vector $\mathbf{x}$ restricted to $T$ from the standard Gaussian distribution; we refer to this type of signal as a Gaussian signal. Or,

set all entries of $\mathbf{x}$ supported on $T$ to ones; we refer to this type of signal as a zero-one signal.

Note that zero-one sparse signals are of special interest for the comparative study, since they represent a particularly challenging case for OMP-type of reconstruction strategies.

Compute the measurement $\mathbf{y}=\mathbf{\Phi}\mathbf{x}$, apply a reconstruction algorithm to obtain $\hat{\mathbf{x}}$, the estimate of $\mathbf{x}$, and compare $\hat{\mathbf{x}}$ to $\mathbf{x}$;

Repeat the process $500$ times for each $K$, and then simulate the same algorithm for different values of $m$ and $N$.

The improved reconstruction capability of the SP method, compared with that of the OMP and ROMP algorithms, is illustrated by two examples shown in Fig. 2. Here, the signals are drawn both according to the Gaussian and zero-one model, and the benchmark performance of the LP reconstruction technique is plotted as well.

Figure 2 depicts the empirical frequency of exact reconstruction. The numerical values on the $x$-axis denote the sparsity level $K$, while the numerical values on the $y$-axis represent the fraction of exactly recovered test signals. Of particular interest is the sparsity level at which the recovery rate drops below 100% - i.e. the critical sparsity - which, when exceeded, leads to errors in the reconstruction algorithm applied to some of the signals from the given class.

The simulation results reveal that the critical sparsity of the SP algorithm by far exceeds that of the OMP and ROMP techniques, for both Gaussian and zero-one inputs. The reconstruction capability of the SP algorithm is comparable to that of the LP based approach: the SP algorithm has a slightly higher critical sparsity for Gaussian signals, but also a slightly lower critical sparsity for zero-one signals. However, the SP algorithms significantly outperforms the LP method when it comes to reconstruction complexity. As we analytically demonstrate in the exposition to follow, the reconstruction complexity of the SP algorithm for both Gaussian and zero-one sparse signals is $O\left(mN\log K\right)$, whenever $K\leq O\left(\sqrt{N}\right)$, while the complexity of LP algorithms based on interior point methods is $O\left(m^{2}N^{3/2}\right)$ in the same asymptotic regime.

IV Recovery of Sparse Signals

A sufficient condition for exact reconstruction of arbitrary sparse signals is stated in the following theorem.

then the SP algorithm is guaranteed to exactly recover $\mathbf{x}$ from $\mathbf{y}$ via a finite number of iterations.

The requirement on RIP constant can be relaxed to

In the original version of this manuscript, we proved the weaker result $\delta_{3K}\leq 0.06$. At the time of revision of the paper, we were given access to the manuscript by Needel and Tropp. Using some of the proof techniques in their work, we managed to improve the results in Theorem 3 and therefore the RIP constant of the original submission. The interested reader is referred to http://arxiv.org/abs/0803.0811v2 for the first version of the theorem. This paper contains only the proof of the stronger result.

This sufficient condition is proved by applying Theorems 2 and 6. The computational complexity is related to the number of iterations required for exact reconstruction, and is discussed at the end of Section IV-C. Before providing a detailed analysis of the results, let us sketch the main ideas behind the proof.

which implies that at each iteration, the SP algorithm identifies a $K$-dimensional space that reduces the reconstruction error of the vector $\mathbf{x}$. See Fig. 3 for an illustration. This observation is formally stated as follows.

Assume that the conditions of Theorem 1 hold. For each iteration of the SP algorithm, one has

The proof of the theorem is postponed to Appendix -D.

The proof of the result is deferred to Appendix -E.

Based on Theorems 3 and 4, one arrives at the result claimed in Equation (6).

Furthermore, according to Lemmas 1 and 2, one has

where the second equality holds by the definition of the residue, while (4) and (6) refer to the labels of the inequalities used in the bounds. In addition,

Upon combining (9) and (10), one obtains the following upper bound

Finally, elementary calculations show that when $\delta_{3K}\leq 0.165$,

Both in the initialization step and during each iteration of the SP algorithm, we select $K$ indices that maximize the correlations between the column vectors and the residual measurement. Henceforth, this step is referred to as correlation maximization (CM). Consider the ideal case where all columns of $\mathbf{\Phi}$ are orthogonalOf course, in this case no compression is possible.. In this scenario, the signal coefficients can be easily recovered by calculating the correlations $\left\langle\mathbf{v}_{i},\mathbf{y}\right\rangle$ - i.e., all indices with non-zero magnitude are in the correct support of the sensed vector. Now assume that the sampling matrix $\mathbf{\Phi}$ satisfies the RIP. Recall that the RIP (see Lemma 1) implies that the columns are locally near-orthogonal. Consequently, for any $j$ not in the correct support, the magnitude of the correlation $\left\langle\mathbf{v}_{j},\mathbf{y}\right\rangle$ is expected to be small, and more precisely, upper bounded by $\delta_{K+1}\left\|\mathbf{x}\right\|_{2}$. This seems to provide a very simple intuition why correlation maximization allows for exact reconstruction. However, this intuition is not easy to analytically justify due to the following fact. Although it is clear that for all indices $j\notin T$, the values of $\left|\left\langle\mathbf{v}_{j},\mathbf{y}\right\rangle\right|$ are upper bounded by $\delta_{K+1}\left\|\mathbf{x}\right\|$, it may also happen that for all $i\in T$, the values of $\left|\left\langle\mathbf{v}_{i},\mathbf{y}\right\rangle\right|$ are small as well. Dealing with maximum correlations in this scenario cannot be immediately proved to be a good reconstruction strategy. The following example illustrates this point.

Without loss of generality, let $T=\left\{1,\cdots,K\right\}$. Let the vectors $\mathbf{v}_{i}$ ($i\in T$) be orthonormal, and let the remaining columns $\mathbf{v}_{j}$, $j\notin T$, of $\mathbf{\Phi}$ be constructed randomly, using i.i.d. Gaussian samples. Consider the following normalized zero-one sparse signal

It is straightforward to envision the existence of an index $j\notin T$, such that

The latter inequality is critical, because achieving very small values for the RIP constant is a challenging task.

This example represents a particularly challenging case for the OMP algorithm. Therefore, one of the major constraints imposed on the OMP algorithm is the requirement that

To meet this requirement, $\delta_{K+1}$ has to be less than $1/\sqrt{K}$, which decays fast as $K$ increases.

In contrast, the SP algorithm allows for the existence of some index $j\notin T$ with

As long as the RIP constant $\delta_{3K}$ is upper bounded by the constant given in (5), the indices in the correct support of $\mathbf{x}$, that account for the most significant part of the energy of the signal, are captured by the CM procedure. Detailed descriptions of how this can be achieved are provided in the proofs of the previously stated Theorems 3 and 5.

Let us first focus on the initialization step. By the definition of the set $T^{0}$ in the initialization stage of the algorithm, the set of the $K$ selected columns ensures that

Now, if we assume that the estimate $T^{0}$ is disjoint from the correct support, i.e., that $T^{0}\bigcap T=\phi$, then by the near orthogonality property of Lemma 1, one has

The last inequality clearly contradicts (11) whenever $\delta_{K}\leq\delta_{2K}<1/2$. Consequently, if $\delta_{2K}<1/2$, then

and at least one correct element of the support of x is in $T^{0}$. This phenomenon is quantitatively described in Theorem 5.

The proof of the theorem is postponed to Appendix -C.

To study the effect of correlation maximization during each iteration, one has to observe that correlation calculations are performed with respect to the vector

instead of being performed with respect to the vector $\mathbf{y}$. As a consequence, to show that the CM process captures a significant part of residual signal energy requires an analysis including a number of technical details. These can be found in the Proof of Theorem 3.

IV-B Identifying Indices Outside of the Correct Support Set

and that it contains at least $K$ zeros. Consequently, the $K$ indices with smallest magnitude - equal to zero - are clearly not in the correct support set.

IV-C Convergence of the SP Algorithm

In this subsection, we upper bound the number of iterations needed to reconstruct an arbitrary $K$-sparse signal using the SP algorithm.

Given an arbitrary $K$-sparse signal $\mathbf{x}$, we first arrange its elements in decreasing order of magnitude. Without loss of generality, assume that

and that $x_{j}=0,\;\forall\,j>K$. Define

The number of iterations of the SP algorithm is upper bounded by

This result is a combination of Theorems 7 and (12),The upper bound in Theorem 7 is also obtained in while the one in Theorem 8 is not. described below.

The proof of Theorem 7 is intuitively clear and presented below, while the proof of Theorem 8 is more technical and postponed to Appendix -F.

A drawback of Theorem 7 is that it sometimes overestimates the number of iterations, especially when $\rho_{\min}\ll 1$. The example to follow illustrates this point.

Let $K=2$, $x_{1}=2^{10},$ $x_{2}=1$, $x_{3}=\cdots=x_{N}=0$. Suppose that the sampling matrix $\mathbf{\Phi}$ satisfies the RIP with $c_{K}=\frac{1}{2}.$ Noting that $\rho_{\min}\lesssim 2^{-10}$, Theorem 6 implies that

Indeed, if we take a close look at the steps of the SP algorithm, we can verify that

After the initialization step, by Theorem 5, it can be shown that

As a result, the estimate $T^{0}$ must contain the index one and $\left\|\mathbf{x}_{T-T^{0}}\right\|_{2}\leq 1$. After the first iteration, since

It is clear that the number of iterations required for exact reconstruction depends on the values of the entries of the sparse signal. We therefore focus our attention on the following three particular classes of sparse signals.

Zero-one sparse signals. As explained before, zero-one signals represent the most challenging reconstruction category for OMP algorithms. However, this class of signals has the best upper bound on the convergence rate of the SP algorithm. Elementary calculations reveal that $\rho_{\min}=1/\sqrt{K}$ and that

Sparse signals with power-law decaying entries (also known as compressible sparse signals). Signals in this category are defined via the following constraint

for some constants $c_{x}>0$ and $p>1$. Compressible sparse signals have been widely considered in the CS literature, since most practical and naturally occurring signals belong to this class . It follows from Theorem 7 that in this case

where $o\left(1\right)\rightarrow 0$ when $K\rightarrow\infty$.

Sparse signals with exponentially decaying entries. Signals in this class satisfy

for some constants $c_{x}>0$ and $p>0$. Theorem 6 implies that

where again $o\left(1\right)\rightarrow 0$ as $K\rightarrow\infty$.

With the bound on the number of iterations required for exact reconstruction at hand, the computational complexity of the SP algorithm can be easily estimated: it equals the complexity of one iteration multiplied by the number of iterations. In each iteration, CM requires $mN$ computations in general. For some measurement matrices with special structures, for example, sparse matrices, the computational cost can be reduced significantly. The cost of computing the projections is of the order of $O\left(K^{2}m\right)$, if one uses the Modified Gram-Schmidt (MGS) algorithm [20, pg. 61]. This cost can be reduced further by “reusing” the computational results of past iterations within future iterations. This is possible because most practical sparse signals are compressible, and the signal support set estimates in different iterations usually intersect in a large number of indices. Though there are many ways to reduce the complexity of both the CM and projection computation steps, we only focus on the most general framework of the SP algorithm, and assume that the complexity of each iteration equals $O\left(mN+mK^{2}\right)$. As a result, the total complexity of the SP algorithm is given by $O\left(m\left(N+K^{2}\right)\log K\right)$ for compressible sparse signals, and it is upper bounded by $O\left(m\left(N+K^{2}\right)K\right)$ for arbitrary sparse signals. When the signal is very sparse, in particular, when $K^{2}\leq O\left(N\right)$, the total complexity of SP reconstruction is upper bounded by $O\left(mNK\right)$ for arbitrary sparse signals and by $O\left(mN\log K\right)$ for compressible sparse signals (we once again point out that most practical sparse signals belong to this signal category ).

The complexity of the SP algorithm is comparable to OMP-type algorithms for very sparse signals where $K^{2}\leq O\left(N\right)$. For the standard OMP algorithm, exact reconstruction always requires $K$ iterations. In each iteration, the CM operation costs $O\left(mN\right)$ computations and the complexity of the projection is marginal compared with the CM. The corresponding total complexity is therefore always $O\left(mNK\right)$. For the ROMP and StOMP algorithms, the challenging signals in terms of convergence rate are also the sparse signals with exponentially decaying entries. When the $p$ in (13) is sufficiently large, it can be shown that both ROMP and StOMP also need $O\left(K\right)$ iterations for reconstruction. Note that CM operation is required in both algorithms. The total computational complexity is then $O\left(mNK\right)$.

The case that requires special attention during analysis is $K^{2}>O\left(N\right)$. Again, if compressible sparse signals are considered, the complexity of projections can be significantly reduced if one reuses the results from previous iterations at the current iteration. If exponentially decaying sparse signals are considered, one may want to only recover the energetically most significant part of the signal and treat the residual of the signal as noise — reduce the effective signal sparsity to $K^{\prime}\ll K$. In both cases, the complexity depends on the specific implementation of the CM and projection operations and is beyond the scope of analysis of this paper.

One advantage of the SP algorithm is that the number of iterations required for recovery is significantly smaller than that of the standard OMP algorithm for compressible sparse signals. To the best of the authors’ knowledge, there are no known results on the number of iterations of the ROMP and StOMP algorithms needed for recovery of compressible sparse signals.

V Recovery of Approximately Sparse Signals from Inaccurate Measurements

We first consider a sampling scenario in which the signal $\mathbf{x}$ is $K$-sparse, but the measurement vector $\mathbf{y}$ is subjected to an additive noise component, $\mathbf{e}$. The following theorem gives a sufficient condition for convergence of the SP algorithm in terms of the RIP constant $\delta_{3K}$, as well as an upper bounds on the recovery distortion that depends on the energy ($l_{2}$-norm) of the error vector $\mathbf{e}$.

Then the reconstruction distortion of the SP algorithm satisfies

The proof of the theorem is given in Section V-A.

We also study the case where the signal $\mathbf{x}$ is only approximately $K$-sparse, and the measurement $\mathbf{y}$ is contaminated by a noise vector $\mathbf{e}$. To simplify the notation, we henceforth use $\mathbf{x}_{K}$ to denote the vector obtained from $\mathbf{x}$ by maintaining the $K$ entries with largest magnitude and setting all other entries in the vector to zero. In this setting, a signal $\mathbf{x}$ is said to be approximately $K$-sparse if $\mathbf{x}-\mathbf{x}_{K}\neq\mathbf{0}$. Based on Theorem 9, we can upper bound the recovery distortion in terms of the $l_{1}$ and $l_{2}$ norms of $\mathbf{x}-\mathbf{x}_{K}$ and $\mathbf{e}$, respectively, as follows.

The proof of this corollary is given in Section V-B. As opposed to the standard case where the input sparsity level of the SP algorithm equals the signal sparsity level $K$, one needs to set the input sparsity level of the SP algorithm to $2K$ in order to obtain the claim stated in the above corollary.

Theorem 9 and Corollary 1 provide analytical upper bounds on the reconstruction distortion of the noisy version of the SP algorithm. In addition to these theoretical bounds, we performed numerical simulations to empirically estimate the reconstruction distortion. In the simulations, we first select the dimension $N$ of the signal $\mathbf{x}$, and the number of measurements $m$. We then choose a sparsity level $K$ such that $K\leq m/2$. Once the parameters are chosen, an $m\times N$ sampling matrix with standard i.i.d. Gaussian entries is generated. For a given $K$, the support set $T$ of size $\left|T\right|=K$ is selected uniformly at random. A zero-one sparse signal is constructed as in the previous section. Finally, either signal or a measurement perturbations are added as follows:

Signal perturbations: the signal entries in $T$ are kept unchanged but the signal entries outside of $T$ are perturbed by i.i.d. Gaussian $\mathcal{N}\left(0,\sigma_{s}^{2}\right)$ samples.

Measurement perturbations: the perturbation vector $\mathbf{e}$ is generated using a Gaussian distribution with zero mean and covariance matrix $\sigma_{e}^{2}\mathbf{I}_{m}$, where $\mathbf{I}_{m}$ denotes the $m\times m$ identity matrix.

We ran the SP reconstruction process on $\mathbf{y}$, $500$ times for each $K$, $\sigma_{s}^{2}$ and $\sigma_{e}^{2}$. The reconstruction distortion $\left\|\mathbf{x}-\hat{\mathbf{x}}\right\|_{2}$ is obtained via averaging over all these instances, and the results are plotted in Fig. 6. Consistent with the findings of Theorem 9 and Corollary 1, we observe that the recovery distortion increases linearly with the $l_{2}$-norm of the measurement error. Even more encouraging is the fact that the empirical reconstruction distortion is typically much smaller than the corresponding upper bounds. This is likely due to the fact that, in order to simplify the expressions involved, many constants and parameters used in the proof were upper bounded.

The first step towards proving Theorem 9 is to upper bound the reconstruction error for a given estimated support set $\hat{T}$, as succinctly described in the lemma to follow.

The proof of the lemma is given in Appendix -G.

The upper bounds in Inequalities (15) and (16) are proved in Appendix -H and -I, respectively. The inequality (17) is obtained by substituting (15) into (16) as shown below:

To complete the proof, we make use of Lemma 2 stated in Section II. According to this lemma, we have

Based on Theorem 10, we conclude that when the SP algorithm terminates, the inequality (18) is violated and we must have

Under this assumption, it follows from Lemma 3 that

V-B Recovery Distortion under Signal and Measurement Perturbations

The proof of Corollary 1 is based on the following two lemmas, which are proved in and , respectively.

To prove the corollary, consider the measurement vector

VI Conclusion

We introduced a new algorithm, termed subspace pursuit, for low-complexity recovery of sparse signals sampled by matrices satisfying the RIP with a constant parameter $\delta_{3K}$. Also presented were simulation results demonstrating that the recovery performance of the algorithm matches, and sometimes even exceeds, that of the LP programming technique; and, simulations showing that the number of iterations executed by the algorithm for zero-one sparse signals and compressible signals is of the order $O(\log\,K)$.

VII Acknowledgment

The authors are grateful to Prof. Helmut Böölcskei for handling the manuscript, and the reviewers for their thorough and insightful comments and suggestions.

We provide next detailed proofs for the lemmas and theorems stated in the paper.

obviously holds if either one of the norms $\left\|\mathbf{a}\right\|_{2}$ and $\left\|\mathbf{b}\right\|_{2}$ is zero. Assume therefore that neither one of them is zero, and define

-B Proof of Lemma 2

The first claim is proved by observing that

To prove the second part of the lemma, let

On the other hand, the left hand side of the above inequality reads as

and $\left\|\mathbf{y}_{p}\right\|_{2}^{2}\geq 0$, we show that

-C Proof of Theorem 5

The first step consists in proving Inequality (11), which reads as

By assumption, $\left|T\right|\leq K$, so that

The second step is to partition the estimate of the support set $T^{0}$ into two subsets: the set $T^{0}\bigcap T$, containing the indices in the correct support set, and $T^{0}-T$, the set of incorrectly selected indices. Then

where the last inequality follows from the near-orthogonality property of Lemma 1.

Combining the two inequalities (24) and (25), one can show that

By invoking Inequality (11) it follows that

-D Proof of Theorem 3

In this section we show that the CM process allows for capturing a significant part of the residual signal power, that is,

For notational convenience, we first define

which is the set of indices “captured” by the CM process. By the definition of $T_{\Delta}$, we have

An upper bound on the left hand side of (31) is given by

A lower bound on the right hand side of (31) can be derived as

Substitute (33) and (32) into (31). We get

-E Proof of Theorem 4

be the projection coefficient vector, and let

This result implies that the energy concentrated in the erroneously removed signal components is small.

On the other hand, by Step 4) of the subspace algorithm, $\Delta T$ is chosen to contain the $K$ smallest projection coefficients (in magnitude). It therefore holds that

Next, we decompose the vector $\left(\mathbf{x}_{p}\right)_{\Delta T}$ into a signal part and a smear part. Then

Combining (43) and (44) and noting that $\mathbf{x}_{\Delta T}=\mathbf{x}_{T\bigcap\Delta T}$ ($\mathbf{x}$ is supported on $T$, i.e., $\mathbf{x}_{T^{c}}=\mathbf{0}$), we have

This completes the proof of the claimed result.

-F Proof of Theorem 8

The following iterative algorithm is employed to create a partition of the support set $T$ that will establish the correctness of the claimed result.

Suppose that after the iterative partition, we have

where $J\leq K$ is the number of the subsets in the partition. Let $s_{j}=\left|T_{j}\right|$, $j=1,\cdots,J$. It is clear that

Then Theorem 8 is proved by invoking the following lemma.

For a given index $j$, let $\left|T_{j}\right|=s$, and let

where $\left\lfloor\cdot\right\rfloor$ denotes the floor function. Then, for any $1\leq j_{0}\leq J$, after

iterations, the SP algorithm has the property that

iterations, the SP algorithm guarantees that $T\subset T^{n}$.

Both parts of this lemma are proved by mathematical induction as follows.

or equivalently, the desired inequality (48) holds for $k=1$. To use mathematical induction, suppose that for an index $1<k\leq s-1$,

This proves Equation (48) of the lemma. Inequality (49) then follows from the observation that

From (50), it is clear that for $1\leq j\leq J$,

According to Theorem 2, after $n_{1}$ iterations,

On the other hand, for any $i\in T_{1}$, it follows from the first part of this lemma that

Let $T_{0}=\bigcup_{j=1}^{j_{0}-1}T_{j}$. Then

Denote the smallest coordinate in $T_{j_{0}}$ by $i$, and the largest coordinate in $T_{j_{0}}$ by $k$. Then

This completes the proof of the last claim (52).

-G Proof of Lemma 3

The claim in the lemma is established through the following chain of inequalities:

where $\left(a\right)$ is a consequence of the fact that

By relaxing the upper bound in terms of replacing $\delta_{2K}$ by $\delta_{3K}$, we obtain

-H Proof of Inequality (15)

The proof is similar to the proof given in Appendix -D. We start with observing that

The left hand side of (57) is upper bounded by

Comparing the above inequality (59) with its analogue for the noiseless case, (30), one can see that the only difference is the $2\sqrt{1+\delta_{K}}\left\|\mathbf{e}\right\|_{2}$ term on the left hand side of (59). Following the same steps as used in the derivations leading from (30) to (33), one can show that

-I Proof of Inequality (16)

This proof is similar to that of Theorem 4. When there are measurement perturbations, one has

Then the smear energy is upper bounded by

Invoking the same technique as used for deriving (41), we have

It is straightforward to verify that (45) still holds, which now reads as

References