Phase Transitions for Greedy Sparse Approximation Algorithms

Introduction

In compressed sensing , one works under the sparse approximation assumption, namely, that signals/vectors of interest can be well approximated by few components of a known basis. This assumption is often satisfied due to constraints imposed by the system which generates the signal. In this setting, it has been proven (originally in and by many others since) that the number of linear observations of the signal, required to guarantee recovery, need only be proportional to the sparsity of the signal’s approximation. This is in stark contrast to the standard Shannon-Nyquist Sampling paradigm where worst-case sampling requirements are imposed.

where $\|\cdot\|_{2}$ denotes the Euclidean norm.

However, solving (1) via a naive exhaustive search is combinatorial in nature and NP-hard . A major aspect of compressed sensing theory is the study of alternative methods to solving (1). Since the system $y=Ax$ is underdetermined, any successful recovery of $x$ will require some form of nonlinear reconstruction. Under certain conditions, various algorithms have been shown to successfully reduce (1) to a tractable problem, one with a computational cost which is a low degree polynomial of the problem dimensions, rather than the exponential cost associated with a direct combinatorial search for the solution of (1). While there are numerous reconstruction algorithms, they each generally fall into one of three categories: greedy methods, regularizations, or combinatorial group testing. For an indepth discussion of compressed sensing recovery algorithms, see and references therein.

More precisely, the main contributions of this article is the derivation of theorems and corollaries of the following form for each of the CoSaMP, SP, and IHT algorithms.

Given a matrix $A$ with entries drawn i.i.d. from ${\cal N}(0,n^{-1})$, for any $x\in\chi^{N}(k)$, let $y=Ax+e$ for some (unknown) noise vector $e$. For any $\epsilon\in(0,1)$, as $(k,n,N)\rightarrow\infty$ with $n/N\rightarrow\delta\in(0,1)$ and $k/n\rightarrow\rho\in(0,1)$, there exists $\mu^{alg}(\delta,\rho)$ and $\rho_{S}^{alg}(\delta)$, the unique solution to $\mu^{alg}(\delta,\rho)=1$. If $\rho<(1-\epsilon)\rho_{S}^{alg}(\delta)$, there is an exponentially high probability on the draw of $A$ that the output of the algorithm at the $l^{th}$ iteration, $\hat{x}$, approximates $x$ within the bound

for some $\kappa^{alg}(\delta,\rho)$ and $\xi^{alg}(\delta,\rho)$.

Given a matrix $A$ with entries drawn i.i.d. from ${\cal N}(0,n^{-1})$, for any $x\in\chi^{N}(k)$, let $y=Ax$. For any $\epsilon\in(0,1)$, with $n/N\rightarrow\delta\in(0,1)$ and $k/n\rightarrow\rho<(1-\epsilon)\rho_{S}^{alg}(\delta)$ as $(k,n,N)\rightarrow\infty$, there is an exponentially high probability on the draw of $A$ that the algorithm exactly recovers $x$ from $y$ and $A$ in a finite number of iterations not to exceed

with $T:=\{i:x_{i}\neq 0\}$ and $\lceil m\rceil$, the smallest integer greater than or equal to $m$.

The factors $\mu^{alg}(\delta,\rho)$ and $\frac{\xi^{alg}}{1-\mu^{alg}}(\delta,\rho)$ for CoSaMP, SP, and IHT are displayed in Figure 2, while formulae for their calculation are deferred to Section 3.

Corollary 2 implies that $\rho^{alg}_{S}(\delta)$ delineates the region in which the algorithm can be guaranteed to converge provided there exists an $x\in\chi^{N}(k)$ such that $y=Ax$. However, if no such $x$ exists, as $\rho$ approaches $\rho^{alg}_{S}(\delta)$ the guarantees on the number of iteraties required and stability factors become unbounded. Further bounds on the convergence factor $\mu^{alg}(\delta,\rho)$ and the stability factor $\frac{\xi^{alg}}{1-\mu^{alg}}(\delta,\rho)$ result in yet lower curves $\rho_{S}^{alg}(\delta;bound)$ for a specified $bound$; recall that $\rho^{alg}_{S}(\delta)$ corresponds to the bound $\mu^{alg}(\delta,\rho)=1$.

In the next section, we recall the three algorithms and introduce necessary notation. Then we present the asymmetric restricted isometry property and formulate weaker restricted isometry conditions on a matrix $A$ that ensure the respective algorithm will successfully recover all $k$-sparse signals. In addition to exact recovery, we study the known bounds on the behavior of the algorithms in the presence of noisy measurements. In order to make quantitative comparisons of these results, we must select a matrix ensemble for analysis. In Section 3, we present the lower bounds on the phase transition for each algorithm when the measurement matrix is a Gaussian random matrix. Phase transitions are developed in the case of exact sparse signals while bounds on the multiplicative stability constants are also compared through associated level curves. Section 4 is a discussion of our interpretation of these results and how to use this phase transition framework for comparison of other algorithms.

Greedy Algorithms and the Asymmetric Restricted Isometry Property

The CoSaMP recovery algorithm is a support recovery algorithm which applies hard thresholding by selecting the $k$ largest entries of a vector obtained by applying a pseudoinverse to the measurement $y$. In CoSaMP, the columns of $A$ selected for the pseudoinverse are obtained by applying hard thresholding of magnitude $2k$ to $A^{*}$ applied to the residual from the previous iteration and adding these indices to the approximate support set from the previous iteration. This larger pseudoinverse matrix of size $2k\times n$ imposes the most stringent aRIP condition of the three algorithms. However, CoSaMP uses one fewer pseudoinverse per iteration than SP as the residual vector is computed with a direct matrix-vector multiply of size $n\times k$ rather than with an additional pseudoinverse. Furthermore, when computing the output vector $\hat{x}$, CoSaMP does not need to apply another pseudoinverse as does SP. See Algorithm 1.

2 Subspace Pursuit

