Generalizing CoSaMP to Signals from a Union of Low Dimensional Linear Subspaces

TL;DR

This paper extends the CoSaMP algorithm to recover signals from a union of low-dimensional subspaces, providing theoretical guarantees and a new greedy reconstruction method, supported by experimental validation.

Contribution

It generalizes CoSaMP for union of subspaces using Gaussian mean width and introduces a new greedy algorithm inspired by this theory.

Findings

The generalized CoSaMP guarantees recovery under Gaussian measurements.

The new greedy method performs well in combined sparse-synthesis and cosparse-analysis models.

Experimental results validate the effectiveness of the proposed approach.

Abstract

The idea that signals reside in a union of low dimensional subspaces subsumes many low dimensional models that have been used extensively in the recent decade in many fields and applications. Until recently, the vast majority of works have studied each one of these models on its own. However, a recent approach suggests providing general theory for low dimensional models using their Gaussian mean width, which serves as a measure for the intrinsic low dimensionality of the data. In this work we use this novel approach to study a generalized version of the popular compressive sampling matching pursuit (CoSaMP) algorithm, and to provide general recovery guarantees for signals from a union of low dimensional linear subspaces, under the assumption that the measurement matrix is Gaussian. We discuss the implications of our results for specific models, and use the generalized algorithm as an inspiration for a new greedy method for signal reconstruction in a combined sparse-synthesis and cosparse-analysis model. We perform experiments that demonstrate the usefulness of the proposed strategy.

Figures

Equations

$$\mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{e},$$

$$\Vert \hat{\mathbf{x}}-\mathbf{x}{\Vert }_2\le c_0+c_1\Vert \mathbf{e}{\Vert }_2$$

$${\mathcal{K}}_1+{\mathcal{K}}_2\triangleq \left\{{\mathbf{x}}_1+{\mathbf{x}}_2:{\mathbf{x}}_1\in {\mathcal{K}}_1,{\mathbf{x}}_2\in {\mathcal{K}}_2\right\}.$$

$${\mathcal{S}}^B\triangleq \left\{\sum _{i=1}^B{\mathcal{V}}_i:{\mathcal{V}}_i\in \mathcal{S}\right\}.$$

$${\mathcal{U}}^B\triangleq \left\{\bigcup _i{\mathcal{V}}_i^B:{\mathcal{V}}_i^B\in {\mathcal{S}}^B\right\}.$$

$${\mathcal{U}}^B=\sum _{i=1}^B{\mathcal{U}}^1=\left\{\sum _{i=1}^B{\mathbf{x}}_i:{\mathbf{x}}_i\in {\mathcal{U}}^1\right\}.$$

$$w(\mathcal{K})\triangleq E_{\mathbf{g}}\left\{\underset{\mathbf{z}\in \mathcal{K}}{\sup }⟨\mathbf{g},\mathbf{z}⟩\right\},$$

$$b_m=E\left\{\Vert \mathbf{g}{\Vert }_2\right\},$$

$$\frac{\Vert \mathbf{A}\mathbf{u}{\Vert }_2}{\Vert \mathbf{u}{\Vert }_2}\ge b_m-w(\mathcal{K}\cap {\mathbb{S}}^{n-1})-\eta$$

$$\frac{\Vert \mathbf{A}\mathbf{u}{\Vert }_2}{\Vert \mathbf{u}{\Vert }_2}\le b_m+w(\mathcal{K}\cap {\mathbb{S}}^{n-1})+\eta$$

$$\Vert {\mathbf{P}}_{\mathcal{V}}({\mathbf{I}}_n-\mu {\mathbf{A}}^{\ast }\mathbf{A}){\mathbf{P}}_{\mathcal{V}}\Vert \le 1-\mu {\left[b_m-w\left({\mathcal{U}}^B\cap {\mathbb{S}}^{n-1}\right)-\eta \right]}^2$$

