R Giryes

S Nam

M Elad

R Gribonval

M E Davies

Greedy-Like Algorithms for the Cosparse Analysis Model

Keywords: Sparse representations, Compressed sensing, Synthesis, Analysis, CoSaMP, Subspace-pursuit, Iterative hard threshodling, Hard thresholding pursuit. 2010 MSC: 94A20, 94A12, 62H12

The cosparse analysis model has been introduced recently as an interesting alternative to the standard sparse synthesis approach. A prominent question brought up by this new construction is the analysis pursuit problem -the need to find a signal belonging to this model, given a set of corrupted measurements of it. Several pursuit methods have already been proposed based on ℓ 1 relaxation and a greedy approach. In this work we pursue this question further, and propose a new family of pursuit algorithms for the cosparse analysis model, mimicking the greedy-like methods -compressive sampling matching pursuit (CoSaMP), subspace pursuit (SP), iterative hard thresholding (IHT) and hard thresholding pursuit (HTP). Assuming the availability of a near optimal projection scheme that finds the nearest cosparse subspace to any vector, we provide performance guarantees for these algorithms. Our theoretical study relies on a restricted isometry property adapted to the context of the cosparse analysis model. We explore empirically the performance of these algorithms by adopting a plain thresholding projection, demonstrating their good performance.

Introduction

Many natural signals and images have been observed to be inherently low dimensional despite their possibly very high ambient signal dimension. It is by now well understood that this phenomenon lies at the heart of the success of numerous methods of signal and image processing. Sparsity-based models for signals offer an elegant and clear way to enforce such inherent low-dimensionality, explaining their high popularity in recent years. These models consider the signal x ∈ R d as belonging to a finite union of subspaces of dimension k ≪ d [START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF]. In this paper we shall focus on one such approach -the cosparse analysis model -and develop pursuit methods for it.

Before we dive into the details of the model assumed and the pursuit problem, let us first define the following generic inverse problem that will accompany us throughout the paper: For some unknown signal x ∈ R d , an incomplete set of linear observations y ∈ R m (incomplete implies m < d) is available via

y = Mx + e, (1)
where e ∈ R m is an additive bounded noise that satisfies e 2 2 ≤ ǫ 2 . The task is to recover or approximate x. In the noiseless setting where e = 0, this amounts to solving y = Mx. Of course, a simple fact in linear algebra tells us that this problem admits infinitely many solutions (since m < d). Therefore, when all we have is the observation y and the measurement/observation matrix M ∈ R m×d , we are in a hopeless situation to recover x.

The Synthesis Approach

This is where 'sparse signal models' come into play. In the sparse synthesis model, the signal x is assumed to have a very sparse representation in a given fixed dictionary D ∈ R d×n . In other words, there exists α with few nonzero entries, as counted by the "ℓ 0 -norm" α 0 , such that x = Dα, and k := α 0 ≪ d.

(

) 2
Having this knowledge we solve (1) using xℓ 0 = D αℓ 0 , and αℓ 0 = argmin α α 0 subject to y -MDα 2 ≤ ǫ.

(

) 3
More details about the properties of this problem can be found in [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF][START_REF] Gribonval | Sparse representations in unions of bases[END_REF]. Since solving ( 3) is an NP-complete problem [START_REF] Davis | Adaptive greedy approximations[END_REF], approximation techniques are required for recovering x. One strategy is by using relaxation, replacing the ℓ 0 with ℓ 1 norm, resulting with the ℓ 1 -synthesis problem xℓ 1 = D αℓ 1 , and αℓ 1 = argmin α 1 s.t. y -MDα 2 ≤ ǫ.

(

) 4
For a unitary matrix D and a vector x with k-sparse representation α, if δ 2k < δ ℓ 1 then

xℓ 1 -x 2 ≤ C ℓ 1 e 2 , (5)
where xℓ 1 = D αℓ 1 , δ 2k is the constant of the restricted isometry property (RIP) of MD for 2k sparse signals, C ℓ 1 is a constant greater than √ 2 and δ ℓ 1 (≃ 0.4931) is a reference constant [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF][START_REF] Foucart | Sparse recovery algorithms: sufficient conditions in terms of restricted isometry constants[END_REF][START_REF] Mo | New bounds on the restricted isometry constant[END_REF]. Note that this result implies a perfect recovery in the absence of noise. The above statement was extended also for incoherent redundant dictionaries [START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF].

Another option for approximating (3) is using a greedy strategy, like in the thresholding technique or orthogonal matching pursuit (OMP) [START_REF] Pati | Orthonormal matching pursuit : recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]. A different related approach is the greedy-like family of algorithms. Among those we have compressive sampling matching pursuit (CoSaMP) [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], subspace pursuit (SP) [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], iterative hard thresholding (IHT) [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF] and hard thresholding pursuit (HTP) [START_REF] Foucart | Hard thresholding pursuit: an algorithm for compressive sensing[END_REF]. CoSaMP and SP were the first greedy methods shown to have guarantees in the form of (5) assuming δ 4k < δ CoSaMP and δ 3k ≤ δ SP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF][START_REF] Foucart | Sparse recovery algorithms: sufficient conditions in terms of restricted isometry constants[END_REF][START_REF] Giryes | RIP-based near-oracle performance guarantees for SP, CoSaMP, and IHT[END_REF]. Following their work, iterative hard thresholding (IHT) and hard thresholding pursuit (HTP) were shown to have similar guarantees under similar conditions [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF][START_REF] Foucart | Hard thresholding pursuit: an algorithm for compressive sensing[END_REF][START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF][START_REF] Foucart | Sparse recovery algorithms: sufficient conditions in terms of restricted isometry constants[END_REF]. Recently, a RIP based guarantee was developed also for OMP [START_REF] Zhang | Sparse recovery with orthogonal matching pursuit under RIP[END_REF].

The Cosparse Analysis Model

Recently, a new signal model called cosparse analysis model was proposed in [START_REF] Nam | The cosparse analysis model and algorithms[END_REF][START_REF] Nam | Cosparse analysis modeling -uniqueness and algorithms[END_REF]. The model can be summarized as follows: For a fixed analysis operator Ω ∈ R p×d referred to as the analysis dictionary, a signal x ∈ R d belongs to the cosparse analysis model with cosparsity ℓ if

ℓ := p -Ωx 0 . (6)
The quantity ℓ is the number of rows in Ω that are orthogonal to the signal. The signal x is said to be ℓ-cosparse, or simply cosparse. We denote the indices of the zeros of the analysis representation as the cosupport Λ and the submatrix that contains the rows from Ω that belong to Λ by Ω Λ . As the definition of cosparsity suggests, the emphasis of the cosparse analysis model is on the zeros of the analysis representation vector Ωx. This contrasts the emphasis on 'few non-zeros' in the synthesis model [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF]. It is clear that in the case where every ℓ rows in Ω are independent, x resides in a subspace of dimension dℓ that consists of vectors orthogonal to the rows of Ω Λ . In the general case where dependencies occur between the rows of Ω, the dimension is d minus the rank of Ω Λ . This is similar to the behavior in the synthesis case where a k-sparse signal lives in a k-dimensional space. Thus, for this model to be effective, we assume a large value of ℓ.

In the analysis model, recovering x from the corrupted measurements is done by solving the following minimization problem [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]:

x A-ℓ 0 = argmin x Ωx 0 subject to y -Mx 2 ≤ ǫ.

Solving this problem is NP-complete [START_REF] Nam | The cosparse analysis model and algorithms[END_REF], just as in the synthesis case, and thus approximation methods are required.

As before, we can use an ℓ 1 relaxation to [START_REF] Mo | New bounds on the restricted isometry constant[END_REF], replacing the ℓ 0 with ℓ 1 in [START_REF] Mo | New bounds on the restricted isometry constant[END_REF], resulting with the ℓ 1 -analysis problem [START_REF] Nam | The cosparse analysis model and algorithms[END_REF][START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF][START_REF] Vaiter | Robust sparse analysis regularization[END_REF]. Another option is the greedy approach. A greedy algorithm called Greedy Analysis Pursuit (GAP) has been developed in [START_REF] Nam | The cosparse analysis model and algorithms[END_REF][START_REF] Nam | Cosparse analysis modeling -uniqueness and algorithms[END_REF][START_REF] Nam | Cosparse analysis modeling, in: 9th International Conference on SamplingTheory and Applications (sampta-2011)[END_REF] that somehow mimics Orthogonal Matching Pursuit [START_REF] Pati | Orthonormal matching pursuit : recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] with a form of iterative reweighted least Squares (IRLS) [START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF]. Other alternatives for OMP, backward greedy (BG) and orthogonal BG (OBG), were presented in [START_REF] Rubinstein | Analysis K-SVD: A dictionary learning algorithm for the analysis sparse model[END_REF] for the case that M is the identity. For the same case, the parallel to the thresholding technique was analyzed in [START_REF] Peleg | Performance guarantees of the thresholding algorithm for the Co-Sparse analysis model[END_REF].

This Work

Another avenue exists for the development of analysis pursuit algorithms -constructing methods that will imitate the family of greedy-like algorithms. Indeed, we have recently presented preliminary and simplified versions of analysis IHT (AIHT), analysis HTP (AHTP), analysis CoSaMP (ACoSaMP) and Analysis SP (ASP) in [START_REF] Giryes | Iterative cosparse projection algorithms for the recovery of cosparse vectors[END_REF][START_REF] Giryes | CoSaMP and SP for the cosparse analysis model[END_REF] as analysis versions of the synthesis counterpart methods. This paper re-introduces these algorithms in a more general form, ties them to their synthesis origins, and analyze their expected performance. The main contribution of the paper is our result on the stability of these analysis pursuit algorithms. We show that after a finite number of iterations and for a given constant c 0 , the reconstruction result x of AIHT, AHTP, ACoSaMP and ASP all satisfy

x -x 2 ≤ c 0 e 2 , (8)
under a RIP-like condition on M and the assumption that we are given a good near optimal projection scheme. A bound is also given for the case where x is only nearly ℓ-cosparse. Similar results for the ℓ 1 analysis appear in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF][START_REF] Vaiter | Robust sparse analysis regularization[END_REF]. More details about the relation between these papers and our results will be given in Section 6. In addition to our theoretical results we demonstrate the performance of the four pursuit methods under a thresholding based simple projection scheme. Both our theoretical and empirical results show that linear dependencies in Ω that result with a larger cosparsity in the signal x, lead to a better reconstruction performance. This suggests that, as opposed to the synthesis case, strong linear dependencies within Ω are desired. This paper is organized as follows:

• In Section 2 we present the notation used along the paper.

• In Section 3 we define a RIP-like property, the Ω-RIP, for the analysis model, proving that it has similar characteristics like the regular RIP. In Section 4 the notion of near optimal projection is proposed and some nontrivial operators for which a tractable optimal projection exists are exhibited. Both the Ω-RIP and the near optimal projection are used throughout this paper as a main force for deriving our theoretical results.

• In Section 5 the four pursuit algorithms for the cosparse analysis framework are defined, adapted to the general format of the pursuit problem we have defined above.

• In Section 6 we derive the success guarantees for all the above algorithms in a unified way. Note that the provided results can be easily adapted to other union-of-subspaces models given near optimal projection schemes for them, in the same fashion done for IHT with an optimal projection scheme in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]. The relation between the obtained results and existing work appears in this section as well.

• Empirical performance of these algorithms is demonstrated in Section 7 in the context of the cosparse signal recovery problem. We use a simple thresholding as the near optimal projection scheme in the greedy-like techniques.

• Section 8 discuss the presented results and concludes our work.

Notations and Preliminaries

We use the following notation in our work:

• σ M is the largest singular value of M, i.e., σ 2 M = M * M 2 .
• • 2 is the euclidian norm for vectors and the spectral norm for matrices. • 1 is the ℓ 1 norm that sums the absolute values of a vector and • 0 , though not really a norm, is the ℓ 0 -norm which counts the number of non-zero elements in a vector.

• Given a cosupport set Λ, Ω Λ is a sub-matrix of Ω with the rows that belong to Λ.

• For given vectors v, z ∈ R d and an analysis dictionary Ω, cosupp(Ωv) returns the cosupport of Ωv and cosupp(Ωz, ℓ) returns the index set of ℓ smallest (in absolute value) elements in Ωz. If more than ℓ elements are zero all of them are returned. In the case where the ℓ-th smallest entry is equal to the ℓ + 1 smallest entry, one of them is chosen arbitrarily.

• In a similar way, in the synthesis case D T is a sub-matrix of D with columns1 corresponding to the set of indices T , supp(•) returns the support of a vector, supp(•, k) returns the set of k-largest elements and ⌈•⌉ k preserves the k-largest elements in a vector. In the case where the k-th largest entry is equal to the k + 1 largest entry, one of them is chosen arbitrarily.

• Q Λ = I -Ω †
Λ Ω Λ is the orthogonal projection onto the orthogonal complement of range(Ω * Λ ).

• P Λ = I -Q Λ = Ω † Λ Ω Λ is the orthogonal projection onto range(Ω * Λ ). • xAIHT /x AHTP /x ACoSaMP /x ASP
are the reconstruction results of AIHT/ AHTP/ ACoSaMP/ ASP respectively. Sometimes when it is clear from the context to which algorithms we refer, we abuse notations and use x to denote the reconstruction result.

• A cosupport Λ has a corank r if rank(Ω Λ ) = r. A vector v has a corank r if its cosupport has a corank r.

• [p] denotes the set of integers [1 . . . p]. • L Ω,ℓ = {Λ ⊆ [p], |Λ| ≥ ℓ} is the set of ℓ-cosparse cosupports and L corank Ω,r = {Λ ⊆ [p]
, rank(Ω Λ ) ≥ r} is the set of all cosupports with corresponding corank r.

• W Λ = span ⊥ (Ω Λ ) = {Q Λ z, z ∈ R d }
is the subspace spanned by a cosparsity set Λ.

• A Ω,ℓ = Λ∈L Ω,ℓ W Λ is the union of subspaces of ℓ-cosparse vectors and A corank Ω,r

= Λ∈L corank Ω,r
W Λ is the union of subspaces of all vectors with corank r. In the case that every ℓ rows of Ω are independent it is clear that A Ω,ℓ = A corank Ω,r . When it will be clear from the context, we will remove Ω from the subscript.