The Subspace Pursuit algorithm is also a support recovery algorithm which applies hard thresholding of magnitude $k$ to a vector obtained by applying a pseudoinverse to the measurements $y$. The submatrix chosen for the pseudoinverse has its columns selected by applying $A^{*}$ to the residual vector from the previous iteration, hard thresholding of magnitude $k$, and adding the indices of the terms to the previous approximate support set. Compared to the other two algorithms, a computational disadvantage of SP is that the aforementioned residual vector is also computed via a pseudoinverse, this time selecting the columns from $A$ by again applying a hard threshold of magnitude $k$. The computation of the approximation to the target signal also requires the application of a pseudoinverse for a matrix of size $n\times k$. See Algorithm 2.

3 Iterative Hard Thresholding

Iterative Hard Thresholding (IHT) is also a support recovery algorithm. However, IHT applies hard thresholding to an approximation of the target signal, rather than to the residuals. This completely eliminates the use of a pseudoinverse, reducing the computational cost per iteration. In particular, hard thresholding of magnitude $k$ is applied to an updated approximation of the target signal, $x$, obtained by matrix-vector multiplies of size $n\times N$ that represent a move by a fixed stepsize $\omega$ along the steepest descent direction from the current iterate for the residual $\|Ax-y\|_{2}^{2}$. See Algorthm 3.

(Stopping criteria for greedy methods) In the case of corrupted measurements, where $y=Ax+e$ for some noise vector $e$, the stopping criteria listed in Algorithms 1-3 may never be achieved. Therefore, a suitable alternative stopping criteria must be employed. For our analysis on bounding the error of approximation in the noisy case, we bound the approximation error if the algorithm terminates after $l$ iterations. For example, we could change the algorithm to require a maximum number of iterations $l$ as an input and then terminate the algorithm if our stopping criteria is not met in fewer iterations. In practice, the user would be better served to stop the algorithm when the residual is no longer improving. For a more thorough discussion of suitable stopping criteria for each algorithm in the noisy case, see the original announcement of the algorithms .

4 The Asymmetric Restricted Isometry Property

In this section we relax the sufficient conditions originally placed on Algorithms 1-3 by employing a more general notion of a restricted isometry. As discussed in , the singular values of the $n\times k$ submatrices of an arbitrary measurement matrix $A$ do not, in general, deviate from unity symmetrically. The standard notion of the restricted isometry property (RIP) has an inherent symmetry which is unneccessarily restrictive. Hence, seeking the best possible conditions for the measurement matrix under which Algorithms 1-3 will provably recovery every $k$ sparse vector, we reformulate the sufficient conditions in terms of the asymmetric restricted isometry property (aRIP) .

For an $n\times N$ matrix $A$, the asymmetric RIP constants $L(k,n,N)$ and $U(k,n,N)$ are defined as:

The more common, symmetric definition of the RIP constants is recovered by defining $R(k,n,N)=\max\{L(k,n,N),U(k,n,N)\}$. In this case, a matrix $A$ of size $n\times N$ has the RIP constant $R(k,n,N)$ if

Observe that $\chi^{N}(k)\subset\chi^{N}(k+1)$ for any $k$ and therefore the constants $L(k,n,N)$, $U(k,n,N)$, and $R(k,n,N)$ are nondecreasing in $k$ .

For all expressions involving $L(\cdot,n,N)$ it is understood, without explicit statement, that the first argument is limited to the range where $L(\cdot,n,N)<1$. Beyond this range of sparsity, there exist vectors which are mapped to zero, and are unrecoverable.

Using the aRIP, we analyze the three algorithms in the case of a general measurement matrix $A$ of size $n\times N$. For each algorithm, the application of Definition 1 results in a relaxation of the conditions imposed on $A$ to provably guarantee recovery of all $x\in\chi^{N}(k)$. We first present a stability result for each algorithm in terms of bounding the approximation error of the output after $l$ iterations. The bounds show a multiplicative stability constant in terms of aRIP contants that amplifies the total energy of the noise. As a corollary, we obtain a sufficient condition on $A$ in terms of the aRIP for exact recovery of all $k$-sparse vectors. The proofs of these results are found in the Appendix. These theorems and corollaries take the same form, differing for each algorithm only by the formulae for various factors. We state the general form of the theorems and corollaries, analogous to Theorem 1 and Corollary 2, and then state the formulae for each of the algorithms CoSaMP, SP, and IHT.

Given a matrix $A$ of size $n\times N$ with aRIP constants $L(\cdot,n,N)$ and $U(\cdot,n,N)$, for any $x\in\chi^{N}(k)$, let $y=Ax+e$, for some (unknown) noise vector $e$. Then there exists $\mu^{alg}(k,n,N)$ such that if $\mu^{alg}(k,n,N)<1$, the output $\hat{x}$ of algorithm “alg” at the $l^{th}$ iteration approximates $x$ within the bound

for some $\kappa^{alg}(k,n,N)$ and $\xi^{alg}(k,n,N)$.

Given a matrix $A$ of size $n\times N$ with aRIP constants $L(\cdot,n,N)$ and $U(\cdot,n,N)$, for any $x\in\chi^{N}(k)$, let $y=Ax$. Then there exists $\mu^{alg}(k,n,N)$ such that if $\mu^{alg}(k,n,N)<1$, the algorithm “alg” exactly recovers $x$ from $y$ and $A$ in a finite number of iterations not to exceed

where $\nu_{\text{min}}(x)$ defined as in (5).

We begin with Algorithm 1, the Compressive Sampling Matching Pursuit recovery algorithm of Needell and Tropp . We relax the sufficient recovery condition in via the aRIP.

Theorem 3 and Corollary 4 are satisfied by CoSaMP, Algorithm 1, with $\kappa^{csp}(k,n,N):=1$ and $\mu^{csp}(k,n,N)$ and $\xi^{csp}(k,n,N)$ defined as

Next, we apply the aRIP to Algorithm 2, Dai and Milenkovic’s Subspace Pursuit . Again, the aRIP provides a sufficient condition that admits a wider range of measurement matrices than admitted by the symmetric RIP condition derived in .

Theorem 3 and Corollary 4 are satisfied by Subspace Pursuit, Algorithm 2, with $\kappa^{sp}(k,n,N)$, $\mu^{sp}(k,n,N)$, and $\xi^{sp}(k,n,N)$ defined as

Finally, we apply the aRIP analysis to Algorithm 3, Iterative Hard Thresholding for Compressed Sensing introduced by Blumensath and Davies . Theorem 7 employs the aRIP to provide a weaker sufficient condition than derived in .