$$\frac{{\Vert \mathbf{A}{\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2}{{\Vert {\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2}\ge {\left[b_m-w({\mathcal{U}}^B\cap {\mathbb{S}}^{n-1})-\eta \right]}^2$$

$$\frac{{\Vert \mathbf{A}{\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2}{{\Vert {\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2}\le {\left[b_m+w({\mathcal{U}}^B\cap {\mathbb{S}}^{n-1})+\eta \right]}^2$$

$$1-\mu {\left[b_m-w\left({\mathcal{U}}^B\cap {\mathbb{S}}^{n-1}\right)-\eta \right]}^2$$

$$=\underset{\mathbf{z}:\Vert \mathbf{z}{\Vert }_2=1}{\sup }{\Vert {\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^{-2}\left({\Vert {\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2-\mu {\Vert \mathbf{A}{\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2\right)$$

$$\ge \underset{\mathbf{z}:\Vert \mathbf{z}{\Vert }_2=1}{\sup }\left({\Vert {\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2-\mu {\Vert \mathbf{A}{\mathbf{P}}_{\mathcal{V}}\mathbf{z}\Vert }_2^2\right)$$

$$=\underset{\mathbf{z}:\Vert \mathbf{z}{\Vert }_2=1}{\sup }⟨\mathbf{z},{\mathbf{P}}_{\mathcal{V}}({\mathbf{I}}_n-\mu {\mathbf{A}}^{\ast }\mathbf{A}){\mathbf{P}}_{\mathcal{V}}\mathbf{z}⟩$$

$$=\Vert {\mathbf{P}}_{\mathcal{V}}({\mathbf{I}}_n-\mu {\mathbf{A}}^{\ast }\mathbf{A}){\mathbf{P}}_{\mathcal{V}}\Vert .$$

$$\hat{\mathcal{V}}=\underset{\mathcal{V}\in {\mathcal{S}}^B}{\mathrm{argmin}}\Vert \tilde{\mathbf{x}}-{\mathbf{P}}_{\mathcal{V}}\tilde{\mathbf{x}}{\Vert }_2.$$

$$\Vert {\mathbf{x}}^t-\mathbf{x}{\Vert }_2\le \frac{4{\rho }_2(\eta )}{\sqrt{1-{\rho }_1^2(\eta )}}\Vert {\mathbf{x}}^{t-1}-\mathbf{x}{\Vert }_2+\left(\frac{4{\xi }_2(\eta )}{\sqrt{1-{\rho }_1^2(\eta )}}+\frac{2{\xi }_1(\eta )}{1-{\rho }_1(\eta )}\right)\Vert \mathbf{e}{\Vert }_2,$$

$${\rho }_1(\eta )$$

$${\xi }_1(\eta )$$

$${\rho }_2(\eta )$$

$${\xi }_2(\eta )$$

$${\mu }_1$$

$${\mu }_2$$

$$\Vert {\mathbf{x}}^t-\mathbf{x}{\Vert }_2\le 2\Vert {\tilde{\mathbf{x}}}^t-\mathbf{x}{\Vert }_2.$$

$$\Vert {\tilde{\mathbf{x}}}^t-\mathbf{x}{\Vert }_2\le \frac{1}{\sqrt{1-{\rho }_1^2(\eta )}}\Vert {\mathbf{Q}}_{{\tilde{\mathcal{V}}}^t}({\tilde{\mathbf{x}}}^t-\mathbf{x}){\Vert }_2+\frac{{\xi }_1(\eta )}{1-{\rho }_1(\eta )}\Vert \mathbf{e}{\Vert }_2$$

$$\Vert {\mathbf{Q}}_{{\tilde{\mathcal{V}}}^t}({\tilde{\mathbf{x}}}^t-\mathbf{x}){\Vert }_2\le 2{\rho }_2(\eta )\Vert \mathbf{x}-{\mathbf{x}}^{t-1}{\Vert }_2+2{\xi }_2(\eta )\Vert \mathbf{e}{\Vert }_2$$

$$b_m\ge \frac{m}{\sqrt{m+1}}=\sqrt{\frac{m}{1+\frac{1}{m}}}\ge \sqrt{\frac{14.5^2{m}_0+1}{1+\frac{1}{m}}}>\sqrt{\frac{14.5^2{m}_0+\frac{1}{m}14.5^2{m}_0}{1+\frac{1}{m}}}=14.5\sqrt{m_0},$$

$$16{\rho }_2^2(\eta )+{\rho }_1^2(\eta )$$

$$\Vert {\mathbf{x}}^t-\mathbf{x}{\Vert }_2\le {\rho }^t(m)\Vert {\mathbf{x}}^0-\mathbf{x}{\Vert }_2+\frac{\xi (m)}{1-\rho (m)}\Vert \mathbf{e}{\Vert }_2,$$

$$\rho (m)$$

$$\xi (m)$$

$$\Vert {\mathbf{x}}^t-\mathbf{x}{\Vert }_2\le {\tilde{\rho }}^t\Vert {\mathbf{x}}^0-\mathbf{x}{\Vert }_2+\frac{1-{\tilde{\rho }}^t}{1-\tilde{\rho }}\tilde{\xi }\Vert \mathbf{e}{\Vert }_2\le {\tilde{\rho }}^t\Vert {\mathbf{x}}^0-\mathbf{x}{\Vert }_2+\frac{\tilde{\xi }}{1-\tilde{\rho }}\Vert \mathbf{e}{\Vert }_2,$$

Full text

I Introduction

In the recent decade many signal and image processing applications have benefited remarkably from adopting low dimensional models for the signals of interest [1, 2, 3, 4, 5, 6]. The most popular form of this low dimensional modeling assumes that signals admit a sparse representation in a given dictionary - this is referred to as the "sparsity model" [1]. Other options for low dimensional modeling include the analysis cosparse framework [2], the low-rank model [3, 4], low dimensional manifolds [5] and Gaussian mixture models [6]. For each of them, vast amounts of algorithms and theory were developed [7, 8, 9]. The common factor of all these models is that they characterize data with a number of parameters, which is significantly smaller than the ambient dimension of the data.

Recently, it was suggested by several authors that the intrinsic low dimensionality of the data can be measured by its Gaussian mean width [10, 11], which is strongly related to the statistical dimension [12]. This new view enables analyzing low dimensional data that do not have an explicit representation and to unify existing results, providing general theory and algorithms for signals residing in low dimensional models.

In this work we consider the problem of recovering an unknown signal $\bm{x}\in\mathbb{R}^{n}$ from $m$ linear noisy measurements of the form

[TABLE]

where $\bm{A}$ is an $m\times n$ matrix with entries independently drawn from the standard normal distribution and $\bm{e}\in\mathbb{R}^{m}$ represents the noise. We assume that the noise does not depend on $\bm{A}$, and that its energy is bounded. Further, we assume that $m<n$. As in this case, recovering $\bm{x}$ from $\bm{y}$ becomes an ill-posed problem, an additional prior on $\bm{x}$ is required. For example, in compressive sampling [13, 14] the signal is assumed to admit a sparse representation, which allows its stable recovery from the measurements. In this work we rely on a more general signal model, a union of low dimensional linear subspaces. More specifically, we assume that a possibly infinite set of finite-dimensional subspaces is given, and that $\bm{x}$ resides in one of these subspaces. However, we do not know a priori in which one. This model subsumes many of the priors that were mentioned above, such as sparse representation; and more.

In the context of compressive sampling, two major strategies are used to recover a sparse signal from incomplete and possibly inaccurate samples. The first approach uses convex relaxation to formulate a convex program whose minimizer is known to approximate the target signal [15, 16, 17, 18]. The second one uses greedy methods to iteratively build an approximation, making locally optimal choices at each iteration. Examples of greedy programs include matching pursuit (MP) [19], orthogonal matching pursuit (OMP) [20, 21], regularized orthogonal matching pursuit (ROMP) [22], iterative hard thresholding (IHT) [23], subspace pursuit (SP) [24] and compressive sampling matching pursuit (CoSaMP) [25]. Recovery guarantees for the methods above (of both strategies) are mostly proved using the assumption that the measurement matrix $\bm{A}$ has the restricted isometry property (RIP) [16], with a sufficiently small RIP-constant. In general, convex relaxation methods have better recovery guarantees than greedy methods, but they are also more computationally demanding. However, IHT, SP, and in particular CoSaMP, were shown to have recovery guarantees that are competitive with those of the convex relaxation methods.

In this work we focus on the greedy strategy for a general union of low dimensional subspaces. We study a generalized version of CoSaMP, and provide recovery guarantees that depend on the Gaussian mean width of a union of subspaces, which is strongly related to the one in which $\bm{x}$ resides. We note that extensions of CoSaMP to different models exist [26, 27, 28, 29, 30]. However, most of these works study specific models, while here we consider a general framework that includes these models and more, as discussed in Section IV. Moreover, Our proof technique, and assumptions on $\bm{A}$ and $\bm{e}$, facilitate obtaining bounds on the reconstruction error that are more robust to noise than previous results, which are based on RIP and its generalizations.

The setting where $\bm{x}$ lies in a union of linear subspaces was studied in several works. Conditions for unique and stable sampling were developed in [31] for unions of infinitely many finite-dimensional subspaces, and in [32] for finite unions of finite-dimensional subspaces. In [33] the sampling theory was extended to models beyond unions of subspaces. However, all these works do not provide concrete recovery algorithms. Two prominent works that use results from [32] and propose recovery algorithms are [34] and [26]. In [34] the modeling of each subspace in the union is rather general (though, less general than in this paper), but the number of possible subspaces is assumed to be finite and known, and the proposed recovery algorithm, which is based on convex relaxation, deeply depends on this assumption. In [26], similarly to the work we present here, a generalized version of CoSaMP was studied. However, the signal model in [26] is restricted to predefined sparsity patterns under a predefined basis, thus the number of possible subspaces is known and the guarantees do not cover other models (such as sparsity under a given redundant dictionary). The model that we consider here does not assume that the union is finite. Hence, it includes specific models that are not included in [34, 26]. One such example is a signal that can be formed as a rank-$k$ matrix. Perhaps the most related work to the current paper is [35], which also considers reconstruction of signals from a possibly infinite union of linear subspaces. However, [35] study the projected Landweber algorithm, which is basically a generalization of IHT to models beyond sparsity. Here we build on CoSaMP, which is stronger than IHT in general. In addition, our analysis shows more robustness to noise than the analysis in [35], which is based on assumptions that can be seen as a generalization of RIP.

Contribution. In this work we make the following contributions. We introduce a generalized CoSaMP (GCoSaMP) algorithm that relies on a general signal model - union of low dimensional subspaces. We show that when the number of measurements $m$ is large enough with respect to some constant $m_{0}$, which scales as the Gaussian mean width of the union of subspaces, we have that the reconstruction result $\hat{\bm{x}}$ satisfies

[TABLE]

with high probability, where $c_{0}$ is a constant that tends to zero as the number of iterations grows (to be more precise $c_{0}=\mathcal{O}(m^{-t/2})$, where $t$ is the number of iterations), and $c_{1}$ is another constant that obeys $c_{1}=\mathcal{O}(m^{-1})$. The significance of this result, except for its generality, is in showing that the effect of any stationary noise, that does not depend on $\bm{A}$, tends to zero as the number of measurements grows. The recovery guarantees hold for any algorithm that is an exact instance of GCoSaMP. Two such examples are CoSaMP and ADMiRA [27]. Another main contribution of the paper is introducing a new practical method for signal reconstruction, and possible decomposition, in a combined sparse synthesis and cosparse analysis model. The proposed technique is a relaxed version of GCoSaMP, specialized to this combined model. The idea of combining two sparse models in order to describe a signal has already been practiced in [36], where two adapted dictionaries (e.g., one dictionary adapted to represent textures, and the other to represent cartoons) were used. However, our scheme is significantly different than the one in [36]. Our solution covers also the case of compressive and limited measurements, uses the recently proposed cosparse analysis model instead of a second sparse synthesis model, and provides a natural way to combine the two models, and any other ones.

The remainder of the paper is organized as follows. In Section II we present the notation we use in this work and some preliminaries. In Section III we present the GCoSaMP algorithm and prove the main theorem of the paper. We also discuss the significance of the results in showing the denoising capabilities of the algorithm. In Section IV we discuss the implications of our results for several specific models that fall into the general union of subspaces framework, and in particular, we present SACoSaMP, a new signal recovery algorithm for sparse synthesis and cosparse analysis combined model. Experiments that support our theoretical results for GCoSaMP, and experiments that examine the performance of SACoSaMP, are given in Section V. Section VI concludes the paper.

II Notations and Preliminaries

In this section we present notations, definitions and basic results that will be used in our work. We write $\|\cdot\|_{2}$ for the Euclidean norm of a vector, and $\|\cdot\|$ for the spectral norm of a matrix. We denote the unit Euclidean ball and sphere in $\mathbb{R}^{n}$ by $\mathbb{B}^{n}$ and $\mathbb{S}^{n-1}$, respectively, and the $n\times n$ identity matrix by $\bm{I}_{n}$. By abuse of notation, $C,c$ denote positive absolute constants whose values may change from instance to instance. Given a linear subspace $\mathcal{V}\subset\mathbb{R}^{n}$, we denote the orthogonal projection onto $\mathcal{V}$ by $\bm{P}_{\mathcal{V}}$, and the orthogonal projection onto the orthogonal complement of $\mathcal{V}$ by $\bm{Q}_{\mathcal{V}}=\bm{I}_{n}-\bm{P}_{\mathcal{V}}$. The sum of sets $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ is defined as

[TABLE]

We denote by $\mathcal{S}$ a known, possibly infinite, set of finite-dimensional linear subspaces $\{\mathcal{V}_{i}\}$. The $B$-order sum for the set $\mathcal{S}$, with an integer $B\geq 1$, is defined as

[TABLE]

Note that $\mathcal{S}^{B}$ is also a known, possibly infinite, set of linear subspaces, and $\mathcal{S}\equiv\mathcal{S}^{1}$. The notation $\mathcal{U}^{B}$ stands for the union of all subspaces in the set $\mathcal{S}^{B}$, namely

[TABLE]

It is easy to show that $\mathcal{U}^{B}$ can also be expressed as

[TABLE]

We assume that the signal belongs to one of the subspaces in $\mathcal{S}$, i.e., $\bm{x}\in\mathcal{V}_{0}$, and $\mathcal{V}_{0}\in\mathcal{S}$. However, this subspace is unknown, and therefore our signal model is $\bm{x}\in\mathcal{U}$, where $\mathcal{U}\equiv\mathcal{U}^{1}$.

To exemplify the definitions above, we consider the traditional sparsity model $\{\bm{x}\in\mathbb{R}^{n}:\|\bm{x}\|_{0}\leq k\}$, where $\|\cdot\|_{0}$ denotes the $\ell_{0}$ pseudo-norm, that counts the number of nonzero elements of a vector. Clearly, if $\mathcal{S}$ is the set of $\binom{n}{k}$ subspaces, where each subspace is spanned by a different selection of $k$ columns of the identity matrix $\bm{I}_{n}$, then $\bm{x}$ resides in one of the subspaces in $\mathcal{S}$, and the sparsity model can be compactly expressed as $\bm{x}\in\mathcal{U}$. Moreover, in this case, the set $\mathcal{S}^{B}$ contains all the $\binom{n}{Bk}$ subspaces spanned by different selections of $Bk$ columns of $\bm{I}_{n}$, and other subspaces of lower dimensions, which are contained in those $\binom{n}{Bk}$ subspaces. Lastly, from either (5) or (6), we have $\mathcal{U}^{B}=\{\bm{x}\in\mathbb{R}^{n}:\|\bm{x}\|_{0}\leq Bk\}$.

Definition 2.1**.**

The Gaussian mean width of a set $\mathcal{K}\in\mathbb{R}^{n}$ is defined as

[TABLE]

where the expectation is taken over $\bm{g}\sim\mathcal{N}(\mathbf{0},\bm{I}_{n})$, a vector of independent zero-mean unit-variance Gaussians.

The following lemma results from [37, Cor 1.2] together with standard concentration of measure for Gaussian random variable [38].

Lemma 2.2**.**

Let $\mathcal{K}\in\mathbb{R}^{n}$ be a closed set and let $\bm{A}$ be an $m\times n$ matrix with entries independently drawn from the standard normal distribution $\mathcal{N}(0,1)$. Also define

[TABLE]

where $\bm{g}\in\mathbb{\bm{R}}^{m}$ is distributed as $\mathcal{N}(\mathbf{0},\bm{I}_{m})$. Then for all $\bm{u}\in\mathcal{K}$ and $\eta>0$ we have

[TABLE]

holds with probability at least $1-\mathrm{e}^{-\frac{\eta^{2}}{2}}$. Furthermore, for all $\bm{u}\in\mathcal{K}$ and $\eta>0$ we have

[TABLE]

holds with probability at least $1-\mathrm{e}^{-\frac{\eta^{2}}{2}}$.

Lemma 2.3**.**

Consider the assumptions and notations of Lemma 2.2. If the positive parameters $\mu$ and $\eta$ satisfy $\mu\leq\left(b_{m}+w(\mathcal{U}^{B}\cap\mathbb{S}^{n-1})+\eta\right)^{-2}$ then for all $\mathcal{V}\in\mathcal{S}^{B}$, we have

[TABLE]

holds with probability at least $1-2\mathrm{e}^{-\frac{\eta^{2}}{2}}$.

Proof.

Note that if $\bm{z}\in\mathbb{R}^{n}$ has unit Euclidean norm, then $\bm{P}_{\mathcal{V}}\bm{z}\in\mathcal{V}\cap\mathbb{B}^{n}$. Since $\mathcal{V}\cap\mathbb{B}^{n}\subset\mathcal{U}^{B}\cap\mathbb{B}^{n}$, using Lemma 2.2 we have

[TABLE]

and

[TABLE]

hold together with probability at least $1-2\mathrm{e}^{-\frac{\eta^{2}}{2}}$. Assuming that this event occurred, we have

[TABLE]

The first inequality uses (10) and $\mu>0$. The second inequality uses $\left\|\bm{P}_{\mathcal{V}}\bm{z}\right\|_{2}\leq\left\|\bm{z}\right\|_{2}=1$ and the fact that the expression in the brackets is nonnegative, due to (11) and the upper bound on $\mu$. The last equality follows from the fact that $\bm{P}_{\mathcal{V}}(\bm{I}_{n}-\mu\bm{A}^{*}\bm{A})\bm{P}_{\mathcal{V}}$ is Hermitian. ∎

III Algorithm and Recovery Guarantees

The generalized CoSaMP (GCoSaMP) algorithm, which is a generalization of the sparsity-based algorithm proposed in [25], is presented in Algorithm 1. The main difference compared to CoSaMP is that support selections are replaced by subspace selections. For a fixed value of $B$, we define a subspace selection problem as finding a subspace $\hat{\mathcal{V}}\in\mathcal{S}^{B}$, such that the orthogonal projection onto $\hat{\mathcal{V}}$ is the closest (in Euclidean distance) to the signal proxy/estimation $\tilde{\bm{x}}$, i.e.,

[TABLE]

Although this problem may be demanding or intractable in general, for many models there exist efficient methods to obtain exact solutions, and for many others, approximation methods may still produce good results. In Section IV we discuss several specific models.

Our main results are given in the following theorem and corollary.

Theorem 3.1**.**

Let $\bm{y}=\bm{A}\bm{x}+\bm{e}$, where $\bm{x}\in\mathcal{U}$ and $\bm{A}$ is an $m\times n$ matrix with entries independently drawn from the standard normal distribution, then with the notation of Algorithm 1, we have

[TABLE]

holds for all $t$, with probability at least $1-6\mathrm{e}^{-\frac{\eta^{2}}{2}}$, where

[TABLE]

and $\mu_{1}$, $\mu_{2}$ and $\eta$ are positive parameters satisfying

[TABLE]

The iterates converge if $m>14.5^{2}(w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})+\eta)^{2}+1$.

Proof.

We divide the proof into several lemmas.

Lemma 3.2**.**

With the notation of Algorithm 1, we have

[TABLE]

Lemma 3.3**.**

With the notation of Algorithm 1 and Theorem 3.1, we have

[TABLE]

holds for all $t$, with probability at least $1-3\mathrm{e}^{-\frac{\eta^{2}}{2}}$.

Lemma 3.4**.**

With the notation of Algorithm 1 and Theorem 3.1, we have

[TABLE]

holds for all $t$, with probability at least $1-3\mathrm{e}^{-\frac{\eta^{2}}{2}}$.

The proofs of Lemmas 3.2, 3.3 and 3.4 appear in Appendices A, B and C, respectively. Combining the above three lemmas, we immediately get (14).

Note that $\mu_{1}$ and $\mu_{2}$ are two arbitrary parameters (with restricted permitted range) that can be used to optimize the bound in (14). Increasing $\mu_{1}$ leads to lower value of $\rho_{1}(\eta)$ and higher value of $\xi_{1}(\eta)$, and the same goes for $\mu_{2}$. However, in order to ensure that the iterates converge, it is enough to consider only the first term in (14) and require that $\frac{4\rho_{2}(\eta)}{\sqrt{1-\rho_{1}^{2}(\eta)}}<1$, which is equivalent to $16\rho_{2}^{2}(\eta)+\rho_{1}^{2}(\eta)<1$.

To simplify the notations we define $m_{0}\triangleq(w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})+\eta)^{2}$. Following the calculation from [10], we have

[TABLE]

where in the last two inequalities we used $m\geq 14.5^{2}m_{0}+1$. Setting both $\mu_{1}$ and $\mu_{2}$ at their highest permitted value $\left(b_{m}+\sqrt{m_{0}}\right)^{-2}$, we have

[TABLE]

where the last inequality follows from (22). ∎

This leads us to the following corollary.

Corollary 3.5**.**

Assuming the conditions of Theorem 3.1, define $m_{0}\triangleq(w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})+\eta)^{2}$ and assume $m\geq 14.5^{2}m_{0}+1$, then with the notation of Algorithm 1, we have

[TABLE]

holds with probability at least $1-6\mathrm{e}^{-\frac{\eta^{2}}{2}}$, where

[TABLE]

Proof.

From Theorem 3.1, using recursion we get

[TABLE]

where

[TABLE]

and for brevity we omitted the dependencies on $\eta$. Similarly to the second part of the proof of Theorem 3.1, we set both $\mu_{1}$ and $\mu_{2}$ at their largest permitted value $\left(b_{m}+\sqrt{m_{0}}\right)^{-2}$. We turn to bounding the basic terms in (28) and (29). We repeatedly use the fact $\sqrt{m}\geq b_{m}\geq\frac{m}{\sqrt{m+1}}>\sqrt{m_{0}}$. First, note that

[TABLE]

Therefore,

[TABLE]

and we get

[TABLE]

Using (27) and the last two inequalities we get (24).

∎

For $m\gg m_{0}$, it is easy to see that $\rho(m)$ and $\xi(m)$ behave like $16\sqrt{\frac{m_{0}}{m}}$ and $6\frac{\sqrt{m_{0}}}{m}$, respectively. Therefore, the noise coefficient in (24) also behaves like $6\frac{\sqrt{m_{0}}}{m}$. Assuming that the noise characteristics are the same at all measurements, it follows that $\|\bm{e}\|_{2}=\mathcal{O}(\sqrt{m})$. Thus, we conclude that the effect of the noise tends to zero as the number of measurements grows. However, please note that since our derivation assumes that $\bm{e}$ does not depend on $\bm{A}$ (see the remarks before (67) and (75)), the bound does not cover the case of adversarial noise $\bm{e}=-c\bm{A}\bm{x}$.

Discussion. Let us compare our bound with previous results. As was mentioned in Section II, by choosing $\mathcal{S}$ as the set of $\binom{n}{k}$ subspaces, where each subspace is spanned by a different selection of $k$ columns of the identity matrix $\bm{I}_{n}$, we get $\mathcal{U}=\{\bm{x}\in\mathbb{R}^{n}:\|\bm{x}\|_{0}\leq k\}$. That is to say, $\mathcal{U}\equiv\mathcal{U}_{k\text{-sparse}}$ is the set of $k$-sparse vectors in $\mathbb{R}^{n}$. Also, the set $\mathcal{S}^{B}$ contains all the $\binom{n}{Bk}$ subspaces spanned by different selections of $Bk$ columns of $\bm{I}_{n}$ (and other subspaces of lower dimensions, which are contained in those $\binom{n}{Bk}$ subspaces). Therefore, in this case GCoSaMP reduces to the original sparsity-based CoSaMP [25]. In order to use the theoretical results above, the evaluation of the quantity $w(\mathcal{U}^{4}_{k\text{-sparse}}\cap\mathbb{S}^{n-1})$ is required. Conveniently, since $\mathcal{U}^{B}_{k\text{-sparse}}=\{\bm{x}\in\mathbb{R}^{n}:\|\bm{x}\|_{0}\leq Bk\}$ we can use the following result from [11]

[TABLE]

where $c$ and $C$ are two positive constants. Before we present the reconstruction error bound that was obtained in [25], recall that we say that $\tilde{\bm{A}}$ satisfies the RIP of order $k$ with constant $\delta_{k}\in[0,1)$ if

[TABLE]

holds for all $\bm{x}$ satisfying $\|\bm{x}\|_{0}\leq k$. In [25, Thm 4.1] it was shown that given the measurements $\bm{y}=\tilde{\bm{A}}\bm{x}+\tilde{\bm{e}}$, where $\bm{x}$ is a $k$-sparse vector and $\tilde{\bm{A}}$ satisfies the RIP of order $4k$ with $\delta_{4k}\leq 0.1$, the solution of the $t$-th iteration $\bm{x}^{t}$ satisfies

[TABLE]

where

[TABLE]

Note that we present the bound without substituting $\delta_{2k}\leq\delta_{3k}\leq\delta_{4k}\leq 0.1$, as done in [25]. Note also that the columns of $\tilde{\bm{A}}$ are assumed to be normalized in order to have $\delta_{1}=0$, while we did not normalize the columns of $\bm{A}$ to have unit variance (by dividing it by $\sqrt{m}$). Therefore, in order to have the same averaged noise level as in our model, we have $\|\tilde{\bm{e}}\|_{2}=\|\bm{e}\|_{2}/\sqrt{m}$, and if we also assume that the noise characteristics are the same at all measurements, it follows that $\|\tilde{\bm{e}}\|_{2}=\mathcal{O}(1)$. Now, it is easy to see that for fixed $k,n$, even if $m$ grows without bound such that the RIP constants tend to zero [32], the effect of the noise does not vanish in (39), contrary to our result in (24). The last remark is not surprising, since (39) covers the case of adversarial noise. However, it emphasizes that our bound (24) is more useful if the denoising capabilities of the algorithm are of interest (and not only its stability).

We turn to comparing (24) with the bound obtained in [39] for CoSaMP under the assumption of Gaussian noise. In [39, Thm 2.3] it was shown that given the measurements $\bm{y}=\tilde{\bm{A}}\bm{x}+\tilde{\bm{e}}$, where $\bm{x}$ is a $k$-sparse vector, $\tilde{\bm{A}}$ satisfies the RIP of order $4k$ with $\delta_{4k}\leq 0.1$ and $\tilde{\bm{e}}\sim\mathcal{N}(\mathbf{0},\tilde{\sigma}^{2}\bm{I}_{m})$, after $t\geq t^{*}$ iterations we have that

[TABLE]

holds with probability exceeding $1-(\sqrt{\pi(1+a)\mathrm{log}n}\cdot n^{a})^{-1}$. The constant $\tilde{C}$ and the value of $t^{*}$ are given in [39]. For simplicity, we examine our bound in (24) only for $m\gg m_{0}$. In this case, as we mentioned above, the noise coefficient behaves like $6\frac{\sqrt{m_{0}}}{m}$ and the first term can be as small as we desire for large enough $t$. Therefore, using (37) and the definition of $m_{0}$, our bound in (24) leads to

[TABLE]

with probability exceeding $1-6\mathrm{e}^{-\frac{\eta^{2}}{2}}$. Once again, since we did not normalize the columns of $\bm{A}$, and assuming that the noise is stationary, we have that $\frac{\|\bm{e}\|_{2}}{m}$ is equivalent to $\tilde{\sigma}$. Both bounds have quite similar dependencies on $k,n$, and they both agree that CoSaMP can recover signals with an effective reduction of the additive noise. However, the bound in [39] is more general in terms of the distribution of the measurement matrix, and the bound obtained here is more general in terms of the noise, which can be drawn from any distribution, or even exhibit nonstationarity. The only restriction on the noise is that it does not depend on $\bm{A}$. As a supporting example, in Section V-A we conduct experiments with noise drawn from Laplace distribution. We also stress that the bound obtained here can be applied to signal models beyond sparsity. Lastly, we note that our assumptions on the measurement matrix and noise vector are similar to those in [40], where bounds on the reconstruction error of projected gradient descent algorithm were developed.

IV Examples

The results in Section III require the evaluation of the quantity $w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})$, which depends on the subspaces included in the set $\mathcal{S}$ (recall the definitions in (4) and (5)). However, note that having an upper bound on $w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})$ is enough for bounding the error in (24). This is true since plugging an upper bound on $m_{0}$ in (25) and (26) yields upper bounds on $\rho(m)$ and $\xi(m)$, that can be used to upper bound the error in (24). The following lemma may facilitate the evaluation of a bound on the Gaussian mean width in some cases.

Lemma 4.1**.**

Let the union of subspaces $\mathcal{U}$ be the sum of $N$ unions of subspaces $\{\mathcal{U}_{i}\}_{i=1}^{N}$, i.e., $\mathcal{U}=\sum\limits_{i=1}^{N}\mathcal{U}_{i}$. Then, we have

[TABLE]

Proof.

Let $\{\mathcal{K}_{i}\}_{i=1}^{N}$ be $N$ arbitrary sets, we have the following fact

[TABLE]

Next, note that $\mathcal{U}\cap\mathbb{B}^{n}\subset\sum\limits_{i=1}^{N}(\mathcal{U}_{i}\cap\mathbb{B}^{n})$, and therefore we have

[TABLE]

Finally, using (45) together with (46), and exploiting the property $w(\mathcal{U}_{i}\cap\mathbb{B}^{n})=w(\mathcal{U}_{i}\cap\mathbb{S}^{n-1})$ (see the explanation below (75) in Appendix C), we get (44). ∎

A direct result of (6) and Lemma 4.1 is

[TABLE]

However, in many specific models, $\mathcal{U}^{B}$ is very similar to $\mathcal{U}^{1}$, and a bound on $w(\mathcal{U}^{B}\cap\mathbb{S}^{n-1})$ which is evaluated directly, can be tighter than (47). Nevertheless, Lemma 4.1 is beneficial when a combined model is being considered, as discussed below.

The computational complexity of GCoSaMP is dominated by the complexity of the subspace selections problems. For certain choices of $\mathcal{S}$, optimal subspace selection methods exist and can be implemented efficiently, while for other choices of $\mathcal{S}$, relaxations of GCoSaMP are required in order to obtain practical algorithms. These relaxed algorithms may not possess the recovery guarantees from Section III. However, they are still closely related to GCoSaMP.

In the sequel, we discuss various specific models that fall into the general union of subspaces framework, and mainly focus on the following two subjects: Evaluation of the Gaussian mean width and relaxations of GCoSaMP.

IV-A Sparsity and Structured Sparsity Models

As discussed in Section III, choosing $\mathcal{S}$ as the set of $\binom{n}{k}$ subspaces, where each subspace is spanned by a different selection of $k$ columns of $\bm{I}_{n}$, we have that $\mathcal{U}^{B}\equiv\mathcal{U}^{B}_{k\text{-sparse}}$ is the set of $Bk$-sparse vectors in $\mathbb{R}^{n}$. The quantity $w(\mathcal{U}^{B}_{k\text{-sparse}}\cap\mathbb{S}^{n-1})$ is bounded in (37). Note that the upper bound in (37) is tighter than the one obtained by (47). For this choice of $\mathcal{S}$, GCoSaMP reduces to CoSaMP [25]: The subspace selection problems, i.e., determining $\tilde{\mathcal{V}}^{t}$ and $\mathcal{V}^{t}$, involve finding the best $2k$-sparse or $k$-sparse approximations to a given vector, and thus reduce to support selections that can be calculated easily using thresholding.

Another instance of union of subspaces is structured sparsity, in which the signal model is restricted to predefined sparsity patterns. Thus, assuming that the total sparsity is $k$, the set $\mathcal{S}$ includes only part of the $\binom{n}{k}$ subspaces associated with the unconstrained $k$-sparse vectors. For certain sparsity patterns, such as block-sparsity or tree-sparsity, there exist efficient methods to obtain a solution to the subspace selection problem [26]. In order to evaluate the Gaussian mean width for these patterns, we make use of the following lemma, which is an extension of Lemma 2.3 in [11].

Lemma 4.2**.**

Let the set $\mathcal{S}$ be a subset of the $\binom{n}{k}$ subspaces associated with the unconstrained $k$-sparse vectors, such that $|\mathcal{S}|\leq\operatorname{exp}(\gamma(k,n))$, where $\gamma(k,n)$ is some function of $k$ and $n$. Then, we have

[TABLE]

Proof.

Following the proof in [11] (which is given in more details in [41]), but replacing $|\mathcal{S}|=\binom{n}{k}\leq\left(\frac{\mathrm{e}n}{k}\right)^{k}=\operatorname{exp}(k+k\operatorname{log}\frac{n}{k})$ with $|\mathcal{S}|\leq\operatorname{exp}(\gamma(k,n))$, it is straightforward to obtain (48). ∎

For tree-sparsity we have $|\mathcal{S}|\leq\frac{(2\mathrm{e})^{k}}{k+1}=\operatorname{exp}(k+k\operatorname{log}\frac{2}{(k+1)^{1/k}})$ [26], i.e., $\gamma(k,n)=k+k\operatorname{log}\frac{2}{(k+1)^{1/k}}$, which yields

[TABLE]

For block-sparsity, where the blocks are disjoint and of size $J$, we have $|\mathcal{S}|=\binom{n/J}{k/J}\leq\left(\frac{\mathrm{e}n}{k}\right)^{k/J}=\operatorname{exp}(k/J+k/J\operatorname{log}(n/k))$, i.e., $\gamma(k,n)=k/J+k/J\operatorname{log}(n/k)$, which yields

[TABLE]

We remind the reader that we may calculate $w(\mathcal{U}^{B}\cap\mathbb{S}^{n-1})$ from $w(\mathcal{U}^{1}\cap\mathbb{S}^{n-1})$ using (47).

IV-B Low-Rank Model

Note that using (column or row stack) vectorization, our theory can be applied immediately to matrices. Assuming that the signal is a rank-$r$ matrix $\bm{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ (now the signal dimension is $n=n_{1}\cdot n_{2}$), using singular value decomposition (SVD), $\bm{X}$ can be written as

[TABLE]

where $\bm{u}_{i},\bm{v}_{i}$ are the columns of the orthonormal matrices $\bm{U}$ and $\bm{V}$, respectively, associated with nonzero singular value $\sigma_{i}$. It is easy to see that all the rank-$r$ matrices that can be written as a linear combination of the same $r$ rank-1 matrices $\{\bm{u}_{i}\bm{v}_{i}^{*}\}_{i=1}^{r}$ lie on the subspace $\left\{\sum\limits_{i=1}^{r}a_{i}\bm{u}_{i}\bm{v}_{i}^{*}:a_{i}\in\mathbb{R}\right\}$. Therefore, we can choose $\mathcal{S}$ as the infinite set of such subspaces, each spanned by $r$ rank-1 matrices, and the set of all rank-$r$ matrices can be seen as an infinite union of these subspaces. As a result from this choice of $\mathcal{S}$, the set $\mathcal{S}^{B}$ is also an infinite set of subspaces, where each subspace is spanned by no more than $Br$ rank-1 matrices, and GCoSaMP reduces to ADMiRA [27]. In this case, the subspace selection problems involve finding the best rank-$2r$ or rank-$r$ approximations to a given matrix, which can be done efficiently by truncating the SVD of the matrix. We turn to evaluate the Gaussian mean width for $\mathcal{U}^{B}\equiv\mathcal{U}^{B}_{\text{rank-}r}$. Note that we have the following identity for the Frobenius norm $\|\bm{X}\|_{F}$

[TABLE]

where $\mathrm{diag}$ denotes the main diagonal of a matrix. Therefore, applying Cauchy-Schwarz inequality to the singular values of $\bm{X}$, we have that the set $\mathcal{U}^{B}_{\text{rank-}r}\cap\mathbb{S}^{n-1}$ is contained in the convex set

[TABLE]

where $\|\bm{X}\|_{*}$ denotes the nuclear norm, i.e., the sum of the singular values of $\bm{X}$. Therefore, the upper bound on $w(\mathcal{K}_{n_{1},n_{2},Br})$ in [11] also upper bounds the Gaussian mean width of the set we are interested in

[TABLE]

IV-C Sparse Synthesis Model

Other choices of $\mathcal{S}$ may require relaxation of GCoSaMP in order to reduce its computational complexity. For example, [28, 29] considered the synthesis model, which assumes that the signal $\bm{x}$ admits a $k$-sparse representation $\bm{\alpha}$ in a given dictionary $\bm{D}\in\mathbb{R}^{n\times d}$, i.e., $\bm{x}=\bm{D}\bm{\alpha}$, where $\bm{\alpha}$ is a $k$-sparse vector in $\mathbb{R}^{d}$ (typically $d>n$). Denote the columns of $\bm{D}$ by $\{\bm{d}_{i}\}_{i=1}^{d}$. For this model, $\mathcal{S}$ consists of $\binom{d}{k}$ subspaces, where each subspace is spanned by a different selection of $k$ vectors from $\{\bm{d}_{i}\}_{i=1}^{d}$. The set $\mathcal{S}^{B}$ consists of subspaces spanned by no more than $Bk$ vectors from $\{\bm{d}_{i}\}_{i=1}^{d}$. In this case, the subspace selection problem is NP-hard in general [42]. For this reason, both works suggested a signal space CoSaMP (SSCoSaMP) algorithm, in which the optimal subspace selections are replaced with approximated methods. While the theoretical results in [28, 29] require that these approximations satisfy strict near-optimality constraints, both works also empirically examine the performance of SSCoSaMP using practical (but not theoretically backed) selections methods such as hard thresholding or OMP. Let us upper bound $w(\mathcal{U}^{B}_{k\text{-synthesis}}\cap\mathbb{S}^{n-1})$ for a dictionary $\bm{D}$ that has the RIP of order $Bk$ with an RIP-constant $\delta_{Bk}$. It is easy to show that if $\bm{\alpha}$ satisfies $\|\bm{\alpha}\|_{0}=Bk$ and $\|\bm{D}\bm{\alpha}\|_{2}=1$, it also satisfies $\|\bm{\alpha}\|_{1}\leq\sqrt{\frac{Bk}{1-\delta_{Bk}}}$. Now, define the following convex set (for fixed $\bm{\alpha}$ and $\bm{D}$)

[TABLE]

where $\mathrm{conv}$ denotes the convex hull of a set. It was shown in [43] that

[TABLE]

where $C$ is a positive constant. Therefore, since $\bm{x}=\sum\limits_{i=1}^{d}\alpha_{i}\bm{d}_{i}$ leads to $\bm{x}\in\mathcal{K}_{\bm{\alpha},\bm{D}}$, from the arguments above, we have

[TABLE]

IV-D Cosparse Analysis Model

Another case in which the subspace selection problems are NP-hard in general is the cosparse analysis model [2, 44]. The cosparse analysis model assumes that the analyzed vector $\bm{\Omega}\bm{x}$ is sparse, where $\bm{\Omega}\in\mathbb{R}^{p\times n}$ is a possibly redundant analysis operator ($p\geq n$). However, it focusses on the cosupport of $\bm{\Omega}\bm{x}$, i.e., the zeros of $\bm{\Omega}\bm{x}$, instead of on its support. Denote the rows of $\bm{\Omega}$ by $\{\bm{\omega}_{i}^{T}\}_{i=1}^{p}$, the cosupport of $\bm{\Omega}\bm{x}$ by $\Lambda$, and the size of the cosupport by $\ell$. We have that $\bm{x}\in\mathcal{V}_{\Lambda}$, where

[TABLE]

However, since it is assumed that only the size of the cosupport is known, and not the cosupport itself, the set $\mathcal{S}$ consists of all the $\binom{p}{\ell}$ possible subspaces $\mathcal{S}=\{\mathcal{V}_{\Lambda}\,:\,|\Lambda|=\ell\}$, and $\bm{x}$ resides in a union of these subspaces. An analysis version of CoSaMP (ACoSaMP) was proposed in [30]. Even though for some analysis operators, such as the 1D finite difference operator or a concatenation of the identity operator and the 1D finite difference operator, the optimal cosupport selection (or equivalently, subspace selection) can be accomplished in polynomial complexity, in general, no efficient method for optimal cosupport selection is known. Therefore, ACoSaMP replaces optimal cosupport selections with approximated methods. However, except this relaxation, ACoSaMP can be seen as an instance of GCoSaMP. In each iteration, it extracts $\Lambda_{\Delta}^{t}$, a near-optimal cosupport of size $2\ell-p$ of the proxy $\bm{A}^{*}\bm{r}$. This is equivalent to finding a near-optimal subspace $\mathcal{V}_{\Lambda_{\Delta}^{t}}\in\mathcal{S}^{2}$, since the minimal size of the cosupport of a sum of two $\ell$-cosparse vectors is $2\ell-p$ (see [30] for more details). Next, the cosupport for the least squares (LS) estimation step is obtained using intersection of the newly extracted cosupport and the cosupport of the previous iteration, i.e., $\tilde{\Lambda}^{t}=\Lambda^{t-1}\cap\Lambda_{\Delta}^{t}$. This is equivalent to constraining the LS estimation to $\mathcal{V}_{\Lambda^{t}}+\mathcal{V}_{\Lambda_{\Delta}^{t}}$. In practice, efficient methods to obtain near-optimal cosupports are also unknown (for most analysis operators). Therefore, [30] empirically examined the performance of ACoSaMP using simple thresholding. Nevertheless, the algorithm exhibited good performance. Assuming that $\bm{\Omega}$ has linearly independent columns and denoting its pseudoinverse by $\bm{\Omega}^{\dagger}$, the quantity $w(\mathcal{U}^{B}_{\ell\text{-analysis}}\cap\mathbb{S}^{n-1})$ can be evaluated by exploiting results for the sparse synthesis model with an $n\times p$ dictionary $\bm{D}=\bm{\Omega}^{\dagger}$, since $\bm{x}=\bm{\Omega}^{\dagger}\bm{\Omega}\bm{x}$ and $\bm{\Omega}\bm{x}$ is a $(p-\ell)$-sparse vector. Therefore, if $\bm{\Omega}^{\dagger}$ has the RIP of order $B(p-\ell)$ with an RIP-constant $\delta_{B(p-\ell)}$, we have

[TABLE]

for the same reasons that yield (57).

IV-E Sparsity and Low-Rank Combined Model

In the following, we discuss other choices of $\mathcal{S}$ for which the optimal subspace selections are unknown in general. Specifically, we consider the case of a combined model. One such example is recovering a matrix $\bm{X}=\bm{L}+\bm{S}$, that is the sum of a rank-$r$ matrix, $\bm{L}$, and a $k$-sparse matrix, $\bm{S}$. In this case, each subspace in the set $\mathcal{S}$ is the sum of two ”componential subspaces”, one associated with $k$-sparse matrices and the other associated with rank-$r$ matrices (as discussed above). As a result, each subspace in the set $\mathcal{S}^{B}$ is a sum of two componential subspaces as well, one associated with matrices whose sparsity does not exceed $Bk$, and the other associated with matrices whose rank does not exceed $Br$. In [45], an algorithm named SpaRCS was proposed for this problem. This algorithm resembles the general algorithm that we present here, except that in SpaRCS the subspace selection steps and the LS estimation step are performed separately for each componential subspace. Although lacking theoretical performance guarantees, SpaRCS has shown good performance in experiments. Splitting the subspace selection allows performing hard thresholding for the sparse part and truncated SVD for the low-rank part, that can be seen as an approximation to the unknown optimal subspace selection. Denoting the union of subspaces for this model by $\mathcal{U}^{B}_{k\text{-sparse, rank-}r}$, we have $\mathcal{U}^{B}_{k\text{-sparse, rank-}r}=\mathcal{U}^{B}_{k\text{-sparse}}+\mathcal{U}^{B}_{\text{rank-}r}$. Therefore, using Lemma 4.1 we get

[TABLE]

where the terms in the right-hand side of the last inequality are bounded in (37) and (54). However, note that in order to determine that an algorithm has the guarantees from Section III, it must be an exact instance GCoSaMP. Such an algorithm is still unknown in this case, as optimal subspace selection methods are unknown.

IV-F Sparse Synthesis and Cosparse Analysis Combined Model

Another interesting combined model is the one that uses both the sparse synthesis model and the cosparse analysis model, and considers a signal $\bm{x}=\bm{x}_{1}+\bm{x}_{2}$, where $\bm{x}_{1}$ admits a $k$-sparse representation $\bm{\alpha}$ in a given dictionary $\bm{D}\in\mathbb{R}^{n\times d}$, and $\bm{x}_{2}$ is $\ell$-cosparse under a given analysis operator $\bm{\Omega}\in\mathbb{R}^{p\times n}$, i.e., $\bm{\Omega}\bm{x}_{2}$ has at least $\ell$ zeros. This combined model may be advantageous for images that contain both texture and piecewise constant (cartoon) parts. The texture part may admit a sparse representation in a suitable dictionary, and the cartoon part is expected to have a cosparse analysis representation under the finite difference operator, even if the original image does not comply with each of these models separately. Clearly, this is the case when the texture and cartoon parts overlap.

Denoting the union of subspaces for this model by $\mathcal{U}^{B}_{k\text{-synthesis, }\ell\text{-analysis}}$, we have $\mathcal{U}^{B}_{k\text{-synthesis, }\ell\text{-analysis}}=\mathcal{U}^{B}_{k\text{-synthesis}}+\mathcal{U}^{B}_{\ell\text{-analysis}}$. Therefore, using Lemma 4.1 we get

[TABLE]

where the terms in the right-hand side of the last inequality are bounded in (57) and (59) for $\bm{D}$ and $\bm{\Omega}^{\dagger}$ that have the RIP. However, as follows from our discussion on each sub-model, the subspace selection problem for the combined model is NP-hard in general. Therefore, we propose a practical reconstruction algorithm for this model, inspired by GCoSaMP.

We propose synthesis-analysis CoSaMP (SACoSaMP) for reconstruction of $\bm{x}$ from $m$ ($m<n$) compressive linear measurements given by $\bm{y}=\bm{A}\bm{x}+\bm{e}$. The algorithm is presented in Algorithm 2. It can also be used to directly recover $\bm{x}_{1}$ and $\bm{x}_{2}$ separately. It combines SSCoSaMP [28] and ACoSaMP [30] in a way that resembles the way SpaRCS combines CoSaMP [25] and ADMiRA [27], except that the LS estimation step in SACoSaMP is unified, a fact that can be considered as being more loyal to GCoSaMP, and is justified by our experiments in Section V-B. Let us clarify the notations being used in Algorithm 2. We use $[1..p]$ to denote the integers between $1$ and $p$, $\operatorname{supp}(\bm{u},k)$ for the indices of the $k$ largest elements (in magnitude) of the vector $\bm{u}$, and $\operatorname{cosupp}(\bm{u},\ell)$ for the indices of the $\ell$ smallest elements (in magnitude) of the vector $\bm{u}$. The notation $\bm{\alpha}|_{T}$ stands for a sub-vector of $\bm{\alpha}$ with elements that belong to the set $T$, and $T^{C}$ denotes the complementary set of $T$. $\bm{D}_{T}$ is a sub-matrix of $\bm{D}$ with columns corresponding to the set $T$, and by abuse of notation, $\bm{\Omega}_{\Lambda}$ is a sub-matrix of $\bm{\Omega}$ with rows that belong to the set $\Lambda$. Finally, $\bm{Q}_{\Lambda}=\bm{I}_{n}-\bm{\Omega}_{\Lambda}^{\dagger}\bm{\Omega}_{\Lambda}$ is the orthogonal projection onto the orthogonal complement of the row space of $\bm{\Omega}_{\Lambda}$. We note that the extracted cosupport $\Lambda_{\Delta}^{t}$ is of size $\ell$ rather than $2\ell-p$, since it was empirically shown in [30] that this modification improves the performance of ACoSaMP.

SACoSaMP is a relaxed version of GCoSaMP. The support and cosupport selections are separated, and use simple thresholding, which is not optimal in general, and thus the resulted subspace may not be optimal. Therefore, even if one can evaluate the Gaussian mean width for the combined model described above, SACoSaMP does not possess the recovery guarantees from Section III. However, it is still inspired by the general scheme of GCoSaMP and by the idea of not being restricted to a single traditional model. The usefulness of the algorithm is demonstrated in Section V-B.

V Experiments

The experiments section is divided into two parts. In the first part, we conduct an experiment that supports the theoretical result on the denoising capabilities of GCoSaMP. Specifically, we focus on the original sparsity-based CoSaMP, which is equivalent to GCoSaMP when $\mathcal{S}$ is the set of subspaces associated with the traditional sparsity model. In the second part, we consider the sparse synthesis and cosparse analysis combined model. We examine the performance of SACoSaMP, focusing on its application to simultaneous image reconstruction and structured noise removal.

V-A The vanishing effect of stationary noise in CoSaMP

As was mentioned above, when $\mathcal{S}$ is the set of subspaces associated with the traditional sparsity model, GCoSaMP reduces to CoSaMP. Using synthetic data, we demonstrate that for fixed $k,n$ and large enough $m$, the reconstruction error (after the convergence of the algorithm) decreases proportionally to $1/\sqrt{m}$. This empirical result supports the theoretical result in Section III, that for $\|\bm{e}\|_{2}=\mathcal{O}(\sqrt{m})$ (stationary noise) and $m\gg m_{0}$, the reconstruction error behaves like $\mathcal{O}(\sqrt{\frac{m_{0}}{m}})$. We use a signal of length $n=2\mathrm{e}5$ with $k=5$ nonzero entries independently drawn from the standard normal distribution, and zero-mean noise drawn from Laplace distribution. The measurement matrix $\bm{A}$ is an $m\times n$ matrix with $m$ varying between $100$ and $1\mathrm{e}5$. For each value of $m$ we perform 30 Monte Carlo experiments, in which the entries of $\bm{A}$ are independently drawn from the standard normal distribution and the noise is scaled such that $\|\bm{e}\|_{2}/\|\bm{A}\bm{x}\|_{2}=0.01$. Fig. 1. shows the Euclidean norm of the reconstruction error. The markers represent empirical results and the solid line is computed according to $cm^{-0.5}$, where the constant $c$ compensates for the shift at the last point ($m=1\mathrm{e}5$). Starting from $m=1\mathrm{e}3$, the results are in good agreement, even though our theoretical convergence guarantees required larger values of $m$, since $14.5^{2}m_{0}$ is about $1.18\mathrm{e}4$ (see (37)).

V-B Reconstruction and structured noise removal using SACoSaMP

We turn now to examine the performance of SACoSaMP. We focus on its application to simultaneous image reconstruction and denoising, when the noise is nonstationary and structured. Such noise may result from intermittent interference between electronic components during the image acquisition process. In our experiments we use the sparse synthesis model for the noise. The dictionary $\bm{D}$ is the (inverse) local discrete cosine transform (DCT) with window size of $64\times 64$ pixels and overlap of 32 pixels. The $16\times 16$ lowest frequencies in each window are not included in the dictionary. The texture image $\bm{x}_{1}$ that represents the noise is obtained using $k=500$ random dictionary words with Gaussian coefficients. The original clean image we consider as $\bm{x}_{2}$ is the $512\times 512$ natural image house given in Fig. 2(a). The Frobenius norm of the noise image is scaled to 0.1 of the Frobenius norm of the house image. The noisy image $\bm{x}_{1}+\bm{x}_{2}$, quantized to 8 bits per pixel (bpp), is given in Fig. 2(b). The sampling operator $\bm{A}$ is a two dimensional Fourier transform that measures only $m/n=0.254$ of the points in the Fourier domain according to the binary mask given in Fig. 2(c). The cosparse operator $\bm{\Omega}$ is the finite difference analysis operator that computes horizontal and vertical discrete derivatives of an image. The na**\̈bm{i}**ve reconstruction obtained by inverse Fourier transform on the measurements is shown in Fig. 2(d). Its peak signal to noise ratio (PSNR) with respect to the noisy image is 24.11 dB. Clearly, it is also not a satisfying reconstruction of $\bm{x}_{2}$.

Figs. 2(e) and 2(f) show the results of SACoSaMP obtained after 6 iterations, using $k=500$ and $\ell=386464$ ($\ell$ is slightly larger than the number of entries with magnitude lower than 3 in the analysis representation of $\bm{x}_{2}$). Fig. 2(e) shows the reconstruction of $\bm{x}_{2}$, quantized to 8 bpp. Its PSNR, with respect to the clean house image, is 35.16 dB. One can see that this cosparse part still includes significant portion of the texture of the original house image, but none of the texture of the noise. Fig. 2(f) shows the (full) reconstruction of $\bm{x}_{1}+\bm{x}_{2}$ quantized to 8 bpp. Its PSNR, with respect to the noisy image, is 35.28 dB.

Regarding the implementation of SACoSaMP, in each iteration we take only the real part of the results and we use conjugate gradients for the minimization process. We stop the algorithm when there is no significant reduction in the residual error $\|\bm{r}\|_{2}$.

Using the same parameters but modifying SACoSaMP by splitting the LS step into two separated LS estimations, one for each sub-model (as done in SpaRCS [45]), we get inferior results. The number of iterations until convergence is 14, and the reconstructions have smaller PSNR. Fig. 2(g) shows the reconstruction of $\bm{x}_{2}$. Its PSNR (with respect to the clean image) is 29.34 dB. The reconstruction of $\bm{x}_{1}+\bm{x}_{2}$ has PSNR of 28.78 dB (with respect to the noisy image). We do not show the latter in the paper due to space limitations. We try other several values of $\ell$, but no significant improvement is obtained.

We also examine the performance of ACoSaMP (which is based only on the cosparse model). We try to reconstruct the clean house image using large value of $\ell$, and to reconstruct the noisy image using small value of $\ell$. In both cases the results are not satisfying. For example, the result of ACoSaMP with $\ell=386464$ is shown in Fig. 2(h).

We repeat the experiments using a $256\times 256$ Shepp-Logan phantom image with modified intensities instead of the house image. The binary mask consists of 25 radial lines that contains $m/n=0.095$ of the points in the Fourier domain. The texture image that represents the noise is created as before. The Frobenius norm of the noise image is scaled to 0.2 of the Frobenius norm of the cartoon image. SACoSaMP uses $k=500$ and $\ell=127064$ ($\ell$ is slightly lower than the number of zeros in the analysis representation of the cartoon image). The PSNR of the reconstructed clean image is 35.07 dB, and the PSNR of the reconstructed noisy image is 23.94 dB. The reason for the relatively low PSNR of the latter is the small number of high frequency points in the binary mask. Regarding the algorithm with split LS step (modified SACoSaMP), the best results are received with $\ell=126514$. The PSNR of the reconstructed clean image is 30.71 dB, and the PSNR of the reconstructed noisy image is 23.29 dB. The number of iterations is 17 comparing to 8 when the LS step is unified. As before, the results of ACoSaMP are not satisfying (we show only its result for $\ell=126514$). The inputs and the various results are shown in Figs. 3(a) - 3(h).

VI Conclusions

Under the assumption of Gaussian measurement matrix, we studied a generalized version of the CoSaMP algorithm (GCoSaMP), and provided general recovery guarantees that depend on the Gaussian mean width of the union of subspaces in which the signal resides. We discussed the evaluation of the Gaussian mean width and necessary relaxations of the generalized algorithm in various specific models. We proved that when no relaxation is done, the effect of any (non-adversarial) stationary measurement noise vanishes as the number of measurements grows. This is the case in CoSaMP and ADMiRA, though the proof techniques used in previous works did not lead to this general conclusion. We proposed a new method, inspired by GCoSaMP, for signal recovery in a combined sparse synthesis and cosparse analysis model, and performed experiments that demonstrate its usefulness.

One direction for future research is evaluating the Gaussian mean width of unions of subspaces in specific models for which this quantity is still unknown. Another direction is designing recovery algorithms for combined models, which follows GCoSaMP more closely than our newly proposed SACoSaMP or other solutions such as SpaRCS. Future research may also focus on extending the recovery guarantees of GCoSaMP to measurement matrix with non-Gaussian distribution. This may be done using the proof technique proposed in [46].

Appendix A Proof of Lemma 3.2

Using the triangle inequality, we have

[TABLE]

Note that $\bm{x}\in\mathcal{V}_{0}\in\mathcal{S}$, and also $\bm{x}^{t}\in\mathcal{V}^{t}\in\mathcal{S}$. By construction $\bm{x}^{t}$ is the nearest to $\tilde{\bm{x}}^{t}$, in the $\ell_{2}$-norm sense, among all vectors that reside in one of the subspaces in the set $\mathcal{S}$. Therefore, $\|\tilde{\bm{x}}^{t}-\bm{x}^{t}\|_{2}\leq\|\tilde{\bm{x}}^{t}-\bm{x}\|_{2}$, which together with (62) leads to (19).

Appendix B Proof of Lemma 3.3

To simplify the notations we define $\bm{z}\triangleq\tilde{\bm{x}}^{t}-\bm{x}$. Since $\bm{A}\tilde{\bm{x}}^{t}$ is the orthogonal projection of $\bm{y}$ onto the subspace defined by $\{\bm{A}\bm{u}:\bm{u}\in\tilde{\mathcal{V}}^{t}\}$, then for any $\bm{u}\in\tilde{\mathcal{V}}^{t}$ we have

[TABLE]

Substituting (1) with simple arithmetics gives

[TABLE]

Letting $\mu_{1}$ and $\eta$ be positive parameters satisfying $\mu_{1}\leq\left(b_{m}+w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})+\eta\right)^{-2}$ and turning to look at $\left\|\bm{P}_{\tilde{\mathcal{V}}^{t}}\bm{z}\right\|_{2}^{2}$, we get

[TABLE]

The third equality follows from (64) with $\bm{u}=\bm{P}_{\tilde{\mathcal{V}}^{t}}\bm{z}$, the fourth equality follows from $\bm{P}_{\tilde{\mathcal{V}}^{t}+\mathcal{V}_{0}}\bm{z}=\bm{z}$ and $\bm{P}_{\tilde{\mathcal{V}}^{t}}=\bm{P}_{\tilde{\mathcal{V}}^{t}+\mathcal{V}_{0}}\bm{P}_{\tilde{\mathcal{V}}^{t}}$, and the inequality follows from Cauchy-Schwarz inequality. Define $\tilde{\bm{z}}\triangleq\frac{\bm{P}_{\tilde{\mathcal{V}}^{t}}\bm{z}}{\left\|\bm{P}_{\tilde{\mathcal{V}}^{t}}\bm{z}\right\|_{2}}$ and $\tilde{\Omega}\triangleq\tilde{\mathcal{V}}^{t}\cap\mathbb{S}^{n-1}$, thus $\tilde{\bm{z}}\in\tilde{\Omega}$. Note that $\tilde{\Omega}$ may differ from iteration to iteration but it is always a subset of $\mathcal{U}^{3}\cap\mathbb{S}^{n-1}$. Dividing (65) by $\left\|\bm{P}_{\tilde{\mathcal{V}}^{t}}\bm{z}\right\|_{2}$ we have

[TABLE]

Focusing on the last term, for a fixed error vector $\bm{e}$, note that $\bm{g}\triangleq\bm{A}^{*}\frac{\bm{e}}{\|\bm{e}\|_{2}}\sim\mathcal{N}(\mathbf{0},\bm{I}_{n})$. Applying standard results on concentration of supremum of Gaussian processes we get that

[TABLE]

holds for all $t$, with probability at least $1-\mathrm{e}^{-\frac{\eta^{2}}{2}}$. Using Lemma 2.3, (66) and (67) we have that

[TABLE]

holds for all t, with probability at least $1-3\mathrm{e}^{-\frac{\eta^{2}}{2}}$, where $\rho_{1}(\eta)$ and $\xi_{1}(\eta)$ are defined in (15) and (16). Plugging (68) in $\left\|\bm{z}\right\|_{2}^{2}=\left\|\bm{Q}_{\tilde{\mathcal{V}}^{t}}\bm{z}\right\|_{2}^{2}+\left\|\bm{P}_{\tilde{\mathcal{V}}^{t}}\bm{z}\right\|_{2}^{2}$ gives

[TABLE]

The left-hand side of the last equation is a second order polynomial of $\left\|\bm{z}\right\|_{2}$ with positive leading coefficient, hence it follows that $\left\|\bm{z}\right\|_{2}$ is upper bounded by the larger root. Equation (20) is obtained by calculating this bound and using the fact that $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for nonnegative $a$ and $b$.

Appendix C Proof of Lemma 3.4

Our proof is inspired by the proof of Lemma A.1 in [28]. The main differences between the proofs are in the later steps, where our arguments follow from concentration of measure, while those in [28] follow from a generalization of the RIP. We start with upper bounding the term $\|\bm{Q}_{\mathcal{V}_{\Delta}^{t}}(\bm{x}-\bm{x}^{t-1})\|_{2}$. To simplify the notations we define $\bm{v}\triangleq\bm{x}-\bm{x}^{t-1}$, $\overline{\mathcal{V}}\triangleq\mathcal{V}_{0}+\mathcal{V}^{t-1}$, and note that $\tilde{\bm{v}}=\bm{A}^{*}(\bm{y}-\bm{A}\bm{x}^{t-1})=\bm{A}^{*}\bm{A}\bm{v}+\bm{A}^{*}\bm{e}$. Letting $\mu_{2}$ and $\eta$ be positive parameters satisfying $\mu_{2}\leq\left(b_{m}+w(\mathcal{U}^{4}\cap\mathbb{S}^{n-1})+\eta\right)^{-2}$, we get

[TABLE]

The first inequality follows from the fact that $\bm{P}_{\mathcal{V}_{\Delta}^{t}}\bm{v}$ is the nearest to $\bm{v}$, in the $\ell_{2}$-norm sense, among all vectors in $\mathcal{V}_{\Delta}^{t}$. The last two inequalities follow from the triangle inequality. Focusing on the last term

[TABLE]

The first inequality follows from the fact that both $\mathcal{V}_{\Delta}^{t}$ and $\overline{\mathcal{V}}$ are subsets of $\mathcal{V}_{\Delta}^{t}+\overline{\mathcal{V}}$, and that by construction $\bm{P}_{\mathcal{V}_{\Delta}^{t}}\tilde{\bm{v}}$ is closer to $\tilde{\bm{v}}$ than $\bm{P}_{\overline{\mathcal{V}}}\tilde{\bm{v}}$ (since $\mathcal{V}_{\Delta}^{t},\overline{\mathcal{V}}\in\mathcal{S}^{2}$). For more details look at the discussion above (13) in [28]. The second inequality results from $\bm{P}_{\overline{\mathcal{V}}}\bm{P}_{\overline{\mathcal{V}}+\mathcal{V}_{\Delta}^{t}}\tilde{\bm{v}}$ being the nearest to $\bm{P}_{\overline{\mathcal{V}}+\mathcal{V}_{\Delta}^{t}}\tilde{\bm{v}}$ in the $\ell_{2}$-norm sense among all vectors in $\overline{\mathcal{V}}$, and the last inequality uses the triangle inequality. Note also that

[TABLE]

where the last inequality is due to fact that $\bm{I}_{n}-\bm{P}_{\overline{\mathcal{V}}}$ is an orthogonal projection and hence its spectral norm is bounded by 1. Using (70), (71) and (72), we get

[TABLE]

Dealing with the first term of (73), using $\bm{P}_{\overline{\mathcal{V}}+\mathcal{V}_{\Delta}^{t}}\bm{v}=\bm{v}$ and Lemma 2.3, we get

[TABLE]

where the last inequality holds for all $t$, with probability at least $1-2\mathrm{e}^{-\frac{\eta^{2}}{2}}$, and $\rho_{2}(\eta)$ is defined in (17). Focusing on the second term in (73), using the fact that $\bm{g}\triangleq\bm{A}^{*}\frac{\bm{e}}{\|\bm{e}\|_{2}}\sim\mathcal{N}(\mathbf{0},\bm{I}_{n})$ for a fixed error vector, and the fact that if $\bm{z}\in\mathbb{R}^{n}$ has unit Euclidean norm, then $\bm{P}_{\overline{\mathcal{V}}+\mathcal{V}_{\Delta}^{t}}\bm{z}\in\mathcal{U}^{4}\cap\mathbb{B}^{n}$, we get

[TABLE]

The inequality follows from standard results on concentration of supremum of Gaussian processes, similarly to (67), and holds for all $t$, with probability at least $1-\mathrm{e}^{-\frac{\eta^{2}}{2}}$. The last equality follows from the fact that for a linear subspace $\mathcal{L}$ we have $\underset{\bm{z}\in\mathcal{L}\cap\mathbb{B}^{n}}{\operatorname{sup}}\langle\bm{g},\bm{z}\rangle=\underset{\bm{z}\in\mathcal{L}\cap\mathbb{S}^{n-1}}{\operatorname{sup}}\langle\bm{g},\bm{z}\rangle$. In details, if $\bm{g}$ is orthogonal to $\mathcal{L}$ then this is trivial, else, the maximum is achieved at an extreme point of $\mathcal{L}\cap\mathbb{B}^{n}$, i.e., a point in $\mathcal{L}\cap\mathbb{S}^{n-1}$. Similarly, for a union of linear subspaces the maximum is achieved at an extreme point in the subspace that is closest to $\bm{g}$ in angular distance. We finish with the observation

[TABLE]

The first inequality uses the triangle inequality, the equality follows from $\tilde{\bm{x}}^{t},\bm{x}^{t-1}\in\tilde{\mathcal{V}}^{t}$, and the last inequality follows from $\mathcal{V}_{\Delta}^{t}\subset\tilde{\mathcal{V}}^{t}$. Combining (73) - (76), we get (21).