• x ∈ R d denotes the original unknown ℓ-cosparse vector and Λ x its cosupport.

• v, u ∈ A ℓ are used to denote general ℓ-cosparse vectors and z ∈ R d is used to denote a general vector.

Ω-RIP Definition and its Properties

We now turn to define the Ω-RIP, which parallels the regular RIP as used in [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF]. This property is a very important property for the analysis of the algorithms which holds for a large family of matrices M as we will see hereafter. Definition 3.1. A matrix M has the Ω-RIP property with a constant δ ℓ , if δ ℓ is the smallest constant that satisfies

(1 -δ ℓ ) v 2 2 ≤ Mv 2 2 ≤ (1 + δ ℓ ) v 2 2 , (9)
whenever Ωv has at least ℓ zeroes.

Note that though δ ℓ is also a function of Ω we abuse notation and use the same symbol for the Ω-RIP as the regular RIP. It will be clear from the context to which of them we refer and what Ω is in use with the Ω-RIP. A similar property that looks at the corank of the vectors can be defined Definition 3.2. A matrix M has the corank-Ω-RIP property with a constant δ corank r , if δ corank r is the smallest constant that satisfies

(1 -δ corank r ) u 2 2 ≤ Mu 2 2 ≤ (1 + δ corank r ) u 2 2 ( 10
)
whenever the corank of u with respect to Ω is greater or equal to r.

The Ω-RIP, like the regular RIP, inherits several key properties. We present only those related to δ ℓ , while very similar characteristics can be derived also for the corank-Ω-RIP. The first property we pose is an immediate corollary of the δ ℓ definition.

Corollary 3.3. If M satisfies the Ω-RIP with a constant δ ℓ then MQ Λ 2 2 ≤ 1 + δ ℓ ( 11
)
for any Λ ∈ L ℓ .

Proof: Any v ∈ A ℓ can be represented as v = Q Λ z with Λ ∈ L ℓ and z ∈ R d . Thus, the Ω-RIP in ( 9) can be reformulated as

(1 -δ ℓ ) Q Λ z 2 2 ≤ MQ Λ z 2 2 ≤ (1 + δ ℓ ) Q Λ z 2 2 ( 12
)
for any z ∈ R d and Λ ∈ L ℓ . Since Q Λ is a projection Q Λ z 2 2 ≤ z 2 2 .
Combining this with the right inequality in ( 12) gives

MQ Λ z 2 2 ≤ (1 + δ ℓ ) z 2 2 ( 13
)
for any z ∈ R d and Λ ∈ L ℓ . The first inequality in [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF] follows from ( 13) by the definition of the spectral norm.

Lemma 3.4. For l ≤ ℓ it holds that δ ℓ ≤ δl.

Proof: Since A ℓ ⊆ Al the claim is immediate.

Lemma 3.5. M satisfies the Ω-RIP if and only if

Q Λ (I -M * M)Q Λ 2 ≤ δ ℓ ( 14
)
for any Λ ∈ L ℓ .

Proof: The proof is similar to the one of the regular RIP as appears in [START_REF] Foucart | Sparse recovery algorithms: sufficient conditions in terms of restricted isometry constants[END_REF]. As a first step we observe that Definition 3.1 is equivalent to requiring

Mv 2 2 -v 2 2 ≤ δ ℓ v 2 2 ( 15
)
for any v ∈ A ℓ . The latter is equivalent to

MQ Λ z 2 2 -Q Λ z 2 2 ≤ δ ℓ Q Λ z 2 2 ( 16
)
for any set Λ ∈ L ℓ and any

z ∈ R d , since Q Λ z ∈ A ℓ . Next we notice that MQ Λ z 2 2 -Q Λ z 2 2 = z * Q Λ M * MQ Λ z -z * Q Λ z = Q Λ (M * M -I)Q Λ z, z . Since Q Λ (M * M -I)Q Λ is Hermitian we have that max z | Q Λ (M * M -I)Q Λ z, z | z 2 = Q Λ (M * M -I)Q Λ 2 . ( 17
)
Thus we have that Definition 3.1 is equivalent to [START_REF] Foucart | Hard thresholding pursuit: an algorithm for compressive sensing[END_REF] for any set Λ ∈ L ℓ .

Corollary 3.6. If M satisfies the Ω-RIP then

Q Λ 1 (I -M * M)Q Λ 2 2 ≤ δ ℓ , (18)
for any Λ 1 and

Λ 2 such that Λ 1 ∩ Λ 2 ∈ L ℓ . Proof: Since Λ 1 ∩ Λ 2 ⊆ Λ 1 and Λ 1 ∩ Λ 2 ⊆ Λ 2 Q Λ 1 (I -M * M)Q Λ 2 2 ≤ Q Λ 2 ∩Λ 1 (I -M * M)Q Λ 2 ∩Λ 1 2 .
Using Lemma 3.5 completes the proof.

As we will see later, we require the Ω-RIP to be small. Thus, we are interested to know for what matrices this hold true. In the synthesis case, where Ω is unitary and the Ω-RIP is identical to the RIP, it was shown for certain family of random matrices, such as matrices with Bernoulli or Subgaussian ensembles, that for any value of [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF][START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF][START_REF] Mendelson | Uniform uncertainty principle for bernoulli and subgaussian ensembles[END_REF], where δ k is the RIP constant and C ǫ k is a constant depending on ǫ k and M. A similar result for the same family of random matrices holds for the analysis case. The result is a special case of the result presented in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF].

ǫ k if m ≥ C ǫ k k log( m kǫ k ) then δ k ≤ ǫ k
Theorem 3.7 (Theorem 3.3 in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]). Let M ∈ R m×d be a random matrix such that for any z ∈ R d and 0 < ǫ ≤ 1 3 it satisfies

P Mz 2 2 -z 2 2 ≥ ǫ z 2 2 ≤ e -C M mǫ 2 , ( 19
)
where C M > 0 is a constant. For any value of

ǫ ℓ > 0, if m ≥ 32 C M ǫ 2 r log( L corank r ) + (d -r) log(9/ǫ r ) + t , (20)
then δ corank r ≤ ǫ r with probability exceeding 1e -t .

The above theorem is important since it shows that the Ω-RIP holds with a small constant for a large family of matrices -the same family that satisfy the RIP property. In a recent work it was even shown that by randomizing the signs of the columns in the matrices that satisfy the RIP we get new matrices that also satisfy the RIP [START_REF] Krahmer | New and improved Johnson-Lindenstrauss embeddings via the restricted isometry property[END_REF]. Thus, requiring the Ω-RIP constant to be small, as will be done hereafter, is legitimate.