Theorem 3 and Corollary 4 are satisfied by Iterative Hard Thresholding, Algorithm 3, with $\kappa^{iht}(k,n,N):=1$ and $\mu^{iht}(k,n,N)$ and $\xi^{iht}(k,n,N)$ defined as

Each of Theorems 5, 6 and 7 are derived following the same recipe as in , and , respectively, using the aRIP rather than the RIP and taking care to maintain the least restrictive bounds at each step (for details, see the Appendix). For Gaussian matrices, the aRIP improves the lower bound on the phase transitions by nearly a multiple of 2 when compared to similar statements using the classical RIP. For IHT, the aRIP is simply a scaling of the matrix so that its RIP bounds are minimal. This is possible for IHT as the factors in $\mu^{iht}(k,n,N)$ involve $L(\alpha k,n,N)$ and $U(\alpha k,n,N)$ for only one value of $\alpha$, here $\alpha=3$. No such scaling interpretation is possible for CoSaMP and SP.

At this point, we digress to mention that the first greedy algorithm shown to have guaranteed exact recovery capability is Needell and Vershynin’s ROMP (Regularized Orthogonal Matching Pursuit) . We omit the algorithm and a rigorous discussion of the result, but state an aRIP condition that will guarantee sparse recovery. ROMP chooses additions to the approximate support sets at each iteration with a regularization step requiring comparability between the added terms. This comparability requires a proof of partitioning a vector of length $N$ into subsets with comparable coordinates, namely the magnitudes of the elements of the subset differ by no more than a factor of 2. The proof that such a partition exists, with each partition having a nonzero energy, forces a pessimistic bound that decays with the problem size.

Let $A$ be a matrix of size $n\times N$ with aRIP constants $L(2k,n,N)$ and $U(2k,n,N)$. Define

If $\mu^{r}(k,n,N)<\left(1+\sqrt{\frac{5n}{n-1}(\log n+2)}\right)^{-1}$, then ROMP is guaranteed to exactly recover any $x\in\chi^{N}(k)$ from the measurements $y=Ax$ in a finite number of iterations.

Unfortunately, this dependence of the bound on the size of the problem instance forces the result to be inadequate for large problem instances. In fact, this result is inferior to the results for the three algorithms stated above which are all independent of problem size and therefore applicable to the most interesting cases of compressed sensing, when $(k,n,N)\rightarrow\infty$ and $\delta=n/N\rightarrow 0$. It is possible that this dependence on the problem size is an artifact of the technique of proof; without removing this dependence, large problem instances will require the measurement matrix to be a true isometry and the phase transition framework of the next section does not apply.

Phase Transitions for Greedy Algorithms with Gaussian Matrices

The quantities $\mu^{alg}(k,n,N)$ and $\xi^{alg}(k,n,N)$ in Theorems 5, 6, and 7 dictate the current theoretical convergence bounds for CoSaMP, SP, and IHT. Although some comparisons can be made between the forms of $\mu^{alg}$ and $\xi^{alg}$ for different algorithms, it is not possible to quantitatively state for what range of $k$ the algorithm will satisfy bounds on $\mu^{alg}(k,n,N)$ and $\xi^{alg}(k,n,N)$ for a specific value of $n$ and $N$. To establish quantitative interpretations of the conditions in Theorems 5, 6 and 7, it is necessary to have quantitative bounds on the behaviour of the aRIP constants $L(k,n,N)$ and $U(k,n,N)$ for the matrix $A$ in question, . Currently, there is no known matrix $A$ for which it has been proven that $U(k,n,N)$ and $L(k,n,N)$ remain bounded above and away from one, respectively, as $n$ grows, for $k$ and $N$ proportional to $n$. However, it is known that for some random matrix ensembles, with exponentially high probability on the draw of $A$, $\frac{1}{1-L(k,n,N)}$ and $U(k,n,N)$ do remain bounded as $n$ grows, for $k$ and $N$ proportional to $n$. The ensemble with the best known bounds on the growth rates of $L(k,n,N)$ and $U(k,n,N)$ in this setting is the Gaussian ensemble. In this section, we consider large problem sizes as $(k,n,N)\rightarrow\infty$, with $\frac{n}{N}\rightarrow\delta$ and $\frac{k}{n}\rightarrow\rho$ for $\delta,\rho\in(0,1)$. We study the implications of the sufficient conditions from Section 2 for matrices with Gaussian i.i.d. entries, namely, entries drawn i.i.d. from the normal distribution with mean and variance $n^{-1}$, $\mathcal{N}(0,n^{-1})$.

Gaussian random matrices are well studied and much is known about the behavior of their eigenvalues. Edelman derived bounds on the probability distribution functions of the largest and smallest eigenvalues of the Wishart matrices derived from a matrix $A$ with Gaussian i.i.d. entries. Select a subset of columns indexed by $I\subset\{1,\dots,N\}$ with cardinality $k$ and form the submatrix $A_{I}$, and the associated Wishart matrix derived from $A_{I}$ is the matrix $A_{I}^{*}A_{I}$. The distribution of the most extreme eigenvalues of all ${N\choose k}$ Wishart matrices derived from $A$ with Gaussian i.i.d. entries is only of recent interest and the exact probability distribution functions are not known. Recently, using Edelman’s bounds , the first three authors derived upper bounds on the probability distribution functions for the most extreme eigenvalues of all ${N\choose k}$ Wishart matrices derived from $A$. These bounds enabled them to formulate upper bounds on the aRIP constants, $L(k,n,N)$ and $U(k,n,N)$, for matrices of size $n\times N$ with Gaussian i.i.d. entries.

Let $A$ be a matrix of size $n\times N$ whose entries are drawn i.i.d. from ${\cal N}(0,n^{-1})$ and let $n\rightarrow\infty$ with $\frac{k}{n}\rightarrow\rho$ and $\frac{n}{N}\rightarrow\delta$. Let $H(p):=p\log(1/p)+(1-p)\log(1/(1-p))$ denote the usual Shannon Entropy with base $e$ logarithms, and let

Define $\lambda_{min}(\delta,\rho)$ and $\lambda_{max}(\delta,\rho)$ as the solution to (20) and (21), respectively:

Define $L(\delta,\rho)$ and $U(\delta,\rho)$ as

For any $\epsilon>0$, as $n\rightarrow\infty$,

The details of the proof of Theorem 9 are found in . The bounds are derived using a simple union bound over all $N\choose k$ of the $k\times k$ Wishart matrices $A_{I}^{*}A_{I}$ that can be formed from columns of $A$. Bounds on the tail behavior of the probability distribution function for the largest and smallest eigenvalues of $A_{I}^{*}A_{I}$ can be expressed in the form $p(n,\lambda)\exp(n\psi(\lambda,\rho))$ with $\psi$ defined in (18) and (19) and $p(n,\lambda)$ a polynomial. Following standard practices in large deviation analysis, the tails of the probability distribution functionals are balanced against the exponentially large number of Wishart matrices (20) and (21) to define upper and lower bounds on the largest and smallest eigenvalues of all $N\choose k$ Wishart matrices, with bounds $\lambda_{max}(\delta,\rho)$ and $\lambda_{min}(\delta,\rho)$, respectively. Overestimation of the union bound over the combinatorial number of $N\choose k$ Wishart matrices causes the bound $\lambda_{max}(\delta,\rho)$ to not be strictly increasing in $\rho$ for $\delta$ large; to utilize the best available bound on the extreme of the largest eigenvalue, we note that any bound $\lambda_{max}(\delta,\nu)$ for $\nu\in[\rho,1]$ is also a valid bound for submatrices of size $n\times k$. The asymptotic bounds of the aRIP constants, $L(\delta,\rho)$ and $U(\delta,\rho)$, follow directly. See Figure 3 for level curves of the bounds.

With Theorem 9, we are able to formulate quantitative statements about the matrices $A$ with Gaussian i.i.d. entries which satisfy the sufficient aRIP conditions from Section 2. A naive replacement of each $L(\cdot,n,N)$ and $U(\cdot,n,N)$ in Theorems 5-7 with the asymptotic aRIP bounds in Theorem 9 is valid in these cases. The properties necessary for this replacement are detailed in Lemma 16, stated in the Appendix. For each algorithm (CoSaMP, SP and IHT) the recovery conditions can be stated in the same format as Theorem 1 and Corollary 2, with only the expressions for $\kappa(\delta,\rho)$, $\mu(\delta,\rho)$ and $\xi(\delta,\rho)$ differing. These recovery factors are stated in Theorems 10-12.

Theorem 1 and Corollary 2 are satisfied for CoSaMP, Algorithm 1, with $\kappa^{csp}(\delta,\rho):=1$ and $\mu^{csp}(\delta,\rho)$ and $\xi^{csp}(\delta,\rho)$ defined as

The phase transition lower bound $\rho_{S}^{csp}(\delta)$ is defined as the solution to $\mu^{csp}(\delta,\rho)=1$. $\rho_{S}^{csp}(\delta)$ is displayed as the black curve in Figure 1(a). $\mu^{csp}(\delta,\rho)$ and $\xi^{csp}(\delta,\rho)/(1-\mu^{csp}(\delta,\rho))$ are displayed in Figure 2 panels (a) and (b) respectively.

Theorem 1 and Corollary 2 are satisfied for Subspace Pursuit, Algorithm 2, with $\kappa^{sp}(\delta,\rho)$, $\mu^{sp}(\delta,\rho)$, and $\xi^{sp}(\delta,\rho)$ defined as

The phase transition lower bound $\rho_{S}^{sp}(\delta)$ is defined as the solution to $\mu^{sp}(\delta,\rho)=1$. $\rho_{S}^{sp}(\delta)$ is displayed as the magenta curve in Figure 1(a). $\mu^{sp}(\delta,\rho)$ and $\xi^{sp}(\delta,\rho)/(1-\mu^{sp}(\delta,\rho))$ are displayed in Figure 2 panels (c) and (d) respectively.

Theorem 1 and Corollary 2 are satisfied for Iterative Hard Thresholding, Algorithm 3, with $\omega:=2/(2+U(\delta,3\rho)-L(\delta,3\rho))$, $\kappa^{iht}(\delta,\rho):=1$, and $\mu^{iht}(\delta,\rho)$ and $\xi^{iht}(\delta,\rho)$ defined as

The phase transition lower bound $\rho_{S}^{iht}(\delta)$ is defined as the solution to $\mu^{iht}(\delta,\rho)=1$. $\rho_{S}^{iht}(\delta)$ is displayed as the red curve in Figure 1(a). $\mu^{iht}(\delta,\rho)$ and $\xi^{iht}(\delta,\rho)/(1-\mu^{iht}(\delta,\rho))$ are displayed in Figure 2 panels (e) and (f) respectively.

Given a matrix $A$ with entries drawn i.i.d. from ${\cal N}(0,n^{-1})$, for any $x\in\chi^{N}(k)$, let $y=Ax+e$ for some (unknown) noise vector $e$. Define

Discussion and Conclusions

We have presented a framework in which recoverability results for sparse approximation algorithms derived using the ubiquitous RIP can be easily compared. This phase transition framework, , translates the generic RIP-based conditions of Theorem 3 into specific sparsity levels $k$ and problem sizes $n$ and $N$ for which the algorithm is guaranteed to satisfy the sufficient RIP conditions with high probability on the draw of the measurement matrix; see Theorem 1. Deriving (bounds on) the phase transitions requires bounds on the behaviour of the measurement matrix’ RIP constants . To achieve the most favorable quantitative bounds on the phase transitions, we used the less restrictive aRIP constants; moreover, we employed the best known bounds on aRIP constants, those provided for Gaussian matrices , see Theorem 9.

This framework was illustrated on three exemplar greedy algorithms: CoSaMP , SP , and IHT . The lower bounds on the phase transitions in Theorems 10-12 allow for a direct comparison of the current theoretical results/guarantees for these algorithms.

Computational Cost of CoSaMP, SP and IHT