For completeness we present a proof for theorem 3.7 in Appendix A based on [START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF][START_REF] Mendelson | Uniform uncertainty principle for bernoulli and subgaussian ensembles[END_REF][START_REF] Baraniuk | A simple proof of the restricted isometry property for random matrices[END_REF]. We include in it also the proof of Theorem 3.8 to follow. In the case that Ω is in general position L corank r = p r ≤ ( ep p-r ) p-r (inequality is by Stirling's formula) and thus m ≥ (pr) log( ep p-r ). Since we want m to be smaller than d we need pℓ to be smaller than d. This limits the size of p for Ω since r cannot be greater than d. Thus, we present a variation of the theorem which states the results in terms of δ ℓ and ℓ instead of δ corank r and r. The following theorem is also important because of the fact that our theoretical results are in terms of δ ℓ and not δ corank r . It shows that δ ℓ is small in the same family of matrices that guarantees δ corank r to be small. Theorem 3.8. Under the same setup of Theorem 3.7, for any

ǫ ℓ > 0 if m ≥ 32 C M ǫ 2 ℓ (p -ℓ) log 9p (p -ℓ)ǫ ℓ + t , (21)
then δ ℓ ≤ ǫ ℓ with probability exceeding 1e -t .

Remark that when Ω is in general position ℓ cannot be greater than d and thus p cannot be greater than 2d [START_REF] Nam | The cosparse analysis model and algorithms[END_REF]. For this reason, if we want to have large values for p we should allow linear dependencies between the rows of Ω. In this case the cosparsity of the signal can be greater than d. This explains why linear dependencies are a favorable thing in analysis dictionaries [START_REF] Rubinstein | Analysis K-SVD: A dictionary learning algorithm for the analysis sparse model[END_REF]. In Section 7 we shall see that also empirically we get a better recovery when Ω contains linear dependencies.

Near Optimal Projection

As we will see hereafter, in the proposed algorithms we will face the following problem: Given a general vector z ∈ R d , we would like to find an ℓ-cosparse vector that is closest to it in the ℓ 2 -norm sense. In other words, we would like to project the vector to the closest ℓ-cosparse subspace. Given the cosupport Λ of this space the solution is simply Q Λ z. Thus, the problem of finding the closest ℓ-cosparse vector turns to be the problem of finding the cosupport of the closest ℓ-cosparse subspace. We denote the procedure of finding this cosupport by

S * ℓ (z) = argmin Λ∈L ℓ z -Q Λ z 2 2 . ( 22
)
In the representation domain in the synthesis case, the support of the closest k-sparse subspace is found simply by hard thresholding, i.e., taking the support of the k-largest elements. However, in the analysis case calculating ( 22) is NP-complete with no efficient method for doing it for a general Ω [START_REF] Gribonval | Projection onto the k-cosparse set is NP-hard[END_REF]. Thus an approximation procedure Ŝℓ is needed. For this purpose we introduce the definition of a near-optimal projection [START_REF] Giryes | Iterative cosparse projection algorithms for the recovery of cosparse vectors[END_REF].

Definition 4.1. A procedure Ŝℓ implies a near-optimal projection Q Ŝℓ (•) with a constant C ℓ if for any z ∈ R d z -Q Ŝℓ (z) z 2 2 ≤ C ℓ z -Q S * ℓ (z) z 2 2 . (23)
A clear implication of this definition is that if Ŝℓ implies a near-optimal projection with a constant C ℓ then for any vector z ∈ R d and an ℓ-cosparse

vector v ∈ R d z -Q Ŝℓ (z) z 2 2 ≤ C ℓ z -v 2 2 . (24)
Similarly to the Ω-RIP, the above discussion can be directed also for finding the closest vector with corank r defining S corank * r and near optimal projection for this case in a very similar way to [START_REF] Vaiter | Robust sparse analysis regularization[END_REF] and Definition 4.1 respectively. Having a near-optimal cosupport selection scheme for a general operator is still an open problem and we leave it for a future work. It is possible that this is also NP-complete. We start by describing a simple thresholding rule that can be used with any operator. Even though it does not have any known (near) optimality guarantee besides the case of unitary operators, the numerical section will show it performs well in practice. Then we present two tractable algorithms for finding the optimal cosupport for two non-trivial analysis operators, the one dimensional finite difference operator Ω 1D-DIF [START_REF] Han | Optimal segmentation of signals and its application to image denoising and boundary feature extraction[END_REF] and the fused Lasso operator Ω FUS [START_REF] Tibshirani | Sparsity and smoothness via the fused Lasso[END_REF].

Later in the paper, we propose theoretical guarantees for algorithms that use operators that has an optimal or a near-optimal cosupport selection scheme. We leave the theoretical study of the thresholding technique for a future work but demonstrate its performance empirically in Section 7 where this rule is used showing that also when nearoptimality is not at hand reconstruction is feasible.

Cosupport Selection by Thresholding

One intuitive option for cosupport selection is the simple thresholding

Ŝℓ (z) = cosupp(Ωz, ℓ), (25)
which selects as a cosupport the indices of the ℓ-smallest elements after applying Ω on z. As mentioned above, this selection method is optimal for unitary analysis operators where it coincides with the hard thresholding used in synthesis. However, in the general case this selection method is not guaranteed to give the optimal cosupport. Its near optimality constant C ℓ is not close to one and is equal to the fraction of the largest and smallest eigenvalues (which are not zero) of the submatrices composed of ℓ rows from Ω [START_REF] Giryes | Iterative cosparse projection algorithms for the recovery of cosparse vectors[END_REF].

One example for an operator for which the thresholding is sub-optimal is the 1D-finite difference operator Ω 1D-DIF . This operator is defined as: Figure 1: Comparison between projection using thresholding cosupport selection and optimal cosupport selection. As it can be seen the thresholding projection error is much larger than the optimal projection error by a factor much larger than 1

Ω 1D-DIF =                     -1 1 • • • . . . -1 1
In this case, given a signal z, applying Ω 1D-DIF on it, result with a vector of coefficients that represents the differences in the signal. The thresholding selection method will select the indices of the ℓ smallest elements in Ωz as the cosupport Λ z . For example, for the signal z ∈ R 201 in Fig 1(a) that contains 100 times one, 100 times minus one and 1.5 as the last element, the thresholding will select the cosupport to be the first 199 coefficients in Ω 1D-DIF z that appears in Fig 1(b) and thus the projected vector will be the one in Fig 1(c). Its error in the ℓ 2 -norm sense is √ 200. However, selecting the cosupport to be the first 99 elements and last 100 elements result with the projected vector in Fig. 1(d), which has a smaller projection error (2.5). Thus, it is clear that the thresholding is sub-optimal for Ω 1D-DIF . In a similar way it is also sub-optimal for the 2D-finite difference operator Ω 2D-DIF that returns the vertical and horizontal differences of a two dimensional signal. Though not optimal, the use of thresholding with this operator is illustrated in Section 7 demonstrating that also when a good projection is not at hand, good reconstruction is still possible.

Optimal Analysis Projection Operators

As mentioned above, in general it would appear that determining the optimal projection is computationally difficult with the only general solution being to fully enumerate the projections onto all possible cosupports. Here we highlight two cases where it is relatively easy (polynomial complexity) to calculate the optimal cosparse projection.

Case 1: 1D finite difference

For the 1D finite difference operator the analysis operator is not redundant (p = d -1) but neither is it invertible. As we have seen, a simple thresholding does not provide us with the optimal cosparse projection. Thus, in order to determine the best ℓ-cosparse approximation for a given vector z we take another route and note that we are looking for the closest (in the ℓ 2 -norm sense to z) piecewise constant vector with pℓ change-points. This problem has been solved previously in the signal processing literature using dynamic programming (DP), see for example: [START_REF] Han | Optimal segmentation of signals and its application to image denoising and boundary feature extraction[END_REF]. Thus for this operator it is possible to calculate the best cosparse representation in O(d 2 ) operations. The existence of a DP solution follows from the ordered localized nature of the finite difference operator. To the best of our knowledge, there is no known extension to 2D finite difference.

Case 2: Fused Lasso Operator

A redundant operator related to the 1D finite difference operator is the so-called fused Lasso operator, usually used with the analysis ℓ 1 -minimization [START_REF] Tibshirani | Sparsity and smoothness via the fused Lasso[END_REF]. This usually takes the form:

Ω FUS = Ω 1D-DIF ǫI . (27)
Like Ω 1D-DIF this operator works locally and therefore we can expect to derive a DP solution to the approximation problem. This is presented below. Remark 4.2. Note that in terms of the cosparsity model the ǫ parameter plays no role. This is in contrast to the traditional convex optimization solutions where the value of ǫ is pivotal [START_REF] Vaiter | Robust sparse analysis regularization[END_REF]. It is possible to mimic the ǫ dependence within the cosparsity framework by considering a generalized fused Lasso operator of the form:

Ω ǫFUS =                       Ω 1D-DIF Ω 1D-DIF . . . Ω 1D-DIF I                       . ( 28
)
where the number of repetitions of the Ω 1D-DIF operator (and possibly the I operator) can be selected to mimic a weight on the number of nonzero coefficients of each type. For simplicity we only consider the case indicated by (27)

A recursive solution to the optimal projector for Ω FUS

Rather than working directly with the operator Ω FUS we make use of the following observation. An ℓ-cosparse vector v (or k-sparse vector) for Ω FUS is a piecewise constant vector with k 1 change points and k 2 non-zero entries such that k 1 + k 2 = k = pℓ, where p = 2d -1. To understand better the relation between k 1 and k 2 , notice that k 1 = 0 implies equality of all entries, so k 2 = 0 or d, hence ℓ = p or d -1. Conversely, considering d ≤ ℓ < p or 0 ≤ ℓ < d -1 implies k 1 0. It also implies that there is at least one nonzero value, hence k 2 0.

Thus, an ℓ-cosparse vector v for Ω FUS can be parameterized in terms of a set of change points, {n i } i=0:k 1 +1 , and a set of constants, {µ i } i=1:k 1 +1 , such that:

v j = µ i , n i-1 < j ≤ n i ( 29
)
with the convention that n 0 = 0 and n k 1 +1 = d, unless stated otherwise. We will also make use of the indicator vector, s, defined as:

s i =        0 if µ i = 0, 1 otherwise for 1 ≤ i ≤ k 1 + 1. (30)
Using this alternative parametrization we can write the minimum distance between a vector z and the set of k-sparse fused Lasso coefficients as:

F k (z) = min 1≤k 1 ≤k min {n i } i=1:k 1 {µ i } i=1:k 1 +1 n k 1 <d k 1 +1 i=1 n i j=n i-1 +1 (z j -µ i ) 2 , subject to k 1 +1 i=1 s i (n i -n i-1 ) = k -k 1 (31)
Although this looks a formidable optimization task we now show that it can be computed recursively through a standard DP strategy, modifying the arguments in [START_REF] Han | Optimal segmentation of signals and its application to image denoising and boundary feature extraction[END_REF].

Let us define the optimal cost, I k (L, ω, k 1 ), for the vector [z 1 , . . . , z L ] T with k 1 change points and s k 1 +1 = ω, as:

I k (L, ω, k 1 ) = min {n i } i=1:k 1 {s i } i=1:k 1 +1 n k 1 <L, n k 1 +1 =L s k 1 +1 =ω k 1 +1 i=1 n i j=n i-1 +1 (z j -µ i ) 2 , subject to k 1 +1 i=1 s i (n i -n i-1 ) = k -k 1 and µ i = s i n i -n i-1 n i l=n i-1 +1 z l (32)
where we have set µ i to the optimal sample means. Notice that calculating I k (L, ω, k 1 ) is easy for k 1 ≤ k ≤ 1. Thus, we calculate it recursively considering two separate scenarios:

Case 1: ω = 0 where the last block of coefficients are zero. This gives:

I k (L, 0, k 1 ) = min n k 1 <L                           L j=n k 1 +1 (z j ) 2 + min {n i } i=1:k 1 -1 {s i } i=1:k 1 -1 n k 1 -1 <n k 1 s k 1 =1 k 1 i=1 n i j=n i-1 +1 (z j -µ i ) 2                           , subject to k 1 i=1 s i (n i -n i-1 ) = (k -1) -(k 1 -1)
and

µ i = s i n i -n i-1 n i l=n i-1 +1 z l , (33)
(noting that if s k 1 +1 = 0 then s k 1 = 1 since otherwise n k 1 would not have been a change point). This simplifies to the recursive formula:

I k (L, 0, k 1 ) = min n k 1 <L          L j=n k 1 +1 (z j ) 2 + I k-1 (n k 1 , 1, k 1 -1)          (34)
Case 2: ω = 1 when the final block of coefficients are non-zero we have:

I k (L, 1, k 1 ) = min n k 1 <L n k 1 +1 =L s k 1                   L j=n k 1 +1 (z j -µ k 1 +1 ) 2 + min {n i } i=1:k 1 -1 {s i } i=1:k 1 -1 n k 1 -1 <n k 1 k 1 i=1 n i j=n i-1 +1 (z j -µ i ) 2                   , subject to k 1 i=1 s i (n i -n i-1 ) = (k -L + n k 1 -1) -(k 1 -1)
and

µ i = s i n i -n i-1 n i l=n i-1 +1 z l . (35)
This simplifies to the recursive relationship:

I k (L, 1, k 1 ) = min n k 1 <L s k 1          L j=n k 1 +1 (z j -µ k 1 +1 ) 2 + I k-L+n k 1 -1 (n k 1 , s k 1 , k 1 -1)          subject to µ k 1 +1 = L l=n k 1 +1 z l / L -n k 1 (36)
Equations ( 34) and ( 36) are sufficient to enable the calculation of the optimal projection in polynomial time,starting with k 1 ≤ k ≤ 1 and recursively evaluating the costs for k

≥ k 1 ≥ 1. Finally, we have F k (z) = min k 1 ≤k,ω∈{0,1} I k (d, ω, k 1 ).
The implementation details are left as an exercise for the reader.

New Analysis algorithms

Quick Review of the Greedy-Like Methods

Before we turn to present the analysis versions of the greedy-like techniques we recall their synthesis versions. These use a prior knowledge about the cardinality k and actually aim at approximating a variant of (3)

argmin α y -MDα 2 2 subject to α 0 ≤ k. ( 37
)
For simplicity we shall present the greedy-like pursuits for the case D = I. In the general case M should be replaced with MD, x with α and the reconstruction result should be x = D α. In addition, in the algorithms' description we do not specify the stopping criterion. Any standard stopping criterion, like residual's size or relative iteration change, can be used. More details can be found in [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF].

IHT and HTP: IHT [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF] and HTP [START_REF] Foucart | Hard thresholding pursuit: an algorithm for compressive sensing[END_REF] are presented in Algorithm 1. Each IHT iteration is composed of two basic steps. The first is a gradient step, with a step size µ t , in the direction of minimizing y -Mx 2 2 . The step size can be either constant in all iterations (µ t = µ) or changing [START_REF] Kyrillidis | Recipes on hard thresholding methods[END_REF]. The result vector x g is not guaranteed to be sparse and thus the second step of IHT projects x g to the closest k-sparse subspace by keeping its largest k elements. The HTP takes a different strategy in the projection step. Instead of using a simple projection to the closest k-sparse subspace, HTP selects the vector in this subspace that minimizes y -Mx 2 2 [START_REF] Foucart | Hard thresholding pursuit: an algorithm for compressive sensing[END_REF][START_REF] Blumensath | Accelerated iterative hard thresholding[END_REF]. CoSaMP and SP: CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF] and SP [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] are presented in Algorithm 2. The difference between these two techniques is similar to the difference between IHT and HTP. Unlike IHT and HTP, the estimate for the support of x in each CoSaMP and SP iteration is computed by observing the residual y t resid = y -Mx t . In each iteration, CoSaMP and SP extract new support indices from the residual by taking the indices of the largest elements in M * y t resid . They add the new indices to the estimated support set from the previous iteration creating a new estimated support T t with cardinality larger than k. Having the updated support, in a similar way to the projection in HTP, an objective aware projection is performed resulting with an estimate w for x that is supported on T t . Since we know that x is k-sparse we want to project w to a k-sparse subspace. CoSaMP does it by simple hard thresholding like in IHT. SP does it by an objective aware projection similar to HTP.

Algorithm

Analysis greedy-like methods

Given the synthesis greedy-like pursuits, we would like to define their analysis counterparts. For this task we need to 'translate' each synthesis operation into an analysis one. This gives us a general recipe for converting algorithms between the two schemes. The parallel lines between the schemes are presented in Table 1. Those become more intuitive and clear when we keep in mind that while the synthesis approach focuses on the non-zeros, the analysis concentrates on the zeros.

For clarity we dwell a bit more on the equivalences. For the cosupport selection, as mentioned in Section 4, computing the optimal cosupport is a combinatorial problem and thus the approximation Ŝℓ is used. Having a selected cosupport Λ, the projection to its corresponding cosparse subspace becomes trivial, given by Q Λ .

Given two vectors

v 1 ∈ A ℓ 1 and v 2 ∈ A ℓ 2 such that Λ 1 = cosupp(Ωv 1 ) and Λ 2 = cosupp(Ωv 2 ), we know that |Λ 1 | ≥ ℓ 1 and |Λ 2 | ≥ ℓ 2 . Denoting T 1 = supp(Ωv 1 ) and T 2 = supp(Ωv 2 ) it is clear that supp(Ω(v 1 + v 1 )) ⊆ T 1 ∪ T 2 . Noticing that supp(•) = cosupp(•) C it is clear that |T 1 | ≤ p -ℓ 1 , |T 2 | ≤ p -ℓ 2 and cosupp(Ω(v 1 + v 2 )) ⊇ (T 1 ∪ T 2 ) C = T C 1 ∩ T C 2 = Λ 1 ∩ Λ 2 .
From the last equality we can also deduce that

|Λ 1 ∩ Λ 2 | = p -|T 1 ∪ T 2 | ≥ p -(p -ℓ 1 ) -(p -ℓ 2 ) = ℓ 1 + ℓ 2 -p.
With the above observations we can develop the analysis versions of the greedy-like algorithms. As in the synthesis case, we do not specify a stopping criterion. Any stopping criterion used for the synthesis versions can be used also for the analysis ones. Initialize the support T 0 = ∅, the residual y 0 resid = y and set t = 0. while halting criterion is not satisfied do t = t + 1.

Find new support elements:

T ∆ = supp(M * y t-

Synthesis operation name

Synthesis operation Analysis operation name Analysis operation

Support selection

Largest k elements:

T = supp(•, k)
Cosupport selection Using a near optimal projection:

Λ = Ŝℓ (•) Orthogonal Projection of z to a k-sparse subspace with support T z T Orthogonal projection of z to an ℓ-cosparse subspace with cosupport Λ Q Λ z
Objective aware projection to a k-sparse subspace with support T

M † T y = argmin v y -Mv 2 2 s.t. v T C = 0 Objective aware projection to an ℓ-cosparse subspace with cosupport Λ argmin v y -Mv 2 2 s.t. Ω Λ v = 0 Support of v 1 + v 2 where supp(v 1 ) = T 1 and supp(v 2 ) = T 2 supp(v 1 + v 2 ) ⊆ T 1 ∪ T 2 Cosupport of v 1 + v 2 where cosupp(v 1 ) = Λ 1 and cosupp(v 2 ) = Λ 2 cosupp(v 1 + v 2 ) ⊇ Λ 1 ∩ Λ 2 Maximal size of T 1 ∪ T 2 where |T 1 | ≤ k 1 and |T 2 | ≤ k 2 |T 1 ∪ T 2 | ≤ k 1 + k 2 Minimal size of Λ 1 ∩ Λ 2 where |Λ 1 | ≥ ℓ 1 and |Λ 2 | ≥ ℓ 2 |Λ 1 ∩ Λ 2 | ≥ ℓ 1 + ℓ 2 -p
Table 1: Parallel synthesis and analysis operations AIHT and AHTP: Analysis IHT (AIHT) and analysis HTP (AHTP) are presented in Algorithm 3. As in the synthesis case, the choice of the gradient stepsize µ t is crucial: If µ t 's are chosen too small, the algorithm gets stuck at a wrong solution and if too large, the algorithm diverges. We consider two options for µ t .

In the first we choose µ t = µ for some constant µ for all iterations. A theoretical discussion on how to choose µ properly is given in Section 6.1.

The second option is to select a different µ in each iteration. One way for doing it is to choose an 'optimal' stepsize µ t by solving the following problem

µ t := argmin µ y -Mx t 2 2 . ( 38
)
Since Λt = Ŝℓ (x t-1 + µ t M * (y -Mx t-1 )) and xt = Q Λt (x g ), the above requires a line search over different values of µ and along the search Λt might change several times. A simpler way is an adaptive step size selection as proposed in [START_REF] Kyrillidis | Recipes on hard thresholding methods[END_REF] for IHT. In a heuristical way we limit the search to the cosupport Λ = Ŝℓ (M * (y -Mx t-1 )) ∩ Λt-1 . This is the intersection of the cosupport of xt-1 with the ℓ-cosparse cosupport of the estimated closest ℓ-cosparse subspace to M * (y -Mx t-1 ). Since xt-1 = Q Λ xt-1 , finding µ turns to be

µ t := argmin µ y -M(x t-1 + µQ ΛM * (y -Mx t-1 )) 2 2 . ( 39
)
Algorithm 3 Analysis Iterative hard thresholding (AIHT) and analysis hard thresholding pursuit (AHTP) Require: ℓ, M, Ω, y where y = Mx + e, ℓ is the cosparsity of x under Ω and e is the additive noise. Ensure: xAIHT or xAHTP : ℓ-cosparse approximation of x.

Initialize estimate x0 = 0 and set t = 0. while halting criterion is not satisfied do t = t + 1. Perform a gradient step: This procedure of selecting µ t does not require a line search and it has a simple closed form solution.

x g = xt-1 + µ t M * (y -Mx t-
To summarize, there are three main options for the step size selection:

• Constant step-size selection -uses a constant step size µ t = µ in all iterations.

• Optimal changing step-size selection -uses different values for µ t in each iterations by minimizing y -Mx t 2 .

• Adaptive changing step-size selection -uses [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF].

ACoSaMP and ASP: analysis CoSaMP (ACoSaMP) and analysis SP (ASP) are presented in Algorithm 4. The stages are parallel to those of the synthesis CoSaMP and SP. We dwell a bit more on the meaning of the parameter a in the algorithms. This parameter determines the size of the new cosupport Λ ∆ in each iteration. a = 1 means that the size is ℓ and according to Table 1 it is equivalent to a = 1 in the synthesis as done in SP in which we select new k indices for the support in each iteration. In synthesis CoSaMP we use a = 2 and select 2k new elements. 2k is the maximal support size of two added k-sparse vectors. The corresponding minimal size in the analysis case is 2ℓp according to Table 1. For this setting we need to choose a = 2ℓ-p ℓ .

The Unitary Case

For Ω = I the synthesis and the analysis greedy-like algorithms become equivalent. This is easy to see since in this case we have

p = d, k = d -ℓ, Λ = T C , Q Λ x = x T and T 1 ∪ T 2 = Λ 1 ∩ Λ 2 for Λ 1 = T C 1 and Λ 2 = T C 2 .
In addition, Ŝℓ = S * ℓ finds the closest ℓ-cosparse subspace by simply taking the smallest ℓ elements. Using similar arguments, also in the case where Ω is a unitary matrix the analysis methods coincide with the synthesis ones. In order to get exactly the same algorithms M is replaced with MΩ * in the synthesis techniques and the output is multiplied by Ω * .

Based on this observation, we can deduce that the guarantees of the synthesis greedy-like methods apply also for the analysis ones in a trivial way. Thus, it is tempting to assume that the last should have similar guarantees based on the Ω-RIP. In the next section we develop such claims.

Relaxed Versions for High Dimensional Problems

Before moving to the next section we mention a variation of the analysis greedy-like techniques. In AHTP, ACoSaMP and ASP we need to solve the constrained minimization problem min x y -Mx 2 2 s.t. Ω Λ x 2 2 = 0. For high dimensional signals this problem is hard to solve and we suggest to replace it with minimizing y -Mx 2 2 + λ Ω Λ x 2 2 , where λ is a relaxation constant. This results in a relaxed version of the algorithms. We refer hereafter to these versions as relaxed AHTP (RAHTP) relaxed ASP (RASP) and relaxed ACoSaMP (RACoSaMP).

Algorithms Guarantees

In this section we provide theoretical guarantees for the reconstruction performance of the analysis greedy-like methods. For AIHT and AHTP we study both the constant step-size and the optimal step-size selections. For ACoSaMP and ASP the analysis is made for a = 2ℓ-p ℓ , but we believe that it can be extended also to other values of a, such as a = 1. The performance guarantees we provide are summarized in the following two theorems. The first theorem, for AIHT and AHTP, is a simplified version of Theorem 6.5 and the second theorem, for ASP and ACoSaMP, is a combination of Corollaries 6.9 and 6.14, all of which appear hereafter along with their proofs. Before presenting the theorems we recall the problem we aim at solving: Definition 6.1 (Problem P). Consider a measurement vector y ∈ R m such that y = Mx + e where x ∈ R d is ℓcosparse, M ∈ R m×d is a degradation operator and e ∈ R m is a bounded additive noise. The largest singular value of M is σ M and its Ω-RIP constant is δ ℓ . The analysis operator Ω ∈ R p×d is given and fixed. A procedure Ŝℓ for finding a cosupport that implies a near optimal projection with a constant C ℓ is assumed to be at hand. Our task is to recover x from y. The recovery result is denoted by x. Theorem 6.2 (Stable Recovery of AIHT and AHTP). Consider the problem P and apply either AIHT or AHTP with a certain constant step-size or an optimal changing step-size, obtaining xt after t iterations. If

(C ℓ -1)σ 2 M C ℓ < 1 ( 40
)
and δ 2ℓ-p < δ 1 (C ℓ , σ 2 M ), where δ 1 (C ℓ , σ 2
M ) is a constant guaranteed to be greater than zero whenever [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF] is satisfied and C ℓ is the near-optimal projection constant for cosparsity ℓ (Definition 4.1), then after a finite number of iterations t

* x -xt * 2 ≤ c 1 e 2 , ( 41
)
implying that these algorithms lead to a stable recovery. The constant c 1 is a function of δ 2ℓ-p , C ℓ and σ 2 M , and the constant step-size used is dependent on δ 1 (C ℓ , σ 2 M ).

Theorem 6.3 (Stable Recovery of ASP and ACoSaMP). Consider the problem P and apply either ACoSaMP or ASP with a = 2ℓ-p ℓ , obtaining xt after t iterations. If

(C 2 Ŝ -1)σ 2 M C 2 Ŝ < 1, ( 42
)
and

δ 4ℓ-3p < δ 2 (C Ŝ, σ 2 M ), 14
where C Ŝ = max(C ℓ , C 2ℓ-p ) and δ 2 (C Ŝ, σ 2 M
) is a constant guaranteed to be greater than zero whenever (42) is satisfied, then after a finite number of iterations t

* x -xt * 2 ≤ c 2 e 2 , ( 43
)
implying that these algorithms lead to a stable recovery. The constant c 2 is a function of δ 4ℓ-3p , C ℓ , C 2ℓ-p and σ2 M . Before we proceed to the proofs, let us comment on the constants in the above theorems. Their values can be calculated using Theorem 6.5, and Corollaries 6.9 and 6.14. In the case where Ω is a unitary matrix, [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF] and ( 42) are trivially satisfied since C ℓ = C 2ℓ-p = 1. In this case the Ω-RIP conditions become δ 2ℓ-p < δ 1 (1, σ 2 M ) = 1/3 for AIHT and AHTP, and δ 4ℓ-3p < δ 2 (1, σ 2 M ) = 0.0156 for ACoSaMP and ASP. In terms of synthesis RIP for MΩ * , the condition δ 2ℓ-p < 1/3 parallels δ 2k (MΩ * ) < 1/3 and similarly δ 4ℓ-3p < 0.0156 parallels δ 4k (MΩ * ) < 0.0156. Note that the condition we pose for AIHT and AHTP in this case is the same as the one presented for synthesis IHT with a constant step size [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF]. Better reference constants were achieved in the synthesis case for all four algorithms and thus we believe that there is still room for improvement of the reference constants in the analysis context.

In the non-unitary case, the value of σ M plays a vital role, though we believe that this is just an artifact of our proof technique. For a random Gaussian matrix whose entries are i.i.d with a zero-mean and a variance 1 m , σ M behaves like

d m 1 + d m .
This is true also for other types of distributions for which the fourth moment is known to be bounded [START_REF] Rudelson | Non-asymptotic theory of random matrices: extreme singular values[END_REF]. For example, for d/m = 1.5 we have found empirically that σ 2 M ≃ 5. In this case we need C ℓ ≤ 5 4 for (40) to hold and C Ŝ ≤ 1.118 for (42) to hold, and both are quite demanding on the quality of the near-optimal projection. For C ℓ = C Ŝ = 1.05 we have the conditions δ 2ℓ-p ≤ 0.289 for AIHT and AHTP, and δ 4ℓ-3p ≤ 0.0049 for ACoSaMP and ASP; and for C ℓ = C Ŝ = 1.1 we have δ 2ℓ-p ≤ 0.24 for AIHT and AHTP, and δ 4ℓ-3p ≤ 0.00032 for ACoSaMP and ASP.

As in the synthesis case, the Ω-RIP requirements for the theoretical bounds of AIHT and AHTP are better than those for ACoSaMP and ASP. In addition, in the migration from the synthesis to the analysis we lost more precision in the bounds for ACoSaMP and ASP than in those of AIHT and AHTP. In particular, even in the case where Ω is the identity we do not coincide with any of the synthesis parallel RIP reference constants. We should also remember that the synthesis bound for SP is in terms of δ 3k and not δ 4k [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF]. Thus, we expect that it will be possible to give a condition for ASP in terms of δ 3ℓ-2p with better reference constants. However, our main interest in this work is to show the existence of such bounds, and in Section 6.5 we dwell more on their meaning.

We should note that here and elsewhere we can replace the conditions on δ 2ℓ-p and δ 4ℓ-3p in the theorems to conditions on δ corank 2r-p and δ corank 4r-3p and the proofs will be almost the same 2 . In this case we will be analyzing a version of the algorithms which is driven by the corank instead of the cosparsity. This would mean we need the near-optimal projection to be in terms of the corank. In the case where Ω is in a general position, there is no difference between the cosparsity ℓ and the corank r. However, when we have linear dependencies in Ω the two measures differ and an ℓ-cosparse vector is not necessarily a vector with a corank r.

As we will see hereafter, our recovery conditions require δ 2ℓ-p and δ 4ℓ-3p to be as small as possible and for this we need 2ℓp and 4ℓ -3p to be as large as possible. Thus, we need ℓ to be as close as possible to p and for highly redundant Ω this cannot be achieved without having linear dependencies in Ω. Apart from the theoretical advantage of linear dependencies in Ω, we also show empirically that an analysis dictionary with linear dependencies has better recovery rate than analysis dictionary in a general position of the same dimension. Thus, we deduce that linear dependencies in Ω lead to better bounds and restoration performance.

Though linear dependencies allow ℓ to be larger than d and be in the order of p, the value of the corank is always bounded by d and cannot be expected to be large enough for highly redundant analysis dictionaries. In addition, we will see hereafter that the number of measurements m required by the Ω-RIP is strongly dependent on ℓ and less effected by the value of r. From the computational point of view we note also that using corank requires its computation in each iteration which increases the overall complexity of the algorithms. Thus, it is more reasonable to have conditions on δ 2ℓ-p and δ 4ℓ-3p than on δ corank 2r-p and δ corank 4r-3p , and our study will be focused on the cosparsity based algorithms.

AIHT and AHTP Guarantees

A uniform guarantee for AIHT in the case that an optimal projection is given, is presented in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]. The work in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF] dealt with a general union of subspaces, A, and assumed that M is bi-Lipschitz on the considered union of subspaces. In our case A = A ℓ and the bi-Lipschitz constants of M are the largest B L and smallest B U where 0 < B L ≤ B U such that for all ℓ-cosparse vectors v 1 , v 2 :

B L v 1 + v 2 2 2 ≤ M(v 1 + v 2 ) 2 2 ≤ B U v 1 + v 2 2 2 . ( 44
)
Under this assumption, one can apply Theorem 2 from [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF] to the idealized AIHT that has access to an optimal projection and uses a constant step size µ t = µ. Relying on Table 1 we present this theorem and replace B L and B U with 1δ 2ℓ-p and 1 + δ 2ℓ-p respectively. Theorem 6.4 (Theorem 2 in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]). Consider the problem P with C ℓ = 1 and apply AIHT with a constant step size µ.

If 1 + δ 2ℓ-p ≤ 1 µ < 1.5(1 -δ 2ℓ-p ) then after a finite number of iterations t * x -xt * 2 ≤ c 3 e 2 , ( 45
)
implying that AIHT leads to a stable recovery. The constant c 3 is a function of δ 2ℓ-p and µ.

In this work we extend the above in several ways: First, we refer to the case where optimal projection is not known, and show that the same flavor guarantees apply for a near-optimal projection3 . The price we seemingly have to pay is that σ M enters the game. Second, we derive similar results for the AHTP method. Finally, we also consider the optimal step size and show that the same performance guarantees hold true in that case. Theorem 6.5. Consider the problem P and apply either AIHT or AHTP with a constant step size µ or an optimal changing step size. For a positive constant η > 0, let b

1 := η 1 + η and b 2 := (C ℓ -1)σ 2 M b 2 1 C ℓ (1 -δ 2ℓ-p ) . Suppose b 2 b 2 1 = (C ℓ -1)σ 2 M C ℓ (1-δ 2ℓ-p ) < 1, 1 + δ 2ℓ-p ≤ 1 µ < 1 + 1 -b 2 b 2 1 b 1 (1 -δ 2ℓ-p ) and 1 µ ≤ σ 2 M . Then for t ≥ t * log η e 2 2 y 2 2 log (1 + 1 η ) 2 ( 1 µ(1-δ 2ℓ-p ) -1)C ℓ + (C ℓ -1)(µσ 2 M -1) + C ℓ η 2 , ( 46
) x -xt 2 2 ≤ (1 + η) 2 1 -δ 2ℓ-p e 2 2 , (47)
implying that AIHT and AHTP lead to a stable recovery. Note that for an optimal changing step-size we set µ = 1 1+δ 2ℓ-p in t * and the theorem conditions turn to be b

2 b 2 1 < 1 and 1 + δ 2ℓ-p < (1 + 1 -b 2 b 2 1 )b 1 (1 -δ 2ℓ-p ).
This theorem is the parallel to Theorems 2.1 in [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF] for IHT. A few remarks are in order for the nature of the theorem, especially in regards to the constant η. One can view that η gives a trade-off between satisfying the theorem conditions and the amplification of the noise. In particular, one may consider that the above theorem proves the convergence result for the noiseless case by taking η to infinity; one can imagine solving the problem P where e → 0, and applying the theorem with appropriately chosen η which approaches infinity. It is indeed possible to show that the iterate solutions of AIHT and AHTP converges to x when there is no noise. However, we will not give a separate proof since the basic idea of the arguments is the same for both cases.

As to the minimal number of iterations t * given in (46), one may ask whether it can be negative. In order to answer this question it should be noted that according to the conditions of the Theorem the term inside the log in the denominator (46) is always greater than zero. Thus, t * will be negative only if y 2 2 < η e 2 2 . Indeed, in this case 0 iterations suffice for having the bound in (47).

The last remark is on the step-size selection. The advantage of the optimal changing step-size over the constant step-size is that we get the guarantee of the optimal constant step-size µ = 1 1+δ 2ℓ-p without computing it. This is important since in practice we cannot evaluate the value of δ 2ℓ-p . However, the disadvantage of using the optimal changing step-size is its additional complexity for the algorithm. Thus, one option is to approximate the optimal selection rule by replacing it with an adaptive one, for which we do not have a theoretical guarantee. Another option is to set µ = 6/5 which meets the theorem conditions for small enough δ 2ℓ-p , in the case where an optimal projection is at hand.