The major computational cost per iteration in these algorithms is the application of one or more pseudoinverses. SP uses two pseudoinverses of dimensions $k\times n$ per iteration and another to compute the output vector $\hat{x}$; see Algorithm 2. CoSaMP uses only one pseudoinverse per iteration but of dimensions $2k\times n$; see Algorithm 1. Consequently, CoSaMP and SP have identical computational cost per iteration, of order $kn^{2}$, if the pseudoinverse is solved using an exact $QR$ factorization. IHT avoids computing a pseudoinverse altogether in internal iterations, but is aided by one pseudoinverse of dimensions $k\times n$ on the final support set. Thus IHT has a substantially lower computational cost per iteration than CoSaMP and SP. Note that pseudoinverses may be computed approximately by an iterative method such as conjugate gradients . As such, the exact application of a pseudoinverse could be entirely avoided, improving the implementation costs of these algorithms, especially of CoSaMP and SP.

Globally, all three algorithms converge linearly; in fact, they converge in a finite number of iterations provided there exists a $k$-sparse solution to $Ax=y$ and a sufficient aRIP condition is satisfied, see Corollary 2. For each algorithm, the upper bound on the required number of iterations grows unbounded as the function $\mu^{alg}(k,n,N)\rightarrow 1$. Hence, according to the bounds presented here, to ensure rapid convergence, it is advantageous to have a matrix that satisfies a more strict condition, such as $\mu^{alg}(k,n,N)<\frac{1}{2}$. Similarly, the factor controlling stability to additive noise, namely the vector $e$ in Theorem 1, blows up as the function $\mu^{alg}(k,n,N)\rightarrow 1$. Again, according to the bounds presented here, in order to guarantee stability with small amplification of the additive noise, it is necessary to restrict the range of $\frac{\xi^{alg}}{1-\mu^{alg}}(k,n,N)$. A phase transition function analogous to the functions $\rho_{S}^{alg}(\delta)$ can be easily computed in these settings as well, resulting in curves lower than those presented in Figure 1(a). This is the standard trade-off of compressed sensing, where one must determine the appropriate balance between computational efficiency, stability, and minimizing the number of measurements.

Comparison of Phase Transitions and Constants of Proportionality

From Figure 1(a), we see that the best known lower bounds on the phase transitions for the three greedy algorithms satisfy the ordering $\rho_{S}^{csp}(\delta)<\rho_{S}^{sp}(\delta)<\rho_{S}^{iht}(\delta)$ for Gaussian measurement matrices. Therefore, we now know that, at least for Gaussian matrices, according to existing thoery, IHT has the largest region where recovery for all signals can be guaranteed; the regions with similar guarantees for SP and CoSaMP are considerably smaller. Moreover, IHT has a lower bound on its computational cost.

Future Improvements and Conclusions

The above bounds on greedy algorithms’ phase transitions could be improved by further refining the algorithms’ theory, namely, deriving less strict aRIP conditions on the measurement matrix that still ensure convergence of the algorithm; as the latter is an active research topic, we expect such developments to take place. The phase transition framework presented here may also be applied to such advances. Alternatively, increasing the lower bounds on the phase transitions could be expected to occur from improving the upper bounds we employed on the aRIP constants of the Gaussian measurement matrices, see Theorem 9. However, extensive empirical calculations of lower estimates of aRIP constants show the latter to be within a factor of $1.83$ of our proven upper bounds . During the revision of this manuscript, improved bounds on the aRIP constants for the Gaussian ensemble were derived , tightening the bound to be within $1.57$ of lower estimates. However, for $\rho\approx 10^{-3}$ both bounds were already very sharp , and the resulting increase of the phase transitions shown here was under 0.5%.

Appendix A Proofs of Main Results

We present a framework by which RIP-based convergence results of the form presented in Theorem 3 can be translated into results of the form of Theorem 1; that is removing explicit dependencies on RIP constants in favour of their bounds.

The proofs of Theorems 5, 6, and 7 rely heavily on a sequence of properties of the aRIP constants, which are summarize in Lemma 15 and proven in Section A.1. Theorems 10, 11, and 12 follow from Theorems 5, 6, and 7 and the form of $\mu^{alg}$ and $\xi^{alg}$ as functions of $L$ and $U$; this latter point is summarized in Lemma 16 which is stated and proven in Section A.1. The resulting Theorems 10, 11, and 12 can then be interpreted in the phase transition framework advocated by Donoho et al. , as we have explained in Section 4.

Throughout the analysis of the algorithms, we repeatedly use implications of the aRIP on a matrix $A$ as outlined in Lemma 15. This lemma has been proven in the symmetric case repeatedly in the literature; we include the proof of the asymmetric variant for completeness.

$\left\|A_{I}^{*}y\right\|_{2}\leq\sqrt{1+U(|I|,n,N)}\|y\|_{2}$

$(1-L(|I|,n,N))\|u\|_{2}\leq\left\|A_{I}^{*}A_{I}u\right\|_{2}\leq(1+U(|I|,n,N))\|u\|_{2}$

$\left\|A_{I}^{\dagger}y\right\|_{2}\leq(1-L(|I|,n,N))^{-\frac{1}{2}}\|y\|_{2}$

$\left|\left\langle A_{I}u,A_{J}v\right\rangle\right|\leq\frac{1}{2}\left(L(|I|+|J|,n,N)+U(|I|+|J|,n,N)\right)\|u\|_{2}\|v\|_{2}$

$\left\|A_{I}^{*}A_{J}v\right\|_{2}\leq U(|I|+|J|,n,N)\|v\|_{2}$.

$\left\|(Id-\omega A_{I}^{*}A_{I})u\right\|_{2}\leq\max\left\{\omega(1+U(|I|,n,N))-1,1-\omega(1-L(|I|,n,N))\right\}\|u\|_{2}.$

From Remark 2, it is clear that the aRIP constants are nondecreasing in the first argument pertaining the sparsity level. Therefore, $A$ must also have the aRIP constants $L(|I|,n,N)\leq L(|I|+|J|,n,N)$ and $U(|I|,n,N)\leq U(|I|+|J|,n,N)$. Also, Definition 1 implies that the singular values of the submatrix $A_{I}$ are contained in the interval $[\sqrt{1-L(|I|,n,N)},\sqrt{1+U(|I|,n,N)}]$. Thus, (i)-(iii) follow from the standard relationships between the singular values of $A_{I}$ and the associated matrix in (i)-(iii).

To prove (iv), let $m=|I|+|J|$. We may assume $\left\|u\right\|_{2}=\|v\|_{2}=1$; otherwise we normalize the vectors. Let $\alpha=A_{I}u$ and $\beta=A_{J}v$. Then, since $I\cap J=\emptyset$,

$\left[A_{I},A_{J}\right]$ is a submatrix of $A$ of size $n\times m$, so applying Definition 1 to the right most portion of (35) and invoking (38), we have

Thus $|\left\langle A_{I}u,A_{J}v\right\rangle|=|\left\langle\alpha,\beta\right\rangle|\leq\left(L(m,n,N)+U(m,n,N)\right)/2$, establishing (iv).

Since $I\cap J=\emptyset$, the matrix $-A_{I}^{*}A_{J}$ is a submatrix of $Id_{I\cup J}-A_{I}^{*}A_{J}$. To establish (v), we observe that the aRIP implies that the eigenvalues of every size $n\times m$ submatrix of $A$ lie in the interval $[1-L(m,n,N),1+U(m,n,N)]$. Thus the eigenvalues of $Id_{I\cup J}-A_{I}^{*}A_{J}$ must lie in $[0,U(m,n,N)]$. Therefore $\left\|A_{I}^{*}A_{J}v\right\|_{2}^{2}=\left\|-A_{I}^{*}A_{J}v\right\|_{2}^{2}$ completes the proof of (v).

To prove (vi), note that $\|(Id-\omega A_{I}^{*}A_{I})\|_{2}$ is bounded above by the maximum magnitude of the eigenvalues of $(\omega A_{I}^{*}A_{I}-Id)$, which lie in the interval with endpoints $\omega(1-L(|I|,n,N))-1$ and $\omega(1+U(|I|,n,N))-1$.

Theorems 10, 11, and 12 follow from Theorems 5, 6, and 7 and the form of $\mu^{alg}$ and $\xi^{alg}$ as functions of $L$ and $U$. We formalize the relevant functional dependencies in the next three lemmas.

Suppose, for all $t\in\mathcal{Z}$, $\left(\nabla F[t]\right)_{i}\geq 0$ for all $i=1,\dots,p+q$ and for any $v\in\mathcal{Z}$ we have $\nabla F[t]\cdot v>0$. Then for any $c\epsilon>0$, as $(k,n,N)\rightarrow\infty$ with $\frac{n}{N}\rightarrow\delta,\frac{k}{n}\rightarrow\rho$, there is an exponentially high probability on the draw of the matrix $A$ that

Suppose, for all $t\in\mathcal{Z}$, $\left(\nabla F[t]\right)_{i}\geq 0$ for all $i=1,\dots,p+q$ and there exists $j\in\{1,\dots,p\}$ such that $\left(\nabla F[t]\right)_{j}>0$. Then there exists $c\in(0,1)$ depending only on $F,\delta$,and $\rho$ such that for any $\epsilon\in(0,1)$

and so there is an exponentially high probability on the draw of $A$ that

Also, $F(z(\delta,\rho))$ is strictly increasing in $\rho$.

To prove (i), suppose $u,v\in\mathcal{Z}$ with $v_{i}>u_{i}$ for all $i=1,\dots,p+q$. From Taylor’s Theorem, $F[v]=F[u+(v-u)]=F[u]+\nabla F[t]\cdot[v-u]$ with $t=u+\lambda[v-u]$ for some $\lambda\in(0,1)$. Then

since, by assumption, $\nabla F[t]\cdot[v-u]>0$.

From Theorem 9, for any $c\epsilon>0$ and any $i=1,\dots,p+q$, as $(k,n,N)\rightarrow\infty$ with $\frac{n}{N}\rightarrow\delta$, $\frac{k}{n}\rightarrow\rho$,

with convergence to $1$ exponential in $n$. Therefore, letting $v_{i}:=z(\delta,\rho)_{i}+c\epsilon$ and $u_{i}:=z(k,n,N)_{i}$, for all $i=1,\dots,p+q$, we conclude from (45) that

again with convergence to $1$ exponential in $n$.

To establish (ii), we take the Taylor expansion of $F$ centered at $z(\delta,\rho)$, namely

Since $L(\delta,\rho)$ is strictly increasing in $\rho$ , then $\left(\left.\frac{\partial}{\partial\rho}z(\delta,\rho)\right|_{\rho=t_{2}^{\star}}\right)_{j}>0$ for all $j=1,\dots,p$. Since $U(\delta,\rho)$ is nondecreasing in $\rho$ , then $\left(\left.\frac{\partial}{\partial\rho}z(\delta,\rho)\right|_{\rho=t_{2}^{\star}}\right)_{i}\geq 0$ for all $i=p+1,\dots,p+q$. Hence, by the hypotheses of (ii),

(48) and (49) imply (43). Since the hypotheses of (ii) imply those of (i), (42) also holds, and so (44) follows. $F(z(\delta,\rho))$ strictly increasing follows from the hypotheses of (ii) and $L(\delta,\rho)$ and $U(\delta,\rho)$ strictly increasing and nondecreasing in $\rho$, respectively .

Assume that $\mu^{alg}(\delta,\rho)$ is strictly increasing in $\rho$ and let $\rho_{S}^{alg}(\delta)$ solve $\mu^{alg}(\delta,\rho)=1$. For any $\epsilon\in(0,1)$, if $\rho<(1-\epsilon)\rho_{S}^{alg}(\delta)$, then $\mu^{alg}(\delta,(1+\epsilon)\rho)<1$.

Let $\rho^{alg}_{\epsilon}(\delta)$ be the solution to $\mu^{alg}(\delta,(1+\epsilon)\rho)=1$. Since by definition, $\rho_{S}^{alg}(\delta)$ denotes a solution to $\mu^{alg}(\delta,\rho)=1$, and this solution is unique as $\mu^{alg}(\delta,\rho)$ is strictly increasing, we must have $(1+\epsilon)\rho^{alg}_{\epsilon}(\delta)=\rho_{S}^{alg}(\delta)$. Since $(1-\epsilon)<(1+\epsilon)^{-1}$ for all $\epsilon\in(0,1)$, we have $(1-\epsilon)\rho_{S}^{alg}(\delta)<\rho_{\epsilon}^{alg}(\delta)$. If $\rho<(1-\epsilon)\rho_{S}^{alg}(\delta)$, then since $\mu^{alg}(\delta,\rho)$ is strictly increasing in $\rho$,