We will prove the theorem by proving two key lemmas first. The proof technique is based on ideas from [START_REF] Garg | Gradient descent with sparsification: an iterative algorithm for sparse recovery with restricted isometry property[END_REF] and [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]. Recall that the two iterative algorithms try to reduce the objective y -Mx t 2 2 over iterations t. Thus, the progress of the algorithms can be indirectly measured by how much the objective y -Mx t 2

2 is reduced at each iteration t. The two lemmas that we present capture this idea. The first lemma is similar to Lemma 3 in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF] and relates y -Mx t 2 2 to y -Mx t-1 2 2 and similar quantities at iteration t -1. We remark that the constraint 1 µ ≤ σ 2 M in Theorem 6.5 may not be necessary and is added only for having a simpler derivation of the results in this theorem.

Furthermore, this is a very mild condition compared to 1

µ < 1 + 1 -b 2 b 2 1 b 1 (1 -δ 2ℓ-p
) and can only limit the range of values that can be used with the constant step size versions of the algorithms. Lemma 6.6. Consider the problem P and apply either AIHT or AHTP with a constant step size µ satisfying 1 µ ≥ 1 + δ 2ℓ-p or an optimal step size. Then, at the t-th iteration, the following holds:

y -Mx t 2 2 -y -Mx t-1 2 2 ≤ C ℓ y -Mx 2 2 -y -Mx t-1 2 2 (48) + C ℓ 1 µ(1 -δ 2ℓ-p ) -1 M(x -xt-1 ) 2 2 + (C ℓ -1)µσ 2 M y -Mx t-1 2 2 .
For the optimal step size the bound is achieved with the value µ = 1 1+δ 2ℓ-p .

The proof of the above lemma appears in Appendix B. The second lemma is built on the result of Lemma 6.6. It shows that once the objective y -Mx t-1 2 2 at iteration t -1 is small enough, then we are guaranteed to have small y -Mx t 2 2 as well. Given the presence of noise, this is quite natural; one cannot expect it to approach 0 but may expect it not to become worse. Moreover, the lemma also shows that if y -Mx t-1 2 2 is not small, then the objective in iteration t is necessarily reduced by a constant factor. Lemma 6.7. Suppose that the same conditions of Theorem 6.5 hold true. If

y -Mx t-1 2 2 ≤ η 2 e 2 2 , then y -Mx t 2 2 ≤ η 2 e 2 2 . Furthermore, if y -Mx t-1 2 2 > η 2 e 2 2 , then y -Mx t 2 2 ≤ c 4 y -Mx t-1 2 2 ( 49
)
where

c 4 := 1 + 1 η 2 1 µ(1 -δ 2ℓ-p ) -1 C ℓ + (C ℓ -1)(µσ 2 M -1) + C ℓ η 2 < 1.
Having the two lemmas above, the proof of the theorem is straightforward.

Proof:[Proof of Theorem 6.5] When we initialize x0 = 0, we have y -Mx 0 2 2 = y 2 2 . Assuming that y 2 > η e 2 and applying Lemma 6.7 repeatedly, we obtain

y -Mx t 2 2 ≤ max(c t 4 y 2 2 , η 2 e 2 2 )
.

Since c t 4 y 2 2 ≤ η 2 e 2 2 for t ≥ t * , we have simply

y -Mx t 2 2 ≤ η 2 e 2 2 ( 50
)
for t ≥ t * . If y -Mx 0 2 = y 2 ≤ η e 2 then according to Lemma 6.7, (50) holds for every t > 0. Finally, we observe

x -xt 2 2 ≤ 1 1 -δ 2ℓ-p M(x -xt ) 2 2
(51) and, by the triangle inequality,

M(x -xt ) 2 ≤ y -Mx t 2 + e 2 . ( 52
)
By plugging (50) into (52) and then the resulted inequality into (51), the result of the Theorem follows.

As we have seen, the above AIHT and AHTP results hold for the cases of using a constant or an optimal changing step size. The advantage of using an optimal one is that we do not need to find µ that satisfies the conditions of the theorem -the knowledge that such a µ exists is enough. However, its disadvantage is the additional computational complexity it introduces. In Section 5 we have introduced a third option of using an approximated adaptive step size. In the next section we shall demonstrate this option in simulations, showing that it leads to the same reconstruction result as the optimal selection method. Note, however, that our theoretical guarantees do not cover this case.

ACoSaMP Guarantees

Having the results for AIHT and AHTP we turn to ACoSaMP and ASP. We start with a theorem for ACoSaMP. Its proof is based on the proof for CoSaMP in [START_REF] Foucart | Sparse recovery algorithms: sufficient conditions in terms of restricted isometry constants[END_REF]. Theorem 6.8. Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . Let C Ŝ = max(C ℓ , C 2ℓ-p ) and suppose that there exists γ > 0 such that

(1 + C Ŝ) 1 - C Ŝ (1 + γ) 2 -(C Ŝ -1)σ 2 M < 1. ( 53
)
Then, there exists δ ACoSaMP (C Ŝ, σ 2 M , γ) > 0 such that, whenever δ 4ℓ-3p ≤ δ ACoSaMP (C Ŝ, σ 2 M , γ), the t-th iteration of the algorithm satisfies

x -xt 2 ≤ ρ 1 ρ 2 x -xt-1 2 + (η 1 + ρ 1 η 2 ) e 2 , ( 54
)
where

η 1 2+C ℓ 1+C ℓ + 2 √ C ℓ + C ℓ 1 + δ 3ℓ-2p 1 -δ 4ℓ-3p , η 2 2 1 + δ 3ℓ-2p γ(1 + α) + (1 + δ 2ℓ-p )C 2ℓ-p γ(1 + α)(1 + γ) + (C 2ℓ-p -1)(1 + γ)σ 2 M (1 + α)(1 + γ)γ , ρ 2 1 1 + 2δ 4ℓ-3p √ C ℓ + C ℓ 1 -δ 2 4ℓ-3p , ρ 2 2 1 -δ 4ℓ-3p - C 2ℓ-p (1 + γ) 2 1 -δ 2ℓ-p 2 -(C 2ℓ-p -1)(1 + δ 2ℓ-p )σ 2 M 2 and α = δ 4ℓ-3p C 2ℓ-p (1+γ) 2 1 -δ 2ℓ-p 2 -(C 2ℓ-p -1)(1 + δ 2ℓ-p )σ 2 M -δ 4ℓ-3p . Moreover, ρ 2 1 ρ 2 2 < 1, i.e.
, the iterates converges.

The constant γ plays a similar role to the constant η of Theorem 6.5. It gives a tradeoff between satisfying the theorem conditions and the noise amplification. However, as opposed to η, the conditions for the noiseless case are achieved when γ tends to zero. An immediate corollary of the above theorem is the following. Corollary 6.9. Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . If (53) holds and δ 4ℓ-3p < δ ACoSaMP (C Ŝ, σ 2 M , γ), where C Ŝ and γ are as in Theorem 6.8 and δ ACoSaMP (C Ŝ, σ 2 M , γ) is a constant guaranteed to be greater than zero whenever (42) is satisfied, then for any

t ≥ t * = log( x 2 / e 2 ) log(1/ρ 1 ρ 2 ) , x -xt * ACoSaMP 2 ≤ 1 + 1 -(ρ 1 ρ 2 ) t * 1 -ρ 1 ρ 2 (η 1 + ρ 1 η 2 ) e 2 , ( 55
)
implying that ACoSaMP leads to a stable recovery. The constants η 1 , η 2 , ρ 1 and ρ 2 are the same as in Theorem 6.8.

Proof: By using (54) and recursion we have that

x -xt * ACoSaMP 2 ≤ (ρ 1 ρ 2 ) t * x -x0 ACoSaMP 2
(56)

+(1 + ρ 1 ρ 2 + (ρ 1 ρ 2 ) 2 + . . . (ρ 1 ρ 2 ) t * -1 ) (η 1 + ρ 1 η 2 ) e 2 .

Since x0

ACoSaMP = 0, after t * iterations, one has

(ρ 1 ρ 2 ) t * x -x0 ACoSaMP 2 = (ρ 1 ρ 2 ) t * x 2 ≤ e 2 . ( 57
)
By using the equation of geometric series with (56) and plugging (57) into it, we get (55).

We turn now to prove the theorem. Instead of presenting the proof directly, we divide the proof into several lemmas. The first lemma gives a bound for xw 2 as a function of e 2 and P Λt (xw) 2 . Lemma 6.10. Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . For each iteration we have

x -w 2 ≤ 1 1 -δ 2 4ℓ-3p P Λt (x -w) 2 + 1 + δ 3ℓ-2p 1 -δ 4ℓ-3p e 2 . ( 58
)
The second lemma bounds xxt ACoSaMP 2 in terms of P Λt (xxt ACoSaMP ) 2 and e 2 using the first lemma.

Lemma 6.11. Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . For each iteration we have

x -xt 2 ≤ ρ 1 P Λt (x -w) 2 + η 1 e 2 , ( 59
)
where η 1 and ρ 1 are the same constants as in Theorem 6.8.

The last lemma bounds P Λt (xw) 2 with xxt-1 ACoSaMP 2 and e 2 .

Lemma 6.12. Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . if

C 2ℓ-p < σ 2 M (1 + γ) 2 σ 2 M (1 + γ) 2 -1 , ( 60
)
then there exists δ ACoSaMP (C 2ℓ-p , σ 2 M , γ) > 0 such that for any δ 2ℓ-p < δ ACoSaMP (C 2ℓ-p , σ 2 M , γ)

P Λt (x -w) 2 ≤ η 2 e 2 + ρ 2 x -xt-1 2 . ( 61
)
The constants η 2 and ρ 2 are as defined in Theorem 6.8.

The proofs of Lemmas 6.10, 6.11 and 6.12 appear in Appendix D, Appendix E and Appendix F respectively. With the aid of the above three lemmas we turn to prove Theorem 6.8.

Proof:[Proof of Theorem 6.8] Remark that since 1 + C Ŝ > 1 we have that (53) implies

C Ŝ (1+γ) 2 -(C Ŝ -1)σ 2 M ≥ 0.
Because of that the condition in (60) in Lemma 6.12 holds. Substituting the inequality of Lemma 6.12 into the inequality of Lemma 6.11 gives (54). The iterates convergence if

ρ 2 1 ρ 2 2 = 1+2δ 4ℓ-3p √ C ℓ +C ℓ 1-δ 2 4ℓ-3p ρ 2 2 < 1. By noticing that ρ 2 2 < 1 it is enough to require 1+C ℓ 1-δ 2 4ℓ-3p ρ 2 2 + 2δ 4ℓ-3p √ C ℓ 1-δ 2 4ℓ-3p < 1. The last is equivalent to (1 + C ℓ )           1 -          δ 4ℓ-3p - C 2ℓ-p (1 + γ) 2 1 -δ 2ℓ-p 2 -(C 2ℓ-p -1)(1 + δ 2ℓ-p )σ 2 M          2           (62) + 2δ 4ℓ-3p C ℓ -1 + δ 2 4ℓ-3p < 0.
It is easy to verify that ζ(C, δ)

C (1+γ) 2 1 - √ δ 2 -(C -1)(1 + δ)σ 2 M is a decreasing function of both δ and C for 0 ≤ δ ≤ 1 and C > 1. Since 1 ≤ C 2ℓ-p ≤ C Ŝ, δ 2ℓ-p ≤ δ 4ℓ-3p and δ ≥ 0 we have that ζ(C Ŝ, δ 4ℓ-3p ) ≤ ζ(C 2ℓ-p , δ 4ℓ-3p ) ≤ ζ(C 2ℓ-p , δ 2ℓ-p ) ≤ ζ(1, 0) = 1 (1+γ) 2 ≤ 1. Thus we have that -1 ≤ -( δ 4ℓ-3p -ζ(C 2ℓ-p , δ 2ℓ-p )) 2 ≤ -δ 4ℓ-3p + 2 δ 4ℓ-3p -ζ(C Ŝ, δ 4ℓ-3p
). Combining this with the fact that C ℓ ≤ C Ŝ provides the following guarantee for

ρ 2 1 ρ 2 2 < 1, ( 1
+ C Ŝ) 1 -δ 4ℓ-3p + 2 δ 4ℓ-3p (63)
- C Ŝ (1 + γ) 2 1 -2 δ 4ℓ-3p + δ 4ℓ-3p + (C Ŝ -1)(1 + δ 4ℓ-3p )σ 2 M + 2δ 4ℓ-3p C Ŝ -1 + δ 2 4ℓ-3p < 0.
Let us now assume that δ 4ℓ-3p ≤ 1 2 . This necessarily means that δ ACoSaMP ≤ 1 2 in the end. This assumption implies δ 2 4ℓ-3p ≤ 1 2 δ 4ℓ-3p . Using this and gathering coefficients, we now consider the condition

(1 + C Ŝ) 1 - C Ŝ (1 + γ) 2 + (C Ŝ -1)σ 2 M -1 + 2(1 + C Ŝ) 1 + C Ŝ (1 + γ) 2 δ 4ℓ-3p (64)
+ (1 + C Ŝ) -1 - C Ŝ (1 + γ) 2 + (C Ŝ -1)σ 2 M + 2 C Ŝ + 1 2 δ 4ℓ-3p < 0.
The expression on the LHS is a quadratic function of δ 4ℓ-3p . Note that since (53) holds the constant term in the quadratic function is negative. This guarantees the existence of a range of values [0, δ ACoSaMP (C Ŝ, σ 2 M , γ)] for δ 4ℓ-3p for which (64) holds, where δ ACoSaMP (C Ŝ, σ 2 M , γ) is the square of the positive solution of the quadratic function. In case of two positive solutions we should take the smallest among them -in this case the coefficient of δ 4ℓ-3p in (64) will be positive.

Looking back at the proof of the theorem, we observe that the value of the constant δ ACoSaMP (C Ŝ, σ 2 M , γ) can potentially be improved: at the beginning of the proof, we have used ρ 2 2 ≤ 1. By the end, we obtained ρ 2 2 ≤ ρ -2 1 ≤ 0.25 since ρ 1 > 2. If we were to use this bound at the beginning, we would have obtained better constant δ ACoSaMP (C Ŝ, σ 2 M , γ).