Note that Lemma 16 ii) with $F:=\mu^{alg}$ will be employed to show the first assumption in Lemma 17; this is but one of several good uses of Lemma 16 that we will make.

Corollaries 2 and 4 are easily derived from Lemma 18. Note that this lemma demonstrates only that the support set has been recovered. The proof of Lemma 18 is a minor generalization of a proof from [12, Theorem 7].

Suppose, after $l$ iterations, algorithm $alg$ returns the $k$-sparse approximation $\hat{x}^{l}$ to a $k$-sparse target signal $x$. Suppose there exist constants $\mu$ and $\kappa$ independent of $l$ and $x$ such that

where $\nu_{\text{min}}(x)$ is defined in (5).

Theorems 5, 6 and 7 define the constants $\mu=\mu^{alg}(k,n,N)$ and $\kappa$ to be used in Lemma 18 for proving Corollary 4. For CoSaMP and IHT, $\kappa=\kappa^{alg}(k,n,N)=1$. For SP, the term involving $\kappa$ is removed by combining Lemmas 23 and 24 (with $e=0$) to obtain

which again, gives $\kappa=1$. Similarly, Theorems 10, 11 and 12 define the constants $\mu=\mu^{alg}(\delta,\rho)$ and $\kappa$ to be used in Lemma 18 for proving Corollary 2, with the above comments on the IHT choice of $\kappa$ also applying in this case.

To ensure exact recovery of the target signal, namely, to complete the proof of Corollaries 2 and 4, we actually need something stronger than recovering the support set as implied by Lemma 18. For CoSaMP and SP, since the algorithms employ a pseudoinverse at an appropriate step, the output is then the exact sparse signal. For IHT, no pseudoinverse has been applied; thus, to recover the signal exactly, one simply determines $T$ from the output vector and then $x=A_{T}^{\dagger}y$. These comments and Lemma 18 now establish Corollaries 2 and 4 for each algorithm, and we will not restate the proof for each individual algorithm.

In each of the following subsections, we first consider the case of general measurement matrices, $A$, and prove the results from Section 2 which establish an aRIP condition for an algorithm. We then proceed to choose a specific matrix ensemble, matrices with Gaussian i.i.d. entries, for which Section 3 establishes lower bounds on the phase transition for exact recovery of all $x\in\chi^{N}(k)$ and then provide probabilistic bounds on the multiplicative stability factors.

A.2 Proofs for CoSaMP

In this section we prove the results from Sections 2 and 3 reported for CoSaMP . The proofs mimic those of Needell and Tropp while employing the aRIP constants. In each proof, the smallest possible support is retained for the aRIP constants in order to acquire from this method of analysis the best possible conditions on the measurement matrix used in the CoSaMP algorithm. This change is in many cases straightforward, requiring only a substitution of $U(ak,n,N)$ or $L(ak,n,N)$ for $R(ak,n,N)$, for some $a\in\{2,3,4\}$. In such cases we simply restate the result. Where there is a more substantial change, we provide fuller details of the proof.

The preceding lemmas facilitate the proof of Theorem 5.

Given our assumption that $\mu^{csp}(k,n,N)<1$, we can now prove a stronger statement, namely that for $l\geq 0$ we have

We proceed by induction. Assume the result holds for some $l\geq 0$. Then, applying the inductive hypothesis and (54), we have

and so the result is also true for $l+1$, and so (55) holds for all $l\geq 0$ by induction. Finally, note that

Having established the results of Section 2 for CoSaMP, we now focus on Gaussian random matrices and prove the results from Section 3 concerning CoSAMP.

Let $x,y,A$ and $e$ satisfy the hypothesis of Theorem 10 and select $\epsilon>0$. Fix $\tau<1$ and let

Clearly, $\left(\nabla F^{csp}[t]\right)_{i}\geq 0$ for all $i=1,\dots,5$ and

Hence the hypotheses of Lemma 16 (ii) are satisfied for $F^{csp}$. By (10), (23) and (56), $F^{csp}[z(k,n,N)]=\mu^{csp}(k,n,N)$ and $F^{csp}[z(\delta,\rho)]=\mu^{csp}(\delta,\rho)$. Thus, by Lemma 16, as $(k,n,N)\rightarrow\infty$ with $\frac{n}{N}\rightarrow\delta$, $\frac{k}{n}\rightarrow\rho$,

Also, $\mu^{csp}(\delta,\rho)$ is strictly increasing in $\rho$ and so Lemma 17 applies.

Similarly, $G^{csp}$ satisfies the hypotheses of Lemma 16 (ii). Likewise, by (11), (24) and (57), $G^{csp}[z(k,n,N)]=\xi^{csp}(k,n,N)$ and $G^{csp}[z(\delta,\rho)]=\xi^{csp}(\delta,\rho)$. Again, by Lemma 16, as $(k,n,N)\rightarrow\infty$ with $\frac{n}{N}\rightarrow\delta$, $\frac{k}{n}\rightarrow\rho$,

Therefore, for any $x\in\chi^{N}(k)$ and any noise vector $e$, as $(k,n,N)\rightarrow\infty$ with $\frac{n}{N}\rightarrow\delta$, $\frac{k}{n}\rightarrow\rho$, there is an exponentially high probability on the draw of a matrix $A$ with Gaussian i.i.d. entries that

Combining (60) with Theorem 5 completes the argument.

A.3 Proofs for Subspace Pursuit

In this section we outline the proofs for the results in Section 2 and then prove the results in Section 3 reported for SP . The proofs mimic those of Dai and Milenkovic while employing the aRIP constants. In each proof, the smallest possible support is retained for the aRIP constants in order to acquire from this method of analysis the best possible conditions on the measurement matrix used in the SP algorithm.

We begin in the setting of an arbitrary measurement matrix $A$ of size $n\times N$ and formulate the aRIP conditions of Theorem 6. A sequence of lemmas leads us to Theorem 6. Lemmas 23 and 24 directly follow the proofs from [12, Theorem 10] with the adaptation that we employ the aRIP constants from Definition 1, Lemma 15, and we maintain the smallest support size in $L(\cdot,n,N),U(\cdot,n,N)$.

For $x\in\chi^{N}(k)$ and $y=Ax+e$, after iteration $l$ of SP

For $x\in\chi^{N}(k)$ and $y=Ax+e$, after iteration $l$ of SP

The following lemma is an adaptation of [12, Lemma 3]. By using Definition 1 and selecting the smallest possible support sizes for the aRIP constants, we arrive at Lemma 25.

Let $x\in\chi^{N}(k)$ and $y=Ax+e$ be the measurement contaminated with noise $e$. If the Subspace Pursuit algorithm terminates after $l$ iterations, the output $\hat{x}$ approximates $x$ within the bounds

After applying Lemma 23 to Lemma 24 we bound the entries of $x$ that have not been captured by Algorithm 2, namely

Applying (64) iteratively, we develop a bound in terms of the norm of $x$, by observing that $\|x_{T-T^{0}}\|_{2}\leq\|x\|_{2}$:

The factor $\frac{\phi^{sp}(k,n,N)}{1-\mu^{sp}(k,n,N)}$ amplifying $\|e\|_{2}$ in (66) is found by induction as in the proof of Theorem 5 in Appendix A.2.

From Lemma 25 with $\kappa^{sp}(k,n,N)=1+\frac{U(2k,n,N)}{1-L(k,n,N)}$, we have

Having established the aRIP conditions for an arbitrary measurement matrix, we again return to the Gaussian random matrix ensemble and establish the quantitative bounds for SP from Section 3.

Let $x,y,A,$ and $e$ satisfy the hypothesis of Theorem 11 and select $\epsilon>0$. Fix $\tau<1$ and let

For each of these functions, the gradient is clearly nonnegative componentwise on $\mathcal{Z}$, with the first entry of each gradient strictly positive which is sufficient to verify the hypotheses of Lemma 16 (ii). Moreover, from (12)–(14) and (25)–(27), we have

Invoking Lemma 16 for each of the functions in (70)–(73) yields that with high probability on the draw of $A$ from a Gaussian distribution,

Combining (74) and (75) with Theorem 6 completes the argument, recalling that Lemma 16 applied to $F^{sp}=\mu^{sp}$ also implies that $\mu^{sp}(\delta,\rho)$ is strictly increasing in $\rho$ and so Lemma 17 holds.

A.4 Proofs for Iterative Hard Thresholding

In this section we first outline a proof of Theorem 7, which follows similar lines to that given by Blumensath and Davies in [7, Corollary 4], while considering a generalization to aRIP bounds, and also incorporating a stepsize $\omega$. Having established this result for arbitrary measurement matrices, we then go on to prove Theorem 12 which gives conditions for high-probability convergence of IHT in the specific case of Gaussian random matrices.

Let $B^{l}=T^{l}\cup\mbox{supp}(x)$. Since $|B^{l}|\leq 2k\leq 3k$, we can deduce from Lemma 15 that

where $\phi^{iht}(3k,n,N)$ is defined to be

Since the eigenvalues of a submatrix are bounded in magnitude by the eigenvalues of the entire matrix, and since $|B^{l}\cup B^{l-1}|\leq 3k$, we can again invoke Lemma 15 to obtain

Now we have from the proof of [7, Corollary 4] that

Substituting (76) and (77) into (78), and applying Lemma 15 to the error term, we obtain

Now $B^{l}$ and $(B^{l-1}-B^{l})$ are disjoint, so we have

with $\mu^{iht}(k,n,N)$ and $\xi^{iht}(k,n,N)$ defined in (15) and (16), respectively. Given our assumption that $\mu^{iht}(k,n,N)<1$, an induction argument analogous to the induction in the proof of Theorem 5 gives the stronger result

We finally note that if IHT terminates after $l$ iterations we have $\hat{x}=x^{l}_{T^{l}}$, from which the results now follows.

Armed with the results of Section 2 for IHT, we return to the family of Gaussian random matrices and prove the quantitative bounds for IHT from Section 3.

Let $x,y,A$ and $e$ satisfy the hypothesis of Theorem 12 and select $\epsilon>0$. Fix $\tau<1$ and let

[Note that $F^{iht}_{\omega}[z(k,n,N)]=\mu^{iht}(k,n,N)$ and $G^{iht}_{\omega}[z(k,n,N)]=\xi^{iht}(k,n,N)/(1-\mu^{iht}(k,n,N))$ due to (15) and (16).] Clearly the functions are nondecreasing so that, with any $t\in\mathcal{Z}$, $\left(\nabla F^{iht}_{\omega}[t]\right)_{i}\geq 0$ and $\left(\nabla G^{iht}_{\omega}[t]\right)_{i}\geq 0$ for $i=1,2,3$; note that $F^{iht}_{\omega}[t]$ and $G^{iht}_{\omega}[t]$ have points of nondifferentiability, but that the left and right derivatives at those points remain nonnegative. Also, and for any $v\in\mathcal{Z}$, since $t_{i},v_{i}>0$ for each $i$, $\nabla F^{iht}_{\omega}[t]\cdot v>0$ and $\nabla G^{iht}_{\omega}[t]\cdot v>0$ as both functions clearly increase when each component of the argument increases. Hence, $F^{iht}_{\omega}$ and $G^{iht}_{\omega}$ satisfy the hypotheses of Lemma 16 (i). Therefore, for any $\omega\in(0,1)$, as $(k,n,N)\rightarrow\infty$ with $\frac{n}{N}\rightarrow\delta$, $\frac{k}{n}\rightarrow\rho$,

Now fix $\omega^{\star}:=\frac{2}{2+U(\delta,3\rho)-L(\delta,3\rho)}$ and define

In (81) and (82), the weight was arbitrary; thus both statements certainly hold for the particular weight $\omega^{\star}$. Therefore, combining (81), (85), (87) and combining (82), (86), (88) imply that with exponentially high probability on the draw of $A$,

Therefore, with the weight $\omega^{\star}$, there is an exponentially high probability on the draw of $A$ from a Gaussian distribution that

References