ASP Guarantees

Having the result of ACoSaMP we turn to derive a similar result for ASP. The technique for deriving a result for ASP based on the result of ACoSaMP is similar to the one we used to derive a result for AHTP from the result of AIHT. Theorem 6.13. Consider the problem P and apply ASP with a = 2ℓ-p ℓ . If (53) holds and δ 4ℓ-3p ≤ δ ASP (C Ŝ, σ 2 M , γ), where C Ŝ and γ are as in Theorem 6.8, and δ ASP (C Ŝ, σ 2 M , γ) is a constant guaranteed to be greater than zero whenever (53) is satisfied, then the t-th iteration of the algorithm satisfies

x -xt ASP 2 ≤ 1 + δ 2ℓ-p 1 -δ 2ℓ-p ρ 1 ρ 2 x -xt-1 ASP 2 + 1 + δ 2ℓ-p 1 -δ 2ℓ-p (η 1 + ρ 1 η 2 ) + 2 1 -δ 2ℓ-p e 2 . ( 65
)
and the iterates converges, i.e., ρ 2 1 ρ 2 2 < 1. The constants η 1 , η 2 , ρ 1 and ρ 2 are the same as in Theorem 6.8.

Proof: We first note that according to the selection rule of xASP we have that

y -Mx t ASP 2 ≤ y -MQ Λt w 2 . (66)
Using the triangle inequality and the fact that y = Mx + e for both the LHS and the RHS we have

M(x -xt ASP ) 2 -e 2 ≤ M(x -Q Λt w) 2 + e 2 .
Using the Ω-RIP property of M with the fact that x, xASP and Q Λt w are ℓ-cosparse we have

x -xt ASP 2 ≤ 1 + δ 2ℓ-p 1 -δ 2ℓ-p x -Q Λt w 2 + 2 1 -δ 2ℓ-p e 2 .
Noticing that Q Λt w is the solution we get in one iteration of ACoSaMP with initialization of xt-1 ASP , we can combine the above with the result of Theorem 6.8 getting (65). For

1+δ 2ℓ-p 1-δ 2ℓ-p ρ 1 ρ 2 < 1 to hold we need that 1 + 2δ 4ℓ-3p √ C ℓ + C ℓ (1 -δ 4ℓ-3p ) 2
Remark that the above differs from what we have for ACoSaMP only in the denominator of the first element in the LHS. In ACoSaMP 1δ 2 4ℓ-3p appears instead of (1δ 4ℓ-3p ) 2 . Thus, Using a similar process to the one in the proof of ACoSaMP we can show that (67) holds if the following holds

(1 + C Ŝ) 1 - C Ŝ (1 + γ) 2 + (C Ŝ -1)σ 2 M -1 + 2(1 + C Ŝ) 1 + C Ŝ (1 + γ) 2 δ 4ℓ-3p (68)
+ (1 + C Ŝ) -1 - C Ŝ (1 + γ) 2 + (C Ŝ -1)σ 2 M + 2 C Ŝ + 2 δ 4ℓ-3p < 0.
Notice that the only difference of the above compared to (64) is that we have +2 instead of +0.5 in the coefficient of δ 4ℓ-3p and this is due to the difference we mentioned before in the denominator in (67). The LHS of (68) is a quadratic function of δ 4ℓ-3p . As before, we notice that if (53) holds then the constant term of the above is positive and thus δ ASP (C Ŝ, σ 2 M , γ) ≥ 0 exists and is the square of the positive solution of the quadratic function.

Having Theorem 6.13 we can immediately have the following corollary which is similar to the one we have for ACoSaMP. The proof resembles the one of Corollary 6.9 and omitted. Corollary 6.14. Consider the problem P and apply ASP with a = 2ℓ-p ℓ . If (53) holds and δ 4ℓ-3p ≤ δ ASP (C Ŝ, σ 2 M , γ), where C Ŝ and γ are as in Theorem 6.8, and δ ASP (C Ŝ, σ 2 M , γ) is a constant guaranteed to be greater than zero whenever (42) is satisfied, then for any

t ≥ t * =             log( x 2 / e 2 ) log(1/ 1+δ 2ℓ-p 1-δ 2ℓ-p ρ 1 ρ 2 )             , x t ASP -x 2 ≤ 1 + 1 - 1+δ 2ℓ-p 1-δ 2ℓ-p ρ 1 ρ 2 t 1 - 1+δ 2ℓ-p 1-δ 2ℓ-p ρ 1 ρ 2 • 1 + δ 2ℓ-p 1 -δ 2ℓ-p (η 1 + ρ 1 η 2 ) + 2 1 -δ 2ℓ-p e 2 . ( 69
)
implying that ASP leads to a stable recovery. The constants η 1 , η 2 , ρ 1 and ρ 2 are the same as in Theorem 6.8.

Non-Exact Cosparse Case

In the above guarantees we have assumed that the signal x is ℓ-cosparse. In many cases, it is not exactly ℓ-cosparse but only nearly so. Denote by

x ℓ = Q S * ℓ (x)
x the best ℓ-cosparse approximation of x, we have the following theorem that provides us with a guarantee also for this case. Similar result exists also in the synthesis case for the synthesis-ℓ 1 minimization problem [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF]. Theorem 6.15. Consider a variation of problem P where x is a general vector, and apply either AIHT or AHTP both with either constant or changing step size; or ACoSaMP or ASP with a = 2ℓ-p ℓ , and all are used with a zero initialization. Under the same conditions of Theorems 6.2 and 6.3 we have for any t

≥ t * x -x 2 ≤ x -x ℓ 2 + c M(x -x ℓ ) 2 + c e 2 , ( 70
)
where t * and c are the constants from Theorems 6.2 and 6.3.

Proof: First we notice that we can rewrite y = Mx ℓ + M(xx ℓ ) + e. Denoting e ℓ = M(xx ℓ ) + e we can use Theorems 6.2 and 6.3 to recover x ℓ and have

x ℓ -x 2 ≤ c e ℓ 2 . (71)
Using the triangle inequality for xx 2 with the above gives

x -x 2 ≤ x -x ℓ 2 + x ℓ -x 2 ≤ x -x ℓ 2 + c e ℓ 2 . (72)
Using again the triangle inequality for e ℓ 2 ≤ e 2 + M(xx ℓ ) 2 provides us with the desired result.

Theorem Conditions

Having the results of the theorems we ask ourselves whether their conditions are feasible. As we have seen in Section 3, the requirement on the Ω-RIP for many non-trivial matrices. In addition, as we have seen in the introduction of this section we need C ℓ and C 2ℓ-p to be one or close to one for satisfying the conditions of the theorems. Using the thresholding in [START_REF] Rubinstein | Analysis K-SVD: A dictionary learning algorithm for the analysis sparse model[END_REF] for cosupport selection with a unitary Ω satisfies the conditions in a trivial way since C ℓ = C 2ℓ-p = 1. This case coincides with the synthesis model for which we already have theoretical guarantees. As shown in Section 4, optimal projection schemes exist for Ω 1D-DIF and Ω FUS which do not belong to the the synthesis framework. For a general Ω, a general projection scheme is not known and if the thresholding method is used the constants in [START_REF] Rubinstein | Analysis K-SVD: A dictionary learning algorithm for the analysis sparse model[END_REF] do not equal one and are not even expected to be close to one [START_REF] Giryes | Iterative cosparse projection algorithms for the recovery of cosparse vectors[END_REF]. It is interesting to ask whether there exists an efficient general projection scheme that guarantees small constants for any given operator Ω, or for specifically structured Ω. We leave these questions as subject for future work. Instead, we show empirically in the next section that a weaker projection scheme that does not fulfill all the requirements of the theorems leads to a good reconstruction result. This suggests that even in the absence of good near optimal projections we may still use the algorithms practically.

Comparison to Other Works

Among the existing theoretical works that studied the performance of analysis algorithms [START_REF] Nam | The cosparse analysis model and algorithms[END_REF][START_REF] Vaiter | Robust sparse analysis regularization[END_REF][START_REF] Peleg | Performance guarantees of the thresholding algorithm for the Co-Sparse analysis model[END_REF], the result that resembles ours is the result for ℓ 1 -analysis in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF]. This work analyzed the ℓ 1 -analysis minimization problem with a synthesis perspective. The analysis dictionary Ω was replaced with the conjugate of a synthesis dictionary D which is assumed to be a tight frame, resulting with the following minimization problem.

xA-ℓ 1 = argmin z D * z 1 s.t. y -Mz 2 ≤ ǫ. (73)
It was shown that if M has the D-RIP [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF][START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF] with δ 7k < 0.6, an extension of the synthesis RIP, then

xA-ℓ 1 -x 2 ≤ Cℓ 1 ǫ + D * x -[D * x] k 1 √ k . ( 74
)
We say that a matrix M has a D-RIP with a constant δ k if for any signal z that has a k-sparse representation under D

(1 -δ k ) z 2 2 ≤ Mz 2 2 ≤ (1 + δ k ) z 2 2 . (75)
The authors in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] presented this result as a synthesis result that allows linear dependencies in D at the cost of limiting the family of signals to be those for which D * x -[D * x] k 1 is small. However, having the analysis perspective, we can realize that they provided a recovery guarantee for ℓ 1 -analysis under the new analysis model for the case that Ω is a tight frame. An easy way to see it is to observe that for an ℓ-cosparse signal x, setting k = pℓ, we have that Ωx -[Ω * x] p-ℓ 1 = 0 and thus in the case ǫ = 0 we get that (74) guarantees the recovery of x by using (73) with D * = Ω. Thus, though the result in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] was presented as a reconstruction guarantee for the synthesis model, it is actually a guarantee for the analysis model.

Our main difference from [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] is that the proof technique relies on the analysis model and not on the synthesis one and that the results presented here are for general operators and not only for tight frames. For instance, the operators Ω 1D-DIF and Ω FUS for which the guarantees hold are not tight frames where Ω 1D-DIF is not even a frame. However, the drawback of our approached compared to the work in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] is that it is still not known how to perform an optimal or a near optimal projection for a tight frame.

In the non-exact sparse case our results differ from the one in (74) in the sense that it looks at the projection error and not at the values of Ωx. It would be interesting to see if there is a connection between the two and whether one implies the other.

A recent work has studied the ℓ 1 -analysis minimization with the 2D-DIF operator, also known as anisotropic two dimensional total-variation (2D-TV) [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF]. It would be interesting to see whether similar results can be achieved for the greedy-like techniques proposed here with 2D-DIF.

Experiments

In this section we repeat some of the experiments performed in [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] for the noiseless case (e = 0) and some of the experiments performed in [START_REF] Nam | Cosparse analysis modeling, in: 9th International Conference on SamplingTheory and Applications (sampta-2011)[END_REF] for the noisy case4 .

Targeted Cosparsity

Just as in the synthesis counterpart of the proposed algorithms, where a target sparsity level k must be selected before running the algorithms, we have to choose the targeted cosparsity level which will dictate the projection steps. In the synthesis case it is known that it may be beneficial to over-estimate the sparsity k. Similarly in the analysis framework the question arises: In terms of recovery performance, does it help to under-estimate the cosparsity ℓ? A tentative positive answer comes from the following heuristic: Let Λ be a subset of the cosupport Λ of the signal x with l := | Λ| < ℓ = |Λ|. According to Proposition 3 in [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] κ Ω ( l) ≤ m 2 (76) is a sufficient condition to identify Λ in order to recover x from the relations y = Mx and Ω Λx = 0. κ Ω ( l) = max Λ∈Ll dim(W Λ) is a function of l. Therefore, we can replace ℓ with the smallest l that satisfies (76) as the effective cosparsity in the algorithms. Since it is easier to identify a smaller cosupport set it is better to run the algorithm with the smallest possible value of l, in the absence of noise. In the presence of noise, larger values of ℓ allows a better denoising. Note, that in some cases the smallest possible value of l will be larger than the actual cosparsity of x. In this case we cannot replace ℓ with l. We take two examples for selecting l. The first is for Ω which is in general position and the second is for Ω 2D-DIF , the finite difference analysis operator that computes horizontal and vertical discrete derivatives of an image which is strongly connected to the total variation (TV) norm minimization as noted before. For Ω that is in general position κ Ω ( l) = max(dℓ, 0) [START_REF] Nam | The cosparse analysis model and algorithms[END_REF]. In this case we choose

l = min d - m 2 , ℓ . (77)
For

Ω DIF we have κ Ω DIF ( l) ≥ d -ℓ 2 -ℓ 2 -1 [18] and l = ⌈min((-1/ √ 2 + √ 2d -m -1.5) 2 , ℓ)⌉. (78)
Replacing ℓ with l is more relevant to AIHT and AHTP than ACoSaMP and ASP since in the last we intersect cosupport sets and thus the estimated cosupport set need to be large enough to avoid empty intersections. Thus, for Ω in general position we use the true cosparsity level for ACoSaMP and ASP. For Ω DIF , where linear dependencies occur, the corank does not equal the cosparsity and we use l instead of ℓ since it will be favorable to run the algorithm targeting a cosparsity level in the middle. In this case ℓ tends to be very large and it is more likely to have non-empty intersections .

Phase Diagrams for Synthetic Signals in the Noiseless Case

We begin with with synthetic signals in the noiseless case. We test the performance of AIHT with a constant stepsize, AIHT with an adaptive changing step-size, AHTP with a constant step-size, AHTP with an adaptive changing step-size, ACoSaMP with a = 2ℓ-p ℓ , ACoSaMP with a = 1, ASP with a = 2ℓ-p ℓ and ASP with a = 1. We compare the results to those of A-ℓ 1 -minimization [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] and GAP [START_REF] Nam | The cosparse analysis model and algorithms[END_REF]. We use a random matrix M and a random tight frame with d = 120 and p = 144, where each entry in the matrices is drawn independently from the Gaussian distribution.

We draw a phase transition diagram [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF] for each of the algorithms. We test 20 different possible values of m and 20 different values of ℓ and for each pair repeat the experiment 50 times. In each experiment we check whether we have a perfect reconstruction. White cells in the diagram denotes a perfect reconstruction in all the experiments of the pair and black cells denotes total failure in the reconstruction. The values of m and ℓ are selected according to the following formula:

m = δd ℓ = d -ρm, (79)
where δ, the sampling rate, is the x-axis of the phase diagram and ρ, the ratio between the cosparsity of the signal and the number of measurements, is the y-axis. Figure 2 presents the reconstruction results of the algorithms. It should be observed that AIHT and AHTP have better performance using the adaptive step-size than using the constant step-size. The optimal step-size has similar reconstruction result like the adaptive one and thus not presented. For ACoSaMP and ASP we observe that it is better to use a = 1 instead of a = 2ℓ-p ℓ . Compared to each other we see that ACoSaMP and ASP achieve better recovery than AHTP and AIHT. Between the last two, AHTP is better. Though AIHT has inferior behavior, we should mention that with regards to running time AIHT is the most efficient. Afterwards we have AHTP and then ACoSaMP and ASP. Compared to ℓ 1 and GAP we observe that ACoSaMP and ASP have competitive results.

With the above observations, we turn to test operators with higher redundancy and see the effect of linear dependencies in them. We test two operators. The first is a random tight frame as before but with redundancy factor of 2. The second is the two dimensional finite difference operator Ω 2D-DIF . In Fig. 3 we present the phase diagrams for both operators using AIHT with an adaptive changing step-size, AHTP with an adaptive changing step-size, ACoSaMP with a = 1, and ASP with a = 1. As observed before, also in this case the ACoSaMP and ASP outperform AIHT and AHTP in both cases and AHTP outperform AIHT. We mention again that the better performance comes at the cost of higher complexity. In addition, as we expected, having redundancies in Ω results with a better recovery.

Reconstruction of High Dimensional Images in the Noisy Case

We turn now to test the methods for high dimensional signals. We use RASP and RACoSaMP (relaxed versions of ASP and ACoSaMP defined in Section 5.4) for the reconstruction of the Shepp-Logan phantom from few number of measurements. The sampling operator is a two dimensional Fourier transform that measures only a certain number of radial lines from the Fourier transform. The cosparse operator is Ω 2D-DIF and the cosparsity used is the actual cosparsity of the signal under this operator (ℓ = 128014). The phantom image is presented in Fig. 4(a). Using the RACoSaMP and RASP we get a perfect reconstruction using only 15 radial lines, i.e., only m = 3782 measurements out of d = 65536 which is less then 6 percent of the data in the original image. The algorithms requires less than 20 iterations for having this perfect recovery. For AIHT and RAHTP we achieve a reconstruction which is only close to Exploring the noisy case, we perform a reconstruction using RASP of a noisy measurement of the phantom with 22 radial lines and signal to noise ratio (SNR) of 20. Figure 4(c) presents the noisy image, the result of applying inverse Fourier transform on the measurements, and Fig. 4(d) presents its reconstruction result. Note that for the minimization process we solve conjugate gradients, in each iteration and take only the real part of the result and crop the values of the resulted image to be in the range of [0, 1]. We get a peak SNR (PSNR) of 36dB. We get similar results using RACoSaMP but using more radial lines (25).

Discussion and Conclusion

In this work we presented new pursuits for the cosparse analysis model. A theoretical study of these algorithms was performed giving guarantees for stable recovery under the assumptions of the Ω-RIP and the existence of an optimal or a near optimal projection. We showed that optimal projections exists for some non-trivial operators, i.e., operators that do not take us back to the synthesis case. In addition, we showed experimentally that using simpler kind of projections is possible in order to get good reconstruction results. We demonstrated both in the theoretical and the empirical results that linear dependencies within the analysis dictionary are favorable and enhance the recovery performance. We are aware that there are still some open questions in this work and we leave them for future research. This should deal with following:

• Our work assumed the existence of a procedure that finds a cosupport that implies a near optimal projection with a constant C ℓ . Two examples for optimal cosupport slection schemes were given. However, the existence of an optimal or a near optimal scheme for a general operator is still an open question. The question is: for which types of Ω and values of C ℓ we can find an efficient procedure that implies a near optimal projection.

• As we have seen in the simulations, the thresholding procedure, though not near optimal with the theorems required constants, provides good reconstruction results. A theoretical study of the analysis greedy-like techniques with this cosupport selection scheme is required.

• A family of analysis dictionaries that deserves a special attention is the family of tight frame operators. In synthesis, there is a parallel between the guarantees of ℓ 1 -synthesis and the greedy like algorithms. The fact that a guarantee with a tight frame Ω exists for ℓ 1 -analysis encourage us to believe that similar guarantees exist also for the analysis greedy-like techniques.

• In this paper, the noise e was considered to be adversarial. Random white Gaussian case was considered for the synthesis case in [START_REF] Giryes | RIP-based near-oracle performance guarantees for SP, CoSaMP, and IHT[END_REF] resulting with near-oracle performance guarantees. It would be interesting to verify whether this is also the case for the analysis framework.

Appendix A. Proofs of Theorem 3.7 and Theorem 3.8

Theorem 3.7 (Theorem 3.3 in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]): Let M ∈ R m×d be a random matrix that satisfies that for any z ∈ R d and 0

< ǫ ≤ 1 3 P Mz 2 2 -z 2 2 ≥ ǫ z 2 2 ≤ e -C M mǫ 2 ,
where C M > 0 is a constant. For any value of ǫ r > 0, if

m ≥ 32 C M ǫ 2 r log( L corank r ) + (d -r) log(9/ǫ r ) + t ,
then δ corank r ≤ ǫ r with probability exceeding 1e -t . Theorem 3.8: Under the same setup of Theorem 3.7, for any ǫ ) d-r . We could stop here and use directly Theorem 2.1 from [START_REF] Mendelson | Uniform uncertainty principle for bernoulli and subgaussian ensembles[END_REF] to get the desired result for Theorem 3.7. However, we present the remaining of the proof using a proof technique from [START_REF] Baraniuk | A simple proof of the restricted isometry property for random matrices[END_REF][START_REF] Rauhut | Compressed sensing and redundant dictionaries[END_REF]. Using union bound and the properties of M we have that with probability exceeding 1 -L corank

ǫ ℓ > 0 if m ≥ 32 C M ǫ 2 ℓ (p -ℓ) log 9p (p -ℓ)ǫ ℓ + t , then δ ℓ ≤ ǫ ℓ with probability exceeding 1 -e -t . Proof: Let ǫ = ǫ r /4, B d-r = {z ∈ R d-r , z 2 ≤ 1}
≤ √ 1 + ǫ + M v 0 -ṽ v 0 -ṽ 2 2 v 0 -ṽ 2 ≤ √ 1 + ǫ + 1 + δ corank r ǫ.
The first inequality is due to the triangle inequality; the second one follows from (A.1) and arithmetics; and the last inequality follows from the definition of δ corank r , the properties of ǫ-net and the fact that v 0 -ṽ v 0 -ṽ 2 2

= 1. Reordering (A.2) gives

1 + δ corank r ≤ 1 + ǫ (1 -ǫ) 2 ≤ 1 + 4ǫ = 1 + ǫ r . (A.3)
where the inequality holds because ǫ r ≤ 0.5 and ǫ = ǫ r 4 ≤ 1 8 . Since we want (A.3) to hold with probability greater than 1e -t it remains to require L corank r (1 + 8 ǫ r ) d-r e -C M mǫ 2 r 32 ≤ e -t . Using the fact that (1 + 8 ǫ r ) ≥ 9 ǫ r and some arithmetics we get [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] and this completes the proof of the theorem. , where for the second inequality we use the fact that δ 4ℓ-3p ≤ 1 combined with the inequality of Lemma 6.10, and for the last inequality we use the fact that (1

+ (1 + δ 4ℓ-3p ) √ C ℓ +C ℓ ) 2 ≤ (1 + 2δ 4ℓ-3p √ C ℓ +C ℓ )( 2+C ℓ 1+C ℓ + 2 √ C ℓ +C ℓ
) together with a few algebraic steps. Taking square-root on both sides of (E.4) provides the desired result.

Appendix F. Proof of Lemma 6.12 Lemma 6.12: Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . if

C 2ℓ-p < σ 2 M (1 + γ) 2 σ 2 M (1 + γ) 2 -1
, then there exists δ ACoSaMP (C 2ℓ-p , σ 2 M , γ) > 0 such that for any δ 2ℓ-p < δ ACoSaMP (C 2ℓ-p , σ 2 M , γ)

P Λt (x -w) 2 ≤ η 2 e 2 + ρ 2 x -xt-1 2 .
The constants η 2 and ρ 2 are as defined in Theorem 6.8.

In the proof of the lemma we use the following Proposition. Proposition E.1: For any two given vectors x 1 , x 2 and any constant c > 0 it holds that

x 1 + x 2 2 2 ≤ (1 + c) x 1 2 2 + 1 + 1 c x 2 2 . (F.1)
The proof of the proposition is immediate using the inequality of arithmetic and geometric means. We turn to the proof of the lemma.

Proof: Looking at the step of finding new cosupport elements one can observe that Q Λ ∆ is a near optimal projection for M * y t-1 resid = M * (y -Mx t-1 ) with a constant C 2ℓ-p . The fact that Λt-1 ∩ Λ ≥ 2ℓp combined with [START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF] gives

(I -Q Λ ∆ )M * (y -Mx t-1 ) 2 2 ≤ C 2ℓ-p (I -Q Λt-1 ∩Λ )M * (y -Mx t-1 ) 2 2 .
Using simple projection properties and the fact that Λt ⊆ Λ ∆ with z = M * (y -Mx t-1 ) we have

Q Λt z 2 2 ≥ Q Λ ∆ z 2 2 = z 2 2 -(I -Q Λ ∆ )z 2 2 ≥ z 2 2 -C 2ℓ-p (I -Q Λt-1 ∩Λ )z 2 2 (F.2) = z 2 2 -C 2ℓ-p z 2 2 -Q Λt-1 ∩Λ z 2 2 = C 2ℓ-p Q Λt-1 ∩Λ z 2 2 -(C 2ℓ-p -1) z 2 2 .
We turn to bound the LHS of (F.2) from above. Noticing that y = Mx + e and using (F.1) with a constant γ 1 > 0 gives Putting (F.4) into (F.3) and using [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] and Corollary 3.3 gives (F.5)

Q Λt M * (y -Mx t-1 ) 2 2 ≤ 1 + 1 γ 1 Q Λt M * e
Q Λt M * (y -Mx t-
+ 1 + α + δ 4ℓ-3p + δ 4ℓ-3p α (1 + γ 1 ) x -xt-1 2 2 .
We continue with bounding the RHS of (F.2) from below. For the first element of the RHS we use an altered version of (F.1) with a constant γ > 0 and have (F.9)

Q Λt-1 ∩Λ M * (y -Mx t-1 ) 2 2 ≥ 1 1 + γ 2 Q Λt-1 ∩Λ M * M(x -xt-1 ) 2 2 - 1 γ 2 Q Λt-1 ∩Λ M *
≤ (1 + γ 3 )(1 + δ 2ℓ-p )σ 2 M (x -xt-1 ) 2 2 + 1 + 1 γ 3 σ 2 M e 2 2 .
By combining (F.5), (F.8) and (F.9) with (F.

- C 2ℓ-p (1 + α)(1 + γ 1 )(1 + γ 2 ) 1 1 + β - δ 2ℓ-p β x -xt-1 2 2 .
The smaller the coefficient of xxt-1 2 2 , the better convergence guarantee we obtain. Thus, we choose β = so that the coefficient is minimized. The values of γ 1 , γ 2 , γ 3 provide a tradeoff between the convergence rate and the size of the noise coefficient. For smaller values we get better convergence rate but higher amplification of the noise. We make no optimization on their values and choose them to be γ 1 = γ 2 = γ 3 = γ for an appropriate γ > 0. Thus, the above yields

P Λt (x -xt-1 ) 2 2 ≤ 1 + δ 3ℓ-2p γ(1 + α) + (1 + δ 2ℓ-p )C 2ℓ-p γ(1 + α)(1 + γ) + (C 2ℓ-p -1)(1 + γ)σ 2 M (1 + α)(1 + γ)γ e 2

  Ω 1D-DIF z.

  Optimal projection. The projection ℓ 2 -norm error is 2.5.

Algorithm 2

 2 Subspace Pursuit (SP) and CoSaMP Require: k, M, y where y = Mx + e, k is the cardinality of x and e is an additive noise. a = 1 (SP), a = 2 (CoSaMP). Ensure: xCoSaMP or xSP : k-sparse approximation of x.

Figure 2 :

 2 Figure 2: Recovery rate for a random tight frame with p = 144 and d = 120. From left to right, up to bottom: AIHT with a constant step-size, AIHT with an adaptive changing step-size, AHTP with a constant step-size, AHTP with an adaptive changing step-size, ACoSaMP with a = 2ℓ-p ℓ , ACoSaMP with a = 1, ASP with a = 2ℓ-p ℓ , ASP with a = 1, A-ℓ 1 -minimization and GAP.

Figure 3 :Figure 4 :

 34 Figure 3: Recovery rate for a random tight frame with p = 240 and d = 120 (up) and a finite difference operator (bottom). From left to right: AIHT and AHTP with an adaptive changing step-size, and ACoSaMP and ASP with a = 1.

2 ǫ

 2 and Ψ an ǫ-net for B d-r with size |Ψ| ≤ 1 + any subspaceW B Λ = W Λ ∩ B d-r such that Λ ∈ L corank r we can build an orthogonal matrix U Ω ∈ R d×(d-r) such that W B Λ = {U Λ z, z ∈ R d-r , z 2 ≤ 1} = U Λ B d-r . It is easy to see that Ψ Λ = U Λ Ψ d-r is an ǫ-net for W B Λ and that Ψ A corank r = ∪ Λ∈L corank r Ψ Λ is an ǫ-net for A corankr ∩ B d , where Ψ A corank r

r( 1 2 every v ∈ Ψ A corank r satisfies ( 1 -ǫ) v 2 2 ≤ Mv 2 2 ≤ ( 1 + ǫ) v 2 2 . (A. 1 )= Mv 0 2 ≤

 121221212 + 2 ǫ ) d-r e -C M mǫ 2According to the definition of δ corank r it holds that 1+ δ corank r = sup v∈A corank r ∩B d Mv 2 . Since A corank r ∩ B d is a compact set there exists v 0 ∈ A corank r ∩ B d thatachieves the supremum. Denoting by ṽ its closest vector in Ψ A corank r and using the definition of Ψ A corank r we have v 0ṽ 2 ≤ ǫ. This yields 1 + δ corank r Mṽ 2 + M(v 0ṽ) 2 (A.2)

2 2 + 2 2≤ 2 2 + 1 + 1 α 2 (F. 4 )≤ ( 1 2 2 + 1 + 1 α

 22211241211 (1 + γ 1 ) Q Λt M * M(xxt-1 ) .1) again, now with a constant α > 0, we haveQ Λt M * M(xxt-1 ) (1 + α) Q Λt (xxt-1 ) Q Λt (I -M * M)(xxt-1 ) 2 + α) xxt-1 2 2 -(1 + α) P Λt (xxt-1 ) Q Λt (I -M * M)(xxt-

  e

6 ) 1 2 2≤ 2 2 + 1 + 1 γ 3 M * e 2 2

 6122132 Using the altered form again, for the first element in the RHS of (F.6), with a constant β > 0 givesQ Λt-1 ∩Λ M * M(xxt-1 ) ∩Λ (M * M -I)(xxt-1 ) .7) in (F.6) and using the RIP properties and (18) provide Q Λt-1 ∩Λ M * (y -Mx t-1 ) .1), with a constant γ 3 > 0,[START_REF] Pati | Orthonormal matching pursuit : recursive function approximation with applications to wavelet decomposition[END_REF], and some basic algebraic steps we have for the second element in the RHS of (F.2)M * (y -Mx t-1 ) (1 + γ 3 ) M * M(xxt-1 )

  1 resid , ak). Update the support: T t = T t-1 ∪ T ∆ . Compute a temporary representation: w = M † T t y. Prune small entries: T t = supp(w, k). Calculate a new representation: xt CoSaMP = w T t for CoSaMP, and xt SP = M † T t y for SP. Update the residual: y t resid = y -Mx t CoSaMP for CoSaMP, and y t resid = y -Mx t SP for SP. end while Form the final solution xCoSaMP = xt CoSaMP for CoSaMP and xSP = xt SP for SP.

  1 ) Find a new cosupport: Λt = Ŝℓ (x g ) Calculate a new estimate: xt AIHT = Q Λt x g for AIHT, and xt AHTP = argmin x y -Mx 2 2 s.t. Ω Λt x = 0 for AHTP. end while Form the final solution xAIHT = xt AIHT for AIHT and xAHTP = xt AHTP for AHTP. Algorithm 4 Analysis Subspace Pursuit (ASP) and Analysis CoSaMP (ACoSaMP) Require: ℓ, M, Ω, y, a where y = Mx + e, ℓ is the cosparsity of x under Ω and e is the additive noise.

Ensure: xACoSaMP or xASP : ℓ-cosparse approximation of x. Initialize the cosupport Λ 0 = {i, 1 ≤ i ≤ p}, the residual y 0 resid = y and set t = 0. while halting criterion is not satisfied do t = t + 1.

Find new cosupport elements:

Λ ∆ = Ŝaℓ (M * y t-1 resid ). Update the cosupport: Λt = Λt-1 ∩ Λ ∆ . Compute a

temporary estimate: w = argmin x y -Mx 2 2 s.t. Ω Λt x = 0. Enlarge the cosupport: Λt = Ŝℓ (w). Calculate a new estimate: xt ACoSaMP = Q Λt w for ACoSaMP, and xt ASP = argmin x y -Mx 2 2 s.t. Ω Λt x = 0 for ASP. Update the residual: y t resid = y -Mx t ACoSaMP for ACoSaMP, and y t resid = y -Mx t ASP for ASP. end while Form the final solution xACoSaMP = xt ACoSaMP for ACoSaMP and xASP = xt ASP for ASP.

  By substituting (E.2) and (E.3) into (E.1) we havexxt 2 2 ≤ (1 + C ℓ ) xw 2 2 + 2δ 4ℓ-3p C ℓ xw 2 2 + 2 1 + δ 3ℓ-2p C ℓ e 2 xw 2 (E.4) ≤ (1 + 2δ 4ℓ-3p C ℓ + C ℓ ) xw 2 + 2 (1 + δ 3ℓ-2p )C ℓ e 2 xw 2

	≤ + ≤	1 + 2δ 4ℓ-3p 1 -δ 2 √ C ℓ + C ℓ 4ℓ-3p 2 1 + δ 3ℓ-2p (1 + (1 + δ 4ℓ-3p ) P Λt (x -w) √ C ℓ + C ℓ ) 2 2 4ℓ-3p (1 -δ 4ℓ-3p ) 1 -δ 2 1 + 2δ 4ℓ-3p √ C ℓ + C ℓ 1 -δ 2 4ℓ-3p P Λt (x -w) 2 +	P Λt (x -w) 2 e 2 + 2+C ℓ 1 -δ 4ℓ-3p (1 + δ 3ℓ-2p )(1 + 2 (1 -δ 4ℓ-3p ) 2 √ C ℓ + C ℓ ) e 2 1+C ℓ + 2 √ C ℓ + C ℓ 1 + δ 3ℓ-2p 2	e 2 2

  2) we have(1 + α)(1 + γ 1 ) P Λt (xxt-1 )Dividing both sides by (1 + α)(1 + γ 1 ) and gathering coefficients give

	P Λt (x -xt-1 )	2 2 ≤	1 + δ 3ℓ-2p γ 1 (1 + α)	+	(1 + δ 2ℓ-p )C 2ℓ-p γ 2 (1 + α)(1 + γ 1 )	+	(C 2ℓ-p -1)(1 + γ 3 )σ 2 M (1 + α)(1 + γ 1 )γ 3	e 2 2	(F.11)
			+ 1 +	δ 4ℓ-3p α		+	M (C 2ℓ-p -1)(1 + γ 3 )(1 + δ 2ℓ-p )σ 2 (1 + α)(1 + γ 1 )
					2 2 ≤	(1 + γ 1 )(1 + δ 3ℓ-2p ) γ 1	e 2 2 + C 2ℓ-p	(1 + δ 2ℓ-p ) γ 2	e 2 2	(F.10)
	+(C 2ℓ-p -1) 1 +	1 γ 3	σ 2 M e 2 2 + 1 + α + δ 4ℓ-3p +	δ 4ℓ-3p α	(1 + γ 1 ) x -xt-1 2 2
	+(C 2ℓ-p -1)(1 + γ 3 )(1 + δ 2ℓ-p )σ 2 M (x -xt-1 )	2 2 -C 2ℓ-p	1 1 + β	-	δ 2ℓ-p β	1 1 + γ 2	x -xt-1 2 2 .

By the abuse of notation we use the same notation for the selection sub-matrices of rows and columns. The selection will be clear from the context since in analysis the focus is always on the rows and in synthesis on the columns.

At a first glance one would think that the conditions should be in terms of δ corank 2r-d and δ corank 4r-3d . However, given two cosparse vectors with coranks r 1 and r 2 the best estimation we can have for the corank of their sum is r 1 + r 2p.

Remark that we even improve the condition of the idealized case in[START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF] to be δ 2ℓ-p ≤ 1 3 instead of δ 2ℓ-p ≤ 1 5 .

A matlab package with code for the experiments performed in this paper is in preparation for an open source distribution.

Acknowledgment

The authors would like to thank Jalal Fadili for fruitful discussion, and the unknown reviewers for the important remarks that helped to improved the shape of the paper. Without both of them, the examples of the optimal projections would not have appeared in the paper. This research was supported by New York Metropolitan Research Fund. R. Giryes is grateful to the Azrieli Foundation for the award of an Azrieli Fellowship. This work was supported in part by the EU FP7, SMALL project under FET-Open grant number 225913, and EPSRC grants EP/J015180/1 and EP/F039697/1. R. Gribonval acknowledges the support of the European Research Council, PLEASE project, under grant ERC-StG-2011-277906. MED acknowledges support of his position from the Scottish Funding Council and their support of the Joint Research Institute with the Heriot-Watt University as a component part of the Edinburgh Research Partnership.

We turn now to the proof of Theorem 3.8. Its proof is almost identical to the previous proof but with the difference that instead of r, L corank r and δ corank r we look at ℓ, L ℓ and δ ℓ . In this case we do not know what is the dimension of the subspace that each cosupport implies. However, we can have a lower bound on it using pℓ. Therefore, we use B p-ℓ instead of B d-r . This change provides us with a condition similar to [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] but with pℓ in the second coefficient instead of dr. By using some arithmetics, noticing that the size of L ℓ is p ℓ and using Stirling's formula for upper bounding it we get [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] and this completes the proof.

Appendix B. Proof of Lemma 6.6 Lemma 6.6: Consider the problem P and apply either AIHT or AHTP with a constant step size µ satisfying 1 µ ≥ 1 + δ 2ℓ-p or an optimal step size. Then, at the t-th iteration, the following holds:

For the optimal step size the bound is achieved with the value µ = 1 1+δ 2ℓ-p . Proof: We consider the AIHT algorithm first. We take similar steps to those taken in the proof of Lemma 3 in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]. Since 1 µ ≥ 1 + δ 2ℓ-p , we have, from the Ω-RIP property of M,

Thus,

Note that by definition, xt = Q Ŝℓ xt-1 + µM * (y -Mx t-1 ) . Hence, by the C ℓ -near optimality of the projection, we get

Putting this into (B.2), we obtain the desired result for the AIHT algorithm.

We can check that the same holds true for the AHTP algorithm as follows: suppose that xt-1 AHTP is the (t -1)-st estimate from the AHTP algorithm. If we now initialize the AIHT algorithm with this estimate and obtain the next estimate xt ÃIHT , then the inequality of the lemma holds true with xt ÃIHT and xt-1 AHTP in place of xt and xt-1 respectively. On the other hand, from the algorithm description, we know that the t-th estimate xt AHTP of the AHTP satisfies

This means that the result holds for the AHTP algorithm as well.

Using a similar argument for the optimal changing step size we note that it selects the cosupport that minimizes Mx -Mx t 2 2 . Thus, for AIHT and AHTP we have that Mx

for any value of µ, where xt

Opt and xt µ are the recovery results of AIHT or AHTP with an optimal changing step-size and a constant step-size µ respectively. This yields that any theoretical result for a constant step-size selection with a constant µ holds true also to the optimal changing-step size selection. In particular this is true also for µ = 1 1+δ 2ℓ-p . This choice is justified in the proof of Lemma 6.7.

Appendix C. Proof of Lemma 6.7 Lemma 6.7: Suppose that the same conditions of Theorem 6.5 hold true. If

where

Proof: First, suppose that y -Mx t-1 2 2 > η 2 e 2 2 . From Lemma 6.6, we have

Remark that all the coefficients in the above are positive because 1

and, by the triangle inequality,

Therefore, from (C.1), y -Mx t 2 2 < c 4 y -Mx t-1 2 2 . This is the second part of the lemma. Now, suppose that y -Mx t-1 2 2 ≤ η 2 e 2 2 . This time we have

Applying this to (C.1), we obtain

Thus, the proof is complete as soon as we show c 4 < 1, or c 4 -1 < 0.

To see c 4 -1 < 0, we first note that it is equivalent to-all the subscripts are dropped from here on for simplicity of notation-

Solving this quadratic equation in 1 µ , we want

Such µ exists only when b 2 b 2 1 < 1. Furthermore, we have already assumed 1+δ ≤ 1 µ and we know (1-δ) b 1b 2 1b 2 < 1 + δ, and hence the condition we require is

which is what we desired to prove.

As we have seen in Lemma 6.6, for changing optimal step-size selection, (49) holds for any value of µ that satisfies the above conditions. Thus, in the bound of changing optimal step-size we put a value of µ that minimizes c 4 . This minimization result with 1 µ = √ b 2 (1δ 2ℓ-p ). However, since we need 1 µ ≥ 1 + δ 2ℓ-p and have that √ b 2 (1δ 2ℓ-p ) < b 1 (1δ 2ℓ-p ) < 1 + δ 2ℓ-p we set 1 µ = 1 + δ 2ℓ-p in c 4 for the bound in optimal changing step-size case.

Appendix D. Proof of Lemma 6.10 Lemma 6.10: Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . For each iteration we have

Proof: Since w is the minimizer of y -Mv where u is a vector satisfying Ω Λt u = 0. Turning to look at Q Λt (xw) 2 2 and using (D.2) with u = Q Λt (xw), we have

where the first inequality follows from the Cauchy-Schwartz inequality, the projection property that Q Λt = Q Λt Q Λt and the fact that xw = Q Λ∩ Λt (xw). The last inequality is due to the Ω-RIP properties, Corollary 3.6 and that according to Table 1

Utilizing the last inequality with the fact that xw 2 2 = P Λt (xw)

xw 2 2 ≤ P Λt (xw)

By moving all terms to the LHS we get a quadratic function of xw 2 . Thus, xw 2 is bounded from above by the larger root of that function; this with a few simple algebraic steps gives the inequality in (58).

Appendix E. Proof of Lemma 6.11 Lemma 6.11: Consider the problem P and apply ACoSaMP with a = 2ℓ-p ℓ . For each iteration we have

where η 1 and ρ 1 are the same constants as in Theorem 6.8.

Proof: We start with the following observation

and turn to bound the second and last terms in the RHS. For the second term, using the fact that xt = Q Ŝℓ (w) w with [START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF] gives

For bounding the last term, we look at its absolute value and use (D

By using the triangle and Cauchy-Schwartz inequalities with the fact that x-w = Q Λ∩ Λt (x-w) and w

where the last inequality is due to the Ω-RIP definition and Corollary 3.

Since P Λt w = P Λt xt-1 = 0 the above inequality holds also for P Λt (xxt-1 ) 2 2 . Inequality (61) follows since the right-hand side of (F.12) is smaller than the square of the right-hand side of (61).

Before ending the proof, we notice that ρ 2 , the coefficient of xxt-1 2 2 is defined only when which is equivalent to (60). By moving the terms in the RHS to the LHS we get a quadratic function of δ 2ℓ-p . The condition in (60) guarantees that its constant term is smaller than zero and thus there exists a positive δ 2ℓ-p for which the function is smaller than zero. Therefore, for any δ 2ℓ-p < δ ACoSaMP (C 2ℓ-p , σ 2 M , γ) (F.13) holds, where δ ACoSaMP (C 2ℓ-p , σ 2 M , γ) > 0 is the square of the positive solution of the quadratic function.