Keywords: Synthesis, Analysis, Sparse Representations, Union of Subspaces, Pursuit Algorithms, Greedy Algorithms, Compressed-Sensing. 6

After a decade of extensive study of the sparse representation synthesis model, we can safely say that this is a mature and stable field, with clear theoretical foundations, and appealing applications. Alongside this approach, there is an analysis counterpart model, which, despite its similarity to the synthesis alternative, is markedly different. Surprisingly, the analysis model did not get a similar attention, and its understanding today is shallow and partial.

In this paper we take a closer look at the analysis approach, better define it as a generative model for signals, and contrast it with the synthesis one. This work proposes effective pursuit methods that aim to solve inverse problems regularized with the analysis-model prior, accompanied by a preliminary theoretical study of their performance. We demonstrate the effectiveness of the analysis model in several experiments, and provide a detailed study of the model associated with the 2D finite difference analysis operator, a close cousin of the TV norm.

Introduction

Situated at the heart of signal and image processing, data models are fundamental for stabilizing the solution of inverse problems, and enabling various other tasks, such as compression, detection, separation, sampling, and more. What are those models? Essentially, a model poses a set of mathematical properties that the data is believed to satisfy. Choosing these properties (i.e. the model) carefully and wisely may lead to a highly effective treatment of the signals in question and consequently to successful applications.

Throughout the years, a long series of models has been proposed and used, exhibiting an evolution of ideas and improvements. In this context, the past decade has been certainly the era of sparse and redundant representations, a novel synthesis model for describing signals [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF][START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF][START_REF] Starck | Sparse Image and Signal Processing -Wavelets, Curvelets, Morphological Diversity[END_REF]. Here is a brief description of this model:

Assume that we are to model the signal x ∈ R d . The sparse and redundant synthesis model suggests that this signal could be described as x = Dz, where D ∈ R d×n is a possibly redundant dictionary (n ≥ d), and z ∈ R n , the signal's representation, is assumed to be sparse. Measuring the cardinality of non-zeros of z using the 'ℓ 0 -norm', such that z 0 is the count of the non-zeros in z, we expect z 0 to be much smaller than n. Thus, the model essentially assumes that any signal from the family of interest could be described as a linear combination of few columns from the dictionary D. The name "synthesis" comes from the relation x = Dz, with the obvious interpretation that the model describes a way to synthesize a signal.

This model has been the focus of many papers, studying its core theoretical properties by exploring practical numerical algorithms for using it in practice (e.g. [START_REF] Candès | The Dantzig selector: statistical estimation when p is much larger than n[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF][START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]), evaluating theoretically these algorithms' performance guarantees (e.g. [START_REF] Ben-Haim | Coherence-based performance guarantees for estimating a sparse vector under random noise[END_REF][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][START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF][START_REF] Tropp | Just relax: Convex programming methods for subset selection and sparse approximation[END_REF]), addressing ways to obtain the dictionary from a bulk of data (e.g. [START_REF] Aharon | An algorithm for designing of overcomplete dictionaries for sparse representation[END_REF][START_REF] Engan | Multi-frame compression: Theory and design[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF][START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF]), and beyond all these, attacking a long series of applications in signal and image processing with this model, demonstrating often state-of-the-art results (e.g. [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (mca)[END_REF][START_REF] Llagostera Casanovas | Blind audiovisual source separation based on sparse representations[END_REF][START_REF] Plumbley | Sparse representations in audio and music: from coding to source separation[END_REF]). Today, after a decade of an extensive study along the above lines, with nearly 4000 papers 1 written on this model and related issues, we can safely say that this is a mature and stable field, with clear theoretical foundations, and appealing applications.

Interestingly, the synthesis model has a "twin" that takes an analysis point of view. This alternative assumes that for a signal of interest, the analyzed vector Ωx is expected to be sparse, where Ω ∈ R p×d is a possibly redundant analysis operator (p ≥ d). Thus, we consider a signal as belonging to the analysis model if Ωx 0 is small enough. Common examples of analysis operators include: the shift invariant wavelet transform Ω WT [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF]; the finite difference operator Ω DIF , which concatenates the horizontal and vertical derivatives of an image and is closely connected to total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]; the curvelet transform [START_REF] Starck | The curvelet transform for image denoising[END_REF], and more. Empirically, analysis models have been successfully used for a variety of signal processing tasks [START_REF] Farsiu | Advances and challenges in super-resolution[END_REF][START_REF] Kluckner | Exploiting redundancy for aerial image fusion using convex optimization[END_REF][START_REF] Starck | The curvelet transform for image denoising[END_REF][START_REF] Starck | Redundant multiscale transforms and their application for morphological component analysis[END_REF][START_REF] Starck | Gray and color image contrast enhancement by the curvelet transform[END_REF] such as denoising, deblurring, and most recently compressed sensing, but this has been done with little theoretical justification.

It is well known by now [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] that for a square and invertible dictionary, the synthesis and the analysis models are the same with D = Ω -1 . The models remain similar for more general dictionaries, although then the gap between them is unexplored. Despite the close-proximity between the two-synthesis and analysis-models, the first has been studied extensively while the second has been left aside almost untouched. In this paper we aim to bring justice to the analysis model by addressing the following set of topics:

1. Cosparsity: In Section 2 we start our discussion with a closer look at the sparse analysis model in order to better define it as a generative model for signals. We show that, while the synthesis model puts an emphasis on the non-zeros of the representation vector z, the analysis model draws its strength from the zeros in the analysis vector Ωx. 2. Union of Subspaces: Section 2 is also devoted to a comparison between the synthesis model and the analysis one. We know that the synthesis model described above is an instance of a wider family of models, built as a finite union of subspaces [START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF]. By choosing all the sub-groups of columns from D that could be combined linearly to generate signals, we get an exponentially large family of low-dimensional subspaces that cover the signals of interest. Adopting this perspective, the analysis model can obtain a similar interpretation. How are the two related to each other? Section 2 considers this question and proposes a few answers. 3. Uniqueness: We know that the spark of the dictionary governs the uniqueness properties of sparse solutions of the underdetermined linear system Dz = x [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF]. Can we derive a similar relation for the analysis case? As a platform for studying the analysis uniqueness properties, we consider an inverse problem of the form y = Mx, where M ∈ R m×d and m < d, and y ∈ R m is a measurement vector. Put roughly (and this will be better defined later on), assuming that x comes from the sparse analysis model, could we claim that there is only one possible solution x that can explain the measurement vector y? Section 3 presents this uniqueness study. 4. Uniqueness: Worked Examples. Based on the study of the analysis uniqueness properties, we characterize the required number of measurements for the uniqueness of the signal that satisfies the analysis model in the case of analysis operator Ω in general position and the 2D one-step finite difference operator Ω DIF . 5. Pursuit Algorithms: Armed with a deeper understanding of the analysis model, we may ask how to efficiently find x for the above-described linear inverse problem. As in the synthesis case, we can consider either relaxation-based methods or greedy ones. In Section 4 we present two numerical approximation algorithms: a greedy algorithm termed "Greedy Analysis Pursuit" (GAP) that resembles the Orthogonal Matching Pursuit (OMP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]-adapted to the analysis model-, and the previously considered ℓ 1 -minimization approach [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF][START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Selesnick | Signal restoration with overcomplete wavelet transforms: Comparison of analysis and synthesis priors[END_REF]. Section 5 accompanies the presentation of GAP with a theoretical study of its performance guarantee, deriving a condition that resembles the ERC obtained for OMP [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. Similarly, we study the terms of success of the ℓ 1 -minimization approach for the analysis model, deriving a condition that is similar to the one obtained for the synthesis sparse model [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. 6. Tests: In Section 6 we demonstrate the effectiveness of the analysis model and the pursuit algorithms proposed in several experiments, starting from synthetic ones (involving random analysis operators) and going all the way to a compressed-sensing test for an image based on the analysis model: the Shepp Logan phantom.

We believe that with the above set of contributions, the cosparse analysis model becomes a well-defined and competitive model to the synthesis counterpart, equipped with all the necessary ingredients for its practical use. Furthermore, this work leads to a series of new questions that are parallel to those studied for the synthesis model-developing novel pursuit methods, a theoretical study of pursuit algorithms for handling other inverse problems, training Ω just as done for D, and more. We discuss these and other topics in Section 7.

Related Work. Several works exist in the literature that are related to the analysis model. The work by Elad et. al. [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] was the first to observe the dichotomy of analysis and synthesis models for signals. Their study, done in the context of the Maximum-A-Posteriori Probability estimation, presented the two alternatives and explored cases of equivalence between the two. They demonstrated a superiority of the analysis-based approach in signal denoising. Further empirical evidence of the effectiveness of the analysis-based approach for signal and image restroation can be found in [START_REF] Portilla | Image restoration through l0 analysis-based sparse optimization in tight frames[END_REF] and [START_REF] Selesnick | Signal restoration with overcomplete wavelet transforms: Comparison of analysis and synthesis priors[END_REF]. In [START_REF] Selesnick | Signal restoration with overcomplete wavelet transforms: Comparison of analysis and synthesis priors[END_REF] it was noted that the nonzero coefficients play a different role in the analysis and synthesis forms but the importance of the zero coefficients for the analysis model-which is reminiscent of signal characterizations through the zero-crossings of their undecimated wavelet transform [START_REF] Mallat | Zero-crossings of a wavelet transform[END_REF]-was not explicitly identified. More recently, Candès et al. [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] provided a theoretical study on the error when the analysis-based ℓ 1 -minimization is used in the context of compressed sensing. Our work is closely related to these contributions in various ways, and we shall return to these papers when diving into the details of our study.

A Closer Look at the Cosparse Analysis Model

We start our discussion with the introduction of the sparse analysis model, and the notion of cosparsity that is fundamental for its definition. We also describe how to interpret the analysis model as a generative one (just like the synthesis counterpart). Finally, we consider the interpretation of the sparse analysis and synthesis models as two manifestations of union-of-subspaces models, and show how they are related.

Introducing Cosparsity

As described in the introduction, a conceptually simple model for data would be to assume that each signal we consider can be expressed (i.e., wellapproximated) as a combination of a few building atoms. Once we take this view, a simple synthesis model can be thought of: First, there is a collection of the atomic signals {d j } n j=1 ∈ R d that we concatenate as the columns of a dictionary, denoted by D ∈ R d×n . Here, typically n ≥ d, implying that the dictionary is redundant. Second, the signal x ∈ R d can be expressed as a linear combination of some atoms of D, thus there exists z ∈ R n such that x = Dz. Third and most importantly, x must lie in a low dimensional subspace, and in order to ensure this, very few atoms are used in the expression x = Dz, i.e., the number of non-zeros z 0 is very small. By the observation that z 0 is small, we say that x has a sparse representation in D. The number k = z 0 is the sparsity of the coefficient vector z and, by extension, of the signal x.

Often, the validity of the above described sparse synthesis model is demonstrated by applying a linear transform to a class of signals to be processed and observing that most of the coefficients are close to zero, exhibiting sparsity. In signal and image processing, discrete transforms such as wavelet, Gabor, curvelet, contourlet, shearlet, and others [START_REF] Do | The contourlet transform: An efficient directional multiresolution image representation[END_REF][START_REF] Labate | Sparse multidimensional representation using shearlets[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF][START_REF] Starck | The curvelet transform for image denoising[END_REF], are of interest, and this empirical observation seems to give a good support for the sparse synthesis model. Indeed, when aiming to claim optimality of a given transform, this is exactly the approach taken -show that for a (theoretically-modeled) class of signals of interest, the transform coefficients tend to exhibit a strong decay. However, one cannot help but noticing that this approach of validating the synthesis model seems to actually validate another 'similar' model; we are considering a model where the signals of interest have sparse analysis representations. This point is especially pronounced when the transform used is over-complete or redundant.

Let us now look more carefully at the above mentioned model that seems to be similar to the sparse synthesis one. First, let Ω ∈ R p×d be a signal transformation or an analysis operator. Its rows are the row vectors {ω j } p j=1 that will be applied to the signals. Applying Ω to x, we obtain the (analysis) representation Ωx of x. To capture various aspects of the information in x, we typically have p ≥ d.

For simplicity, unless stated otherwise, we shall assume hereafter that all the rows of Ω are in general position, i.e., there are no non-trivial linear dependencies among the rows. 2 Note that this assumption is used for the purpose of contrasting the analysis model to the synthesis model, and that we will study a case when Ω is not in general position in later sections. As a matter of fact, there is some indication that linear dependecies among the rows of Ω can be a 'blessing.' (See, e.g., uniqueness results in Sections 3.3 and 3.4.)

Clearly, unless x = 0, no representation Ωx can be 'very sparse', since at least p -d of the coefficients of Ωx are necessarily non-zeros. We shall put our emphasis on the number of zeros in the representation, a quantity we will call cosparsity.

Definition 1. The cosparsity of a signal x ∈ R d with respect to Ω ∈ R p×d (or simply the cosparsity of x) is defined to be:

Cosparsity : ℓ := p -Ωx 0 (1)
The index set of the zero entries of Ωx is called the cosupport of x. We say that x has cosparse representation or x is cosparse when the cosparsity of x is large, where by large we mean that ℓ is close to d. We will see that, while ℓ ≤ d for an analysis operator in general position, there are specific examples where ℓ may exceed d.

At first sight the replacement of sparsity by cosparsity might appear to be mere semantics. However we will see that this is not the case. In the synthesis model it is the columns d j , j ∈ T associated with the index set T of nonzero coefficients that define the signal subspace. Removing columns from D not in T leaves this subspace unchanged. In contrast, it is the rows ω j associated with the index set Λ such that ω j , x = 0, j ∈ Λ that define the analysis subspace. In this case removing rows from Ω for which ω j , x = 0 leaves the subspace unchanged.

From this perspective, the cosparse model is more related to signal characterizations from the zero-crossings of their undecimated wavelet transform [START_REF] Mallat | Zero-crossings of a wavelet transform[END_REF] than to sparse wavelet expansions.

Sparse Analysis Model as a Generative Model

In a Bayesian context, one can think of data models as generators for random signals from a pre-specified probability density function. In that context, the signals that satisfy the k-sparse synthesis model can be generated as follows: First, choose k distinct columns of the dictionary D at random (e.g. assuming a uniform probability). We denote the index set chosen by T , and clearly |T | = k. Second, form a coefficient vector z that is k-sparse, with zeros outside the support T . The k non-zeros in z can be chosen at random as well (e.g. Gaussian iid entries). Finally, the signal is created by multiplying D to the resulting sparse coefficient vector z.

Could we adopt a similar view for the cosparse analysis model? The answer is positive. Similar to the above, one can produce an ℓ-cosparse signal in the following way: First, choose ℓ rows of the analysis operator Ω at random, and those are denoted by an index set Λ (thus, |Λ| = ℓ). Second, form an arbitrary signal v in R d -e.g., a random vector with Gaussian iid entries. Then, project v onto the orthogonal complement of the subspace generated by the rows of Ω that are indexed by Λ, this way getting the cosparse signal x. Explicitly,

x = (Id -Ω T Λ (Ω Λ Ω T Λ ) -1 Ω Λ )v.
Alternatively, one could first find a basis for the orthogonal complement and then generate a random coefficient vector for the basis.

This way, both models can be considered as generators of signals that have a special structure, and clearly, the two signal generators are different. It is now time to ask how those two families of signals inter-relate. In order to answer this question, we take the union-of-subspaces point of view.

Union-of-Subspaces Models

It is well known that the sparse synthesis model is a special instance of a wider family of models called union-of-subspaces [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF][START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF]. Given a dictionary D, a vector z that is exactly k-sparse with support T leads to a signal x = Dz = D T z T , a linear combination of k columns from D. The notation D T denotes the sub-matrix of D containing only the columns indexed by T . Denoting the subspace spanned by these columns by V T := span(d j , j ∈ T ), the sparse synthesis signals belong to the union of all n k possible subspaces of dimension k,

Sparse Synthesis Model:

x

∈ ∪ T :|T |=k V T . (2)
Similarly, the analysis model is associated to a union of subspaces model as well. Given an analysis operator Ω, a signal that is exactly ℓ-cosparse with respect to the rows Λ from Ω is simply in the orthogonal complement to these ℓ rows. Thus, we have3 Ω Λ x = 0, which implies that x ∈ W Λ , where W Λ := span(ω j , j ∈ Λ) ⊥ = {x, ω j , x = 0, ∀j ∈ Λ}. Put differently, we may write

W Λ = Range(Ω T Λ ) ⊥ = Null(Ω Λ ).
Hence, cosparse analysis signals x belong to the union of all the p ℓ possible such subspaces of dimension d -ℓ, Cosparse Analysis Model:

x ∈ ∪ Λ:|Λ|=ℓ W Λ . (3)
The following table summarizes these two unions of subspaces, where we recall that we assume Ω and D in general position.

Model Subspaces

No. of Subspaces Subspace dimension Synthesis V T := span(d j , j ∈ T )

n k k Analysis W Λ := span(ω j , j ∈ Λ) ⊥ p ℓ d -ℓ
What is the relation between these two union of subspaces, as described in Equations ( 2) and (3)? In general, the answer is that the two are different. An interesting way to compare the two models is to consider an ℓ-cosparse analysis model and a corresponding (d -ℓ)-sparse synthesis model, so that the two have the same dimension in their subspaces.

Following this guideline, we consider first a special case where ℓ = d -1. In such a case, the dimension of the analysis subspaces is d -ℓ = 1, and there are p ℓ of those. An equivalent synthesis union of subspaces can be created, where k = 1. We should construct a dictionary D with n = p ℓ atoms d j , where each atom is the orthogonal complement to one of the sets of ℓ rows from Ω. While the two models become equivalent in this case, clearly n ≫ p in general, implying that the sparse synthesis model becomes intractable since D becomes too large.

By further assuming that p = d, we get that there are exactly p ℓ = d d-1 = d subspaces in the analysis union, and in this case n = p = d as well. Furthermore, it is not hard to see that in this case the synthesis atoms are obtained directly by a simple inversion, D = Ω -1 .

Adopting a similar approach, considering the general case where ℓ is a general value (and not necessarily d -1), one could always construct a synthesis model that is equivalent to the analysis one. We can compose the synthesis dictionary by simply concatenating all the bases for the orthogonal complements to the subspaces W Λ . The obtained dictionary will have at most (d -ℓ) p ℓ atoms. However, not all supports of size k are allowed in the obtained synthesis model, since otherwise the new sparse synthesis model will strictly contain the cosparse analysis one. As such, the cosparse analysis model may be viewed as a sparse synthesis model with some structure.

Further on the comparison between the two models, it would be of benefit to consider again the case d -ℓ = k (i.e., having the same dimensionality), assume that p = n (i.e., having the same overcompleteness, for example with Ω = D T ), and compare the number of subspaces amalgamated in each model. For the sake of simplicity we consider a mild overcompleteness of p = n = 2d. Denoting H(t) := -t log 2 t -(1 -t) log 2 (1 -t), 0 < t < 1, the number of subspaces of low dimension k ≪ d = n/2 in each data model, from Stirling's approximation, roughly satisfies for large d:

Synthesis: log 2 n k ≈ n • H k n ≈ k • log 2 n k Analysis: log 2 p ℓ ≈ n • H d -k n ≈ n • H(0.5) = n.
More generally, unless d/n ≈ 1, there are far fewer low-dimensional synthesis subspaces than there are analysis subspaces of the same dimension. This is illustrated on Figure 2.3 when n = p = 2d. This indicates a strong difference in the structure of the two models: The synthesis model includes very few low-dimensional subspaces, and an increasingly large number of subspaces of higher dimension, whereas the analysis model contains a combinatorial number of low-dimensional subspaces, with fewer high dimensional subspaces. sults in Section 6.1). This stems from the fact that, in general, it is not possible to get cosparsity d ≤ ℓ < p: any vector x that is orthogonal to d linearly independent rows of Ω must be the zero vector, leading to an uninformative model. One may, however, get cosparsities in the range d ≤ ℓ < p when the analysis operator Ω displays certain linear dependencies. Therefore it appears to be desirable, in the cosparse analysis model, to have analysis operators that exhibit high linear dependencies among their rows. We will see in Section 3.4 that a leading example of such operators is the finite difference analysis operator.

Comment
Another interesting point of view towards the difference between the two models is the following: While a synthesis signal is characterized by the support of the non-zeros in its representation in order to define the subspace it belong to, a signal from the analysis model is characterized by the locations of the zeros in its representation Ωx. The fact that this representation may contain many non-zeroes (and especially so when p ≫ d) should be of no consequence to the efficiency of the analysis model.

Comparison with the Traditional Sparse Analysis model

Previous work using analysis representations, both theoretical and algorithmic, has focussed on gauging performance in terms of the more traditional sparsity perspective. For example, in the context of compressed sensing, recent theoretical work [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] has provided performance guarantees for minimum ℓ 1 -norm analysis representations in this light.

The analysis operator is generally viewed as the dual frame for a redundant synthesis dictionary so that Ω = D † . This means that the analysis coefficients Ωx provide a consistent synthesis representation for x in terms of the dictionary D, implying that the representation Ωx is a feasible solution to the linear system of equations Dz = x.

Furthermore, if Ωx 0 = p -ℓ, then Ωx must be an element of the k-sparse

synthesis model, T :|T |=k V T , with k = p -ℓ. Hence: {0} ⊆ Λ:|Λ|=p-k W Λ ⊆ T :|T |=k V T ⊆ R d . (4)
Of course, Ωx is not guaranteed to be the sparsest representation of x in terms of D. Hence the two subspace models are not equivalent. Note that while in Section 2.3 the sparsity k was matched to d -ℓ, here it is matched to p -ℓ. The former was used to get the same dimensions in the resulting subspaces, while the match discussed here considers the vector Ωx as a candidate k-sparse representation.

Such a perspective treats the analysis operator as a poor man's sparse synthesis representation. That is, for certain signals x, the representation Ωx may be reasonably sparse but is unlikely to be as sparse as, for example, the minimum ℓ 1 -norm synthesis representation 4 .

In the context of linear inverse problems, it is tempting to try to exploit the nesting property (4) in order to derive identifiability guarantees in terms of the sparsity of the analysis coefficients Ωx. For example, in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF], the compressed sensing recovery guarantees exploit the nesting property (4) by assuming a sufficient number of observations to achieve a stable embedding (restricted isometry property) for the k-sparse synthesis union of subspaces, which in turn implies a stable embedding of the (p -k)-cosparse analysis union of subspaces.

While such an approach is of course valid, it misses a crucial difference between the analysis and synthesis representations: they do not correspond to equivalent signal models. Treating the two models as equivalent hides the fact that they may be composed of subspaces with markedly different dimensions. The difference between these models is highlighted in the following examples.

Example: generic analysis operators, p = 2d

Assuming the rows of Ω are in general position, then when p ≥ 2d the nesting property (4) is trivial but rather useless! Indeed, if k < d, then the only analysis signal for which Ωx 0 = k = p -ℓ is x = 0. Alternatively, if k ≥ d, the synthesis model is trivially the full space:

T :|T |=k V T = R d .

Example: shift invariant wavelet transform

The shift invariant wavelet transform is a popular analysis transform in signal processing. It is particularly good for processing piecewise smooth signals.

Its inverse transform has a synthesis interpretation as the redundant wavelet dictionary consisting of wavelet atoms with all possible shifts.

The shift invariant wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF] provides a nice example of an analysis operator that has significant dependencies due to the finite support of the individual wavelets. Such nontrivial dependencies within the rows of Ω WT mean that the dimensions of the (analysis or synthesis) signal subspaces are not easily characterised by either the sparsity k or the cosparsity ℓ. However the behaviour of the model is still driven by the zero coefficients not the nonzero ones, i.e., by the zero-crossings of the wavelet transform [START_REF] Mallat | Zero-crossings of a wavelet transform[END_REF]. By considering a particular support set of an analysis representation Ω WT x with the shift invariant wavelet transform we can illustrate the dramatic difference between the analysis and synthesis interpretations of the coefficients. Figure 2 shows the support set of the nonzero analysis coefficients (white region), associated with the cone of influence around a discontinuity in a piecewise polynomial signal of length 128-samples [START_REF] Dragotti | Wavelet footprints: theory, algorithms, and applications[END_REF], using a shift-invariant Daubechies wavelet transform with s = 3 vanishing moments [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF]. For such a signal, the cone of influence at level J in a shift invariant wavelet transform contains L j -1 nonzero coefficients where L j is the length of the wavelet filter at level j. Note though, the nonzero coefficients are not linearly independent and can be elegantly described through the notion of wavelet footprints [START_REF] Dragotti | Wavelet footprints: theory, algorithms, and applications[END_REF].

Synthesis perspective. Interpreting the support set within the synthesis model implies that the signal is not particularly sparse and needs a significant number of wavelet atoms to describe it: in Figure 2 the size of the support set, excluding coefficients of scaling functions, is 122. Could the support set be significantly reduced by using a better support selection strategy such as ℓ 1 minimization? In practice, using ℓ 1 minimization, a support set of 30 can be obtained, again ignoring scaling coefficients.

Analysis perspective. The analysis interpretation of the shift invariant wavelet representation relies on the examination of the size of the analysis subspace associated with the cosupport set. From the theory of wavelet footprints, the dimension of this subspace is equal to the number of vanishing moments of the wavelet filter, which in this example is only . . . 3, providing a much lower dimensional signal model.

We therefore see that the analysis model has a much lower number of degrees of freedom for this support set, leading to a significantly more parsimonious model.

Hybrid Analysis/Synthesis models?

In this section we have demonstrated that while both the cosparse analysis model and the sparse synthesis model can be described by a union of subspaces these models are typically very different. We do not argue that one is inevitably better than the other. The value of the model will very much depend on the problem instance. Indeed the intrinsic difference between the models also suggests that it might be fruitful to explore building other union of subspace models from hybrid compositions of analysis and synthesis operators. For example, one could imagine a signal model where x = Dz through a redundant synthesis dictionary but instead of imposing sparsity on z we restrict z through an additional analysis operator: Ωz 0 ≤ k. In such a case there will still be an underlying union of subspace model but with the subspaces defined by a combination of atoms and analysis operator constraints. A special case of this is the split analysis model suggested in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF].

Uniqueness Properties

In the synthesis model, if a dictionary D is redundant, then a given signal x can admit many synthesis representations z, i.e., z with Dz = x. This makes the following type of problem interesting in the context of the sparse signal recovery: When a signal has a sparse representation z, can there be another representation that is equally sparse or sparser? This problem is well-understood in terms of the so-called spark of D [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF], the smallest number of columns from D that are linearly dependent.

Unlike in the synthesis model, if the signal is known, then its analysis representation Ωx with respect to an analysis operator Ω is completely determined. Hence, there is no inherent question of uniqueness for the cosparse analysis model. The uniqueness question we want to consider in this paper is in the context of the noiseless linear inverse problem,

y = Mx, (5)
where M ∈ R m×d , and m < d, implying that the measurement vector y ∈ R m is not sufficient to fully characterize the original signal x ∈ R d . For this problem we ask: when can we assert that a solution x with cosparsity ℓ is the only solution with that cosparsity or more? The problem (5) (especially, with additive noise) arises ubiquitously in many applications, and we shall focus on this problem throughout this paper as a platform for introducing the cosparse analysis model, its properties and behavior. Not to complicate matters unnecessarily, we assume that all the rows of M are linearly independent, and we omit noise, leaving robustness analysis to further work. For completeness of our discussion, let us return for a moment to the synthesis model and consider the uniqueness property for the inverse problem posed in Equation [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. Assuming that the signal's sparse representation satisfies x = Dz, we have that y = Mx = MDz. Had we known the support T of z, this linear system would have reduced to y = MD T z T , a system of m equations with k unknowns. Thus, recovery of x from y is possible only if k ≤ m.

When the support of z is unknown, it is the spark of the compound matrix MD that governs whether the cardinality of z T is sufficient to ensure uniqueness-if k = z 0 is smaller than half the spark of MD, then necessarily z is the signal's sparsest representation. At best, Spark(MD) = m + 1, and then we require that the number of measurements is at least twice the cardinality k. Put formally, we require

k = z 0 < 1 2 Spark(MD) ≤ m + 1 2 . ( 6
)
It will be interesting to contrast this requirement with the one we will derive hereafter for the analysis model.

Uniqueness When the Cosupport is Known

Before we tackle the uniqueness problem for the analysis model, let us consider an easier question: Given the observations y obtained via a measurement matrix M, and assuming that the cosupport Λ of the signal x is known, what are the sufficient conditions for the recovery of x? The answer to this question is straightforward since x satisfies the linear equation

y 0 = M Ω Λ x = Ax. (7)
To be able to uniquely identify x from Equation ( 7), the matrix A must have a zero null space. This is equivalent to the requirement

Null(Ω Λ ) ∩ Null(M) = W Λ ∩ Null(M) = {0}. (8)
Let us now assume that M and Ω are mutually independent, in the sense that there are no nontrivial linear dependencies among the rows of M and Ω; this is a reasonable assumption because first, one should not be measuring something that may be already available from Ω, and second, for a fixed Ω, mutual independency holds true for almost all M (in the Lebesgue measure). Then, (8) would be satisfied as soon as dim

(W Λ ) + dim(Null(M)) ≤ d, or dim(W Λ ) ≤ m, since dim(Null(M)) = d -m. This motivates us to define κ Ω (ℓ) := max |Λ|≥ℓ dim(W Λ ). (9)
The quantity κ Ω (ℓ) plays an important role in determining the necessary and sufficient cosparsity level for the identification of cosparse signals. Indeed, under the assumption of the mutual independence of Ω and M, a necessary and sufficient condition for the uniqueness of every cosparse signal given the knowledge of its cosupport Λ of size ℓ is κ Ω (ℓ) ≤ m.

(10)

Uniqueness When the Cosupport is Unknown

The uniqueness question that we answered above refers to the case where the cosupport is known, but of course, in general this is not the case. We shall assume that we may only know the cosparsity level ℓ, which means that our uniqueness question now becomes: what cosparsity level ℓ guarantees that there can be only one signal x matching a given observation y?

As we have seen, the cosparse analysis model is a special case of a general union of subspaces model. Uniqueness guarantees for missing data problems such as [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF] with general union of subspace models are covered in [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF][START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF]. In particular [START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF] shows that M is invertible on the union of subspaces ∪ γ∈Γ S γ if and only if M is invertible on all subspaces S γ + S θ for all γ, θ ∈ Γ. In the context of the analysis model this gives the following result whose proof is a direct consequence of the results in [START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF]:

Proposition 2 ([33]). Let ∪ Λ W Λ , |Λ| = ℓ
be the union of ℓ-cosparse analysis subspaces induced by the analysis operator Ω. Then the following statements are equivalent:

1. If the linear system y = Mx admits an ℓ-cosparse solution, then this is the unique ℓ-cosparse solution;

2. M is invertible on ∪ Λ W Λ ; 3. (W Λ1 + W Λ2 ) ∩ Null(M) = 0 for any |Λ 1 |, |Λ 2 | ≥ ℓ;
Proposition 2 answers the question of uniqueness for cosparse signals in the context of linear inverse problems. Unfortunately, the answer we obtained still leaves us in the dark in terms of the necessary cosparsity level or necessary number of measurements. In order to pose a clearer condition, we use Proposition 2 from [START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF] that poses a sharp condition on the number of measurements to guarantee uniqueness (when M and Ω are mutually independent): m ≥ κΩ (ℓ), where κΩ (ℓ) := max {dim(W Λ1 + W Λ2 ) :

|Λ i | ≥ ℓ, i = 1, 2} (11)
Interestingly, a sufficient condition can also be obtained using the quantity κ Ω defined in (9) above, which was observed to play an important role in the uniqueness result when the cosupport is assumed to be known. Namely, we have the following result.

Proposition 3. Assume that κ Ω (ℓ) ≤ m 2 .
Then for almost all M (wrt the Lebesgue measure), the linear inverse problem y = Mx has at most one ℓcosparse solution.

Proof. Assuming the mutual independence of Ω and M, which holds for almost all M, we note that the uniqueness of ℓ cosparse solutions holds if and only

if: dim (W Λ1 + W Λ2 ) ≤ m, whenever |Λ i | ≥ ℓ, i = 1, 2. Assume that κ Ω (ℓ) ≤ m/2. By definition of κ Ω , if |Λ i | ≥ ℓ, i = 1, 2, then dim(W Λi ) ≤ m 2 , hence dim (W Λ1 + W Λ2 ) ≤ m.
In the synthesis model the degree to which columns are interdependent can be partially characterized by the spark of D [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF] defined as the the smallest number of columns of D that are linearly dependent. Here the function κ Ω plays a similar role in quantifying the interdependence between rows in the analysis model.

Remark 4. The condition κ Ω (ℓ) ≤ m
2 is in general not necessary while condition [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] is.

There are two classes of analysis operators for which the function κ Ω is well-understood: analysis operators in general position and the finite difference operators. We discuss the uniqueness results for these two classes in the following subsections.

Analysis Operators in General Position

It can be easily checked that κ Ω (ℓ) = max(d -ℓ, 0). This enables us to quantify the exact level of cosparsity necessary for the uniqueness guarantees: Corollary 5. Let Ω ∈ R p×d be an analysis operator in general position. Then, for almost all m × d matrices M, the following hold:

• Based on Eq. ( 10), if m ≥ d -ℓ, then the equation y = Mx has at most one solution with known cosupport Λ (of cosparsity at least ℓ);

• Based on Proposition 2, if m ≥ 2(d -ℓ)
, then the equation y = Mx has at most one solution with cosparsity at least ℓ.

The Finite Difference Operator

An interesting class of analysis operators with significant linear dependencies is the family of finite difference operators on graphs, Ω DIF . These are strongly related to TV norm minimization, popular in image processing applications [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and has the added benefit that we are able to quantify the function κ Ω and hence the uniqueness properties of the cosparse signal model under Ω DIF .

We begin by considering Ω DIF on an arbitrary graph before restricting our discussion to the 2D lattice associated with image pixels. Consider a nonoriented graph with vertices V and edges E ⊂ V 2 . An edge e is a pair e = (v 1 , v 2 ) of connected vertices. For any vector of coefficients defined on the vertices, x ∈ R V , the finite difference analysis operator Ω DIF computes the collection of differences (x(v 1 ) -x(v 2 )) between end-points, for all edges in the graph. Thus, an edge e ∈ E may be viewed as a finite difference on R V .

Can we estimate the function κ ΩDIF (ℓ)?

The following shows that it is intimately related to topological properties of the graph. For each sub-collection Λ ⊂ E of edges, we can define its vertex-set V (Λ) ⊂ V as the collection of vertices covered by at least one edge in Λ. The support set V (Λ) of Λ can be decomposed into J(Λ) connected components (a connected component is a set of vertices connected to one another by a walk through vertices in Λ). It is easy to check that a vector x belongs to the space W Λ = Null(Ω Λ ) if and only if its values are constant on each connected component. As a result, the dimension of this subspace is given by dim

(W Λ ) = |V | -|V (Λ)| + J(Λ)
where the |V | -|V (Λ)| vertices out of V are associated to arbitrary values in x that are distinct from all their neighbors, while all entries from each of the J(Λ) connected components have an arbitrary common value. It follows that

κ Ω (ℓ) = max |Λ|≥ℓ |V | -|V (Λ)| + J(Λ) = |V | -min |Λ|≥ℓ |V (Λ)| -J(Λ) (12)
Because of the nesting of the subspaces W Λ , the minimum on the right hand side is achieved when |Λ| = ℓ.

Uniqueness Condition for Cosparse Images with respect to the 2D Ω DIF . In the abstract context of general graphs, the characterization (12) may remain obscure, but can we get more concrete estimates by specializing to the 2D regular graph associated to the pixels of an N × N image? It turns out that one can obtain relatively simple upper and lower bounds for κ ΩDIF and hence derive an easily interpretable uniqueness condition (see Appendix C for a proof):

Proposition 6.
Let Ω DIF be the 2D finite difference analysis operator that computes horizontal and vertical discrete derivatives of a d = N × N image.

For any ℓ we have

d - ℓ 2 - ℓ 2 -1 ≤ κ ΩDIF (ℓ) ≤ d - ℓ 2 - ℓ 2 + 1 2 (13)
As a result, assuming that M is 'mutually independent' from Ω DIF , and if

m ≥ 2d -ℓ - √ 2ℓ + 1 2 ≥ 2κ ΩDIF (ℓ) ( 14
)
then the equation y = Mx has at most one solution with cosparsity at least ℓ.

The 2D Ω DIF , Piecewise Constant Images, and the TV norm. The 2D finite difference operator is closely related to the TV norm [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]: the discrete TV norm of x is essentially a mixed ℓ 2 -ℓ 1 norm of Ω DIF x. Just like its close cousin TV norm minimization, the minimization of Ωx 0 is particularly good at inducing piecewise constant images. We illustrate this through a worked example. Consider the popular Shepp Logan phantom image shown in the left-hand side of Figure 3. We denote the cosupport of the image by Λ in line with the discussion in this section: An edge belongs to Λ if a pair of horizontally or vertically neighboring pixels v 1 and v 2 have the same value. This particular image has 14 distinct connected regions of constant intensity. The number of non-zero coefficients in the finite difference representation is determined by the total length (Manhattan distance) of the boundaries between these regions. For the Shepp Logan phantom this length is 2546 pixel widths and thus the cosparsity is ℓ = 130560 -2546 = 128014. Furthermore, as there are no isolated pixels with any other intensity, all pixels belong to a constant intensity region so that |V (Λ)| = |V | and the cosupport has an associated subspace dimension of: In order to determine when the Shepp Logan image is the unique solution to y = Mx with maximum cosparsity it is necessary to consider the maximum subspace dimension of all possible support sets with the same cosparsity. This is the quantity measured by κ ΩDIF (ℓ).

dim(W Λ ) = (|V | -|V (Λ)|) + J(Λ) = 14
Following the arguments used in the proof of Proposition 6 we need to find an image for which the Ω DIF cosupport, Λ, is a single connected subgraph that is as close to a square as possible. Such an image is shown in the right hand side of Figure 3. For this image we have dim(W Λ ) = 1276. Comparing this to the bounds given in (13) of Proposition 6 1275 ≤ κ ΩDIF (ℓ) ≤ 1276, we see that in this instance the upper bound has been achieved. The uniqueness result from Proposition 6 then tells us that a sufficient number of measurements to uniquely determine the Shepp Logan image is given by m = 2κ ΩDIF (128014) = 2552.

We will revisit this again in Subsection 6.2 where we investigate the empirical recovery performance of some practical reconstruction algorithms.

Overview of cosparse vs sparse models for inverse problems

To conclude this section, Figure 4 provides a schematic overview of analysis cosparse models vs synthesis sparse models in the context of linear inverse problems such as compressed sensing. In the synthesis model, the signal x is a projection (through the dictionary D) of a high-dimensional vector z living in the union of sparse coefficient subspaces; in the analysis model, the signal lives in the pre-image by the analysis operator Ω of the intersection between the range of Ω and this union of subspaces. For a given sparsity of z, this is usually a set of much smaller dimensionality.

Pursuit algorithms

Having a theoretical foundation for the uniqueness of the problem x = arg min

x Ωx 0 subject to Mx = y, (15)
we now turn to the question of how to solve it: algorithms. We present two algorithms, both targeting the solution of problem [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF]. As in the uniqueness discussion, we assume that M ∈ R m×d , where m < d. This implies that the equation Mx = y has infinitely many possible solutions, and the term Ωx 0 introduces the analysis model to regularize the problem.

The Cosparse Signal Recovery Problem is NP-complete

Related to [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF], we can consider a cosparse signal recovery problem COSPARSE consisting of a quintuplet (y, M, Ω, ℓ, ǫ) in which we seek to find a vector x * that satisfies

y -Mx * 2 ≤ ǫ, Ωx * 0 ≤ p -ℓ ( 16
)
where p is the number of rows of Ω as before. It is easy to see that the decision problem "given (y, M, Ω, ℓ, ǫ), does there exist x * satisfying (16)?" is NP [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF]:

given a candidate solution, one can check in polynomial time whether it satisfies the constraints [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF]. Moreover, every instance of the classical NP-complete (ǫ, k) SPARSE approximation problem [START_REF] Davis | Adaptive greedy approximations[END_REF][START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF] can trivially be reduced to an instance of the above decision problem with Ω = Id, hence COSPARSE is indeed NPcomplete.

The above consideration prompts us to look for ways to solve [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF] in an 'approximate' way. We discuss two possibilities, a convex relaxation and a greedy approach, with an emphasis on the latter. Of course, there can be many more possibilities to solve [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF] or to find approximate solutions of it. We mention a few works where some of such methods can be found: [START_REF] Cai | Split bregman methods and frame based image restoration[END_REF][START_REF] Portilla | Image restoration through l0 analysis-based sparse optimization in tight frames[END_REF][START_REF] Selesnick | Signal restoration with overcomplete wavelet transforms: Comparison of analysis and synthesis priors[END_REF].

The Analysis ℓ 1 -minimization

A natural convex relaxation of ( 15) is to solve:

x = arg min x Ωx 1 subject to Mx = y. (17)
The analysis ℓ 1 -minimization is well-known and widely used already in practice (see e.g. [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (mca)[END_REF][START_REF] Starck | Sparse Image and Signal Processing -Wavelets, Curvelets, Morphological Diversity[END_REF]). The attractiveness of this approach comes from that the convexity of ( 17) admits computationally tractable algorithms to solve the problem, and that the ℓ 1 -norm promotes high cosparsity in the solution x. An algorithm that targets the solution of ( 17) and its convergence analysis can be found in [START_REF] Cai | Split bregman methods and frame based image restoration[END_REF]. There are many other papers that have introduced algorithms to solve problems of the form [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF] or variants thereof. To give just an example, [START_REF] Raguet | Generalized forward-backward splitting[END_REF] proposes a general form of forward-backward splitting that can be exploited to deal with such problems.

The Greedy Analysis Pursuit Algorithm (GAP)

The algorithm we present in this section is named Greedy Analysis Pursuit (GAP). It is a variant of well-known greedy pursuit algorithm used for the synthesis model-the Orthogonal Matching Pursuit (OMP) algorithm. Similar to the synthesis case, our goal is to detect the informative support of Ωx-as discussed in Section 3.1, in the analysis case, this amounts to the locations (cosupport) of the zeros in the vector Ωx, so as to introduce additional constraints to the underdetermined system Mx = y. Note that for obtaining a solution, one needs to detect at least d-m of these zeros, and thus if ℓ > d-m, detection of the complete set of zeros is not mandatory.

An obvious way to find the cosupport of a cosparse signal would proceed as follows: First, obtain a reasonable estimate of the signal from the given information. Using the initial estimate, select a location as belonging to the cosupport. Having this estimated part of the cosupport, we can obtain a new estimate. One can now see that by alternating the two previous steps, we will have estimated enough locations of the cosupport to get the final estimate.

However, the GAP works in an opposite direction and aims to detect the elements outside the set Λ, this way carving its way towards the detection of the desired cosupport. Therefore, the cosupport Λ is initialized to be the whole set {1, 2, 3, . . . , p}, and through the iterations it is reduced towards a set of size ℓ (or less, d -m).

Let us discuss the algorithm with some detail. First, the GAP uses the following initial estimate:

x0 = arg min x Ωx 2 2 subject to y = Mx. (18)
Not knowing the locations of the cosupport but knowing that many entries of Ωx 0 are zero, this is a reasonable first estimate of x 0 . Once we have x0 , we can view Ωx 0 as an estimate of Ωx 0 . Hence, we find the location of the largest entries (in absolute value) of Ωx 0 and regard them as not belonging to the cosupport. After this, we remove the corresponding rows from Ω and work with a reduced Ω. A detailed description of the algorithm is given in Figure 5.

• Task: Approximate the solution of (15).

• Parameters: Given are the matrices M, Ω, the vector y, the target number of zeros ℓ, and a selection factor t ∈ (0, 1].

• Initialization: Set k = 0 and perform the following steps:

-Initialize Cosupport: Λ0 = {1, 2, 3, . . . , p}, -Initialize Solution:

x0 = arg min x Ω Λ0 x 2 2
subject to y = Mx.

• GAP Iterations: Increment k by 1 and perform the following steps: Some readers may notice that the GAP has similar flavors to the FOCUSS [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] and the IRLS [START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF]. This is certainly true in the sense that the GAP solves constrained least squares problems and adjusts weights as it iterates. However, the weight adjustment in the GAP is more aggressive (removal of rows) and binary in nature. We also note that the use of the selection factor t in the GAP is related to Weak Greedy Algorithms [START_REF] Temlyakov | Weak greedy algorithms[END_REF] for the synthesis model.

-Project: Compute α = Ωx k-1 , -Find largest entries: Γ k = {i : |αi| ≥ t maxj
Stopping criterion / targeted sparsity. In GAP, we have a range of choices between using the full ℓ zeros in the product Ωx versus a minimal and sufficient set of d -m zeros. In between these two values, and assuming that the proper elements of Λ have been detected, we expect the solution obtained by the algorithms to be the same, with a slightly better numerical stability for a larger number of zeros.

Thus, an alternative stopping criterion for the GAP could be to detect whether the solution is static or the analysis coefficients of the solution are small. This way, even if the GAP made an error and removed from Λk an index that belongs to the true cosupport Λ, the tendency of the solution to stabilize could help in preventing the algorithm to incorporate this error into the solution. In fact, this criterion is used in the experiment in Section 6.

Multiple selections. The selection factor 0 < t ≤ 1 allow the selection of multiple rows at once, to accelerate the algorithm by reducing the number of iterations.

Solving the required least squares problems. The solution of Eq. ( 18) (and of the adjusted problems with reduced Ω at subsequent steps of the algorithm)under some suitable conditions-is given analytically (see Appendix E for the derivation) by

x0 = (M T M + (Id -M T (MM T )M)Ω T Ω) -1 M T y. (19)
In practice, instead of ( 18), we compute

x0 = arg min x y -Mx 2 2 + λ Ωx 2 2 = arg min x y 0 - M √ λΩ x 2 2
for a small λ > 0, yielding the solution

x0 = M √ λΩ † y 0 = (M T M + λΩ T Ω) -1 M T y.
We point out that the use of the parameter λ is for convenience in this work but it becomes more useful when dealing with noisy observation. Furthermore, the small value of λ can have adverse effects on computational cost for some numerical algorithms (e.g., the conjugate gradient method) to solve the above minimization. In our implementation, we have used λ values roughly from 10 -4 to 10 -6 .

Theoretical analysis

So far, we have introduced the cosparse analysis data model, provided uniqueness results in the context of linear inverse problems for the model, and described some algorithms that may be used to solve such linear inverse problems to recover cosparse signals. Before validating the algorithms and the model proposed with experimental results, we first investigate theoretically under what conditions the proposed algorithms to solve cosparse signal recovery [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF] are guaranteed to work. After that, we discuss the nature of the condition derived by contrasting it to that for the synthesis model. Further discussion including some desirable properties of Ω and M can be found in Appendix D.

A Sufficient Condition for the Success of the ℓ 1 -minimization

In the sparse synthesis framework, there is a well-known necessary and sufficient condition called the null space property (NSP) [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF] that guarantees the success of the synthesis ℓ 1 -minimization ẑ0 := argmin z z 1 subject to y = Φz [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] to recover the sparsest solution, say z 0 , to y = Φz. To elaborate, in the case of a fixed support T , the ℓ 1 -minimization ( 20) recovers every sparse coefficient vector z 0 supported on T if and only if

z T 1 < z T c 1 , ∀z ∈ Null(Φ), z = 0. ( 21
)
The NSP [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] cannot easily be checked but some 'simpler' sufficient conditions can be derived from it; for example, one can get a recovery condition of [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] called the Exact Recovery Condition (ERC):

|Φ † T Φ T c | 1→1 < 1, (22)
where the notation |A| p→q denotes the operator norm of A from ℓ p to ℓ q , i.e.,

|A| p→q := sup x =0
Ax q x p .

The ERC ( 22) also implies the success of greedy algorithms such as OMP [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF].

Note that here we used the symbol Φ for an object which may be viewed as a dictionary or a measurement matrix. Separating the data model and sampling, we can write Φ = MD as was done in Section 3. One may naturally wonder: is there a condition for the cosparse analysis model that is similar to [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] and ( 22)? The answer to this question seems to be affirmative with some qualification as the following two results show (the proofs are in Appendix A): Theorem 7. Let Λ be a fixed cosupport. The analysis ℓ 1 -minimization

x0 := argmin x Ωx 1 subject to y := Mx 0 = Mx ( 23
)
recovers every x 0 with cosupport Λ as a unique minimizer if, and only if,

sup xΛ:ΩΛxΛ=0 | Ω Λ c z, sign(Ω Λ c x Λ ) | < Ω Λ z 1 , ∀z ∈ Null(M), z = 0. ( 24
)
Corollary 8. Let N T be any d×(d-m) basis matrix for the null space Null(M), and Λ be a fixed cosupport such that the ℓ ×

(d -m) matrix Ω Λ N T is of full rank d -m. If sup xΛ:ΩΛxΛ=0 (NΩ T Λ ) † NΩ T Λ c sign(Ω Λ c x Λ ) ∞ < 1, ( 25
)
then the analysis ℓ 1 -minimization [START_REF] Engan | Multi-frame compression: Theory and design[END_REF] recovers every x 0 with cosupport Λ. Moreover, if

|(NΩ T Λ ) † NΩ T Λ c | ∞→∞ = |Ω Λ c N T (Ω Λ N T ) † | 1→1 < 1 ( 26
)
then condition (25) holds true.

There is an apparent similarity between the analysis ERC condition (26) above and its standard synthesis counterpart ( 22), yet there are some subtle differences between the two that will be highlighted in Section 5.3.

A Sufficient Condition for the Success of the GAP

There is an interesting parallel between the synthesis ERC ( 22) and its analysis version in Corollary 8; namely, the analysis ERC condition (26) also implies the success of the GAP algorithm when the selection factor t of the GAP is 1

(in fact, |Ω Λ c N T (Ω Λ N T ) † | 1→1 < t ≤ 1)
, as we will now show.

From the way GAP algorithm works, we can guarantee that it will perform a correct elimination at the first step if the largest analysis coefficients of Ω Λ c x0 of the first estimate x0 are larger than the largest of Ω Λ x0 where Λ denotes the true cosupport of x 0 . This observation suggests that we can hope to find a condition for success if we can find some relation between Ω Λ c x0 and Ω Λ x0 .

The following result provides such a relation: Lemma 9. Let N T be any d × (d -m) basis matrix for the null space Null(M) and Λ be a fixed cosupport such that the ℓ × (d -m) matrix Ω Λ N T is of full rank d -m. Let a signal x 0 with Ω Λ x 0 = 0 and its observation y = Mx 0 be given. Then the estimate x0 in [START_REF] Dragotti | Wavelet footprints: theory, algorithms, and applications[END_REF] satisfies

Ω Λ x0 = -(NΩ T Λ ) † NΩ T Λ c Ω Λ c x0 . ( 27
)
Having obtained a relation between Ω Λ x0 and Ω Λ c x0 , we can derive a sufficient condition which guarantees the success of GAP for recovering the true target signal x 0 : Theorem 10. Let N T be any d×(d-m) basis matrix for the null space Null(M) and Λ be a fixed cosupport such that the ℓ × (d -m) matrix Ω Λ N T is of full rank d -m. Let a signal x 0 with Ω Λ x 0 = 0 and an observation y = Mx 0 be given. Suppose that the analysis ERC (26) holds true. Then, when applied to solve [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via ℓ 1 minimization[END_REF], GAP with selections factor t > |(NΩ T Λ ) † NΩ T Λ c | ∞→∞ will recover x 0 after at most |Λ c | iterations.

Proof. At the first iteration, GAP is doing the correct thing if it removes a row from Ω Λ c . Clearly, this happens when

Ω Λ x0 ∞ < t Ω Λ c x0 ∞ . ( 28
)
In view of [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming[END_REF], if [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] holds and t > |(NΩ T Λ ) † NΩ T Λ c | ∞→∞ , then ( 28) is guaranteed. Therefore, GAP successfully removes a row from Ω Λ c at the first step.

Now suppose that ( 26) was true and GAP has removed a row from Ω Λ c at the first iteration. Then, at the next iteration, we have the same Ω Λ and, in the place of Ω Λ c , a submatrix ΩΛ c of Ω Λ c (with one fewer row). Thus, we can invoke Lemma 9 again and we have

Ω Λ x1 = -NΩ T Λ † N ΩT Λ c ΩΛ c x1 . Let R 0 := NΩ T Λ † NΩ T Λ c and R 1 := NΩ T Λ † N ΩT Λ c
. We observe that R 1 is a submatrix of R 0 obtained by removing one column. Therefore,

|R 1 ∞→∞ ≤ |R 0 | ∞→∞ < t.
By the same logic as for the first step, the success of the second step is guaranteed. Repeating the same argument for the subsequent steps, we obtain the conclusion.

Clearly, at least one row from Λ c is removed at each iteration. Therefore, x 0 is recovered after at most |Λ c | iterations. [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] for the cosparse signal recovery guarantees the success of both the GAP and the analysis ℓ 1 -minimization.

Remark 11. As pointed out at the beginning of the subsection, the Exact Recovery Condition

Analysis vs synthesis exact recovery conditions

When Φ is written as MD, the exact recovery condition [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (mca)[END_REF] for the sparse synthesis model is equivalent to

|(MD T ) † MD T c | 1→1 < 1. ( 29
)
Here, T is the support of the sparsest representation of the target signal. At first glance, the two conditions ( 29) and ( 26):

|Ω Λ c N T (Ω Λ N T ) † | 1→1 < 1
look similar; that is, for both cases, one needs to understand the characteristics of a single matrix, ΩN T for the cosparse model, and MD for the sparse model. Moreover, the expressions involving these matrices have similar forms. However, upon closer inspection, there is a crucial difference in the structures of the two expressions. In the synthesis case, the operator norm in question depends only on how the columns of MD are related, since a more explicit writing of the pseudo-inverse shows that the matrix to consider is

(D T T M T MD T ) -1 (MD T ) T MD T c
This fact allows us to obtain more easily characterizable conditions like incoherence assumptions [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] that ensure condition [START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF].

To the contrary, in the analysis case, more complicated relations among the rows and the columns of ΩN T have to be taken into account. The matrix to consider being

Ω Λ c N T NΩ T Λ Ω Λ N T -1 NΩ T Λ , the inner expression NΩ T
Λ Ω Λ N T is connected with how the columns of ΩN T are related. However, because the matrices Ω Λ c N T and NΩ T Λ appear outside, it also becomes relevant how the rows of ΩN T are related.

There is also an interesting distinction in terms of the sharpness of these exact recovery conditions. Namely, the violation of ( 29) implies the failure of the OMP in the sense that there exist a sparse vector x = D T z T for which the first step of OMP picks up an atom which is not indexed by T . To the opposite, the violation of ( 26) does not seem to imply the necessary "failure" of GAP in a similar sense.

Note however that both conditions are not essential for the success of the algorithms. One of the reasons is that the violation of the conditions does not guarantee that the algorithms would select wrong atoms. Furthermore, even if the GAP or the OMP "fails" in one step, that does not necessarily mean that the algorithms fail in the end: further steps may still enable them to achieve an accurate estimate of the vector x 0 .

Relation to the Work by Candès et. al. [7]

Before moving on to experimental results, we discuss the recovery guarantee result of Candès et al. [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] for the algorithm x = argmin

x∈R d D T x 1 subject to Mx -y 2 ≤ ǫ (30)
when partial noisy observation y = Mx + w with w 2 ≤ ǫ is given for an unknown target signal x.

In order to derive the result, the concept of D-RIP is introduced [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF]:

A measurement matrix M satisfies D-RIP adapted to D with constant δ D s if (1 -δ D s ) v 2 2 ≤ Mv 2 2 ≤ (1 + δ D s ) v 2 2
holds for all v that can be expressed as a linear combination of s columns of D.

With this definition of D-RIP, the main result of [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] can be stated as follows:

For an arbitrary tight frame D and a measurement matrix M satisfying D-RIP with δ D 7s < 0.6, the solution x to (30) satisfies

x -x 2 ≤ C 0 ǫ + C 1 D T x -(D T x) s 1 √ s (31)
where the constants C 0 and C 1 may depend only on δ D 7s , and the notation (c) s represents a sequence obtained from a sequence c by keeping the s-largest values of c in magnitude (and setting the others to zero).

The above recovery guarantee is one of the few-very likely the only-results existing in the literature on [START_REF] Kluckner | Exploiting redundancy for aerial image fusion using convex optimization[END_REF]. However, we observe that there is much room for improving the result. We now discuss why we hold this view. For clarity and for the purpose of comparison to our result, we consider only the case ǫ = 0 for [START_REF] Kluckner | Exploiting redundancy for aerial image fusion using convex optimization[END_REF].

First, we note that [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] implicitly uses the estimate of type Ω Λ c z 1 < Ω Λ z 1 for [START_REF] Farsiu | Advances and challenges in super-resolution[END_REF]. Hence, the main result of [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] cannot be sharp in general due to the fact that the sign patterns of (24) are ignored. 5Second, the quality of the bound

D T x -(D T x) s 1 / √ s in (31) is measured
in terms of how effective D T x is in sparsifying the signal x with respect to the dictionary D. To explain, let us consider the synthesis ℓ 1 -minimization

∆ 1 (x) := argmin z∈R n z 1 subject to MDz = Mx (32)
and let ∆ 0 (x) be the sparsest representation of x with D:

∆ 0 (x) := argmin z∈R n z 0 subject to Dz = x.
Applying the standard result for the synthesis ℓ 1 -minimization, we have

∆ 1 (x) -∆ 0 (x) 2 ≤ C 2 ∆ 0 (x) -(∆ 0 (x)) s 1 √ s
provided that MD satisfies the standard RIP with, e.g.,

δ 2s < √ 2 -1 ≈ 0.414. Since D is a tight frame, it implies D∆ 1 (x) -x 2 ≤ C 2 ∆ 0 (x) -(∆ 0 (x)) s 1 √ s . (33)
Note that both ∆ 0 (x) and D T x are legitimate representations of x since D∆ 0 (x) = x = DD T x. Thus, ∆ 0 (x) is sparser than D T x in general; in this sense, D T x is not effective in sparsifying x. Given this, we expect that ∆ 0 (x) -

(∆ 0 (x)) s 1 / √ s is smaller than D T x -(D T x) s 1 / √ s.
We now see that [START_REF] Labate | Sparse multidimensional representation using shearlets[END_REF] with ǫ = 0 and (33) are of the same form. Furthermore, given the degree of restriction on the RIP constants (δ D 7s < 0.6 vs. δ 2s < 0.414), we can only expect that the constant C 2 is smaller than C 1 . From these considerations, [START_REF] Labate | Sparse multidimensional representation using shearlets[END_REF] suggests to us that analysis ℓ 1 -minimization (17) performs on par with synthesis ℓ 1 -minimization [START_REF] Llagostera Casanovas | Blind audiovisual source separation based on sparse representations[END_REF], or tends to perform worse.

Third, the only way for [START_REF] Labate | Sparse multidimensional representation using shearlets[END_REF] to explain that the cosparse signals are perfectly recovered by analysis ℓ 1 -minimization is to show that D T x is exactly s-sparse for some s > 0 with D-RIP constant δ D 7s < 0.6. Unfortunately, we can quickly observe that the situation becomes hopeless even for moderately overcomplete D; for example, let D be a 1.15-times overcomplete random tight frame for R d and consider recovering (d -1)-cosparse signals for the operator D T . Note that (d -1)-cosparse signals x lead to (0.15d + 1)-sparse representation D T x. This means that we need δ D 7(0.15d+1) = δ D 1.05d+7 to be smaller than 0.6 to show that x can be recovered with analysis ℓ 1 , which of course cannot happen since δ D d ≥ 1 (unless every element x in the span of D is uniquely characterized by its projection Mx, an uninteresting situation that can only occur if either M is one to one, or D does not span the signal space). By somehow taking the synthesis view of the signals, [START_REF] Labate | Sparse multidimensional representation using shearlets[END_REF] cannot explain the recovery of the simplest cosparse signals (cosparsity d -1) no matter what M is (as long as it is underdetermined).

We also observe that the result of [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] cannot say much about the recovery of cosparse signals with respect to the finite difference operators Ω DIF discussed in Section 3. This is due to the fact that Ω T DIF is not a tight frame. How does our recovery result [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] fare in this regard? For illustration, we took Ω to be the finite difference operator Ω DIF for 32 × 32 images (thus, d = 1024). As a test image, we took x to be constant in the region {(i, j) : i, j = 1, . . . , 16} and {(i, j) : i, j = 1, . . . , 16} c . For this admittedly simple test image, using the same notational convention as in Section 3.4, we computed the operator norm in [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] for random measurement matrices M ∈ R 640×1024 , with N a basis of the null space of M (computed with an SVD), and Λ the cosupport of the test image. When the operator norm was computed for 100 instances M, it was observed to be less than 0.726. Hence, our result does give the guarantee of cosparse signal recovery in simple cases.

Experiments

Empirical performance of the proposed algorithms is presented in this section. First, we show how the algorithms perform in synthetic cosparse recovery problems. Second, experimental results for an analysis-based compressed sensing are presented.

Performance of analysis algorithms

In this section, we apply the algorithms described in Section 4 to synthetic cosparse recovery problems. In the experiment, the entries of M ∈ R m×d were drawn independently from the normal distribution. The analysis operator Ω ∈ R p×d was constructed so that its transpose is a random tight frame with unit norm columns for R d -we will simply say that Ω is a random tight frame in this case. 6 A random (almost) tight frame B with unit columns was generated starting from a d × p Gaussian matrix by alternating the following two steps: 1) Singular value decomposition was performed on B to yield USV T = B, and then S was replaced with a matrix of the form [αId, 0] for α = p/d. This gives us a new B which is a tight frame. 2) The columns of B were normalized to unit length. Next, the cosparsity ℓ was chosen, and the true or target signal x was generated randomly as described in Section 2.2. The observation was obtained by y = Mx.

Matlab cvx package [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming[END_REF] with the default solver SeDuMi [START_REF] Sturm | Using sedumi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] was used for the analysis-ℓ 1 . The precision was set to best. For the final results, we used the estimate x from the ℓ 1 solver to obtain an estimate of the cosupport-the cosupport estimate was obtained by taking the indices for which the corresponding analysis coefficient is of size less than 10 -6 -and then using this cosupport and the observation y to compute the final estimate of x (this process can be considered as de-biasing.). No visually noticeable changes due to the last de-biasing step were noted in the results. Figure 6 shows the results. In all cases, the signal dimension d is set to 200. We then varied the number m of measurements, the cosparsity ℓ of the target signal, and the operator size p according to the following formulae:

σ = 1 σ = 1.2 σ = 2 GAP ρ 0.
m = δd, ℓ = d -ρm, p = σd.
which is consistent with Donoho & Tanner's notations for phase transition diagrams [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF]: δ = m/d is the undersampling ratio, and ρ = (d -ℓ)/m measures the relative dimension of the ℓ-cosparse subspaces compared to the number of measures. For every fixed parameter triplet (σ, δ, ρ), the experiment was repeated 50 times. A relative error of size less than 10 -6 was counted as perfect recovery. Each pixel in the diagrams corresponds to a triplet (σ, δ, ρ) and the pixel intensity represents the ratio of the signals recovered perfectly with white being the 100% success.

The figures show that the GAP can be a viable option when it comes to the cosparse signal recovery. GAP performs better than ℓ 1 -minimization, especially for overcomplete Ω's. It is clear from its description that GAP has polynomial complexity. In practice, computational cost can be high when the size of the problem is very large; to give a rough picture, for the experiment of Section 6.2, (a super greedy version of) GAP was observed to take twice or three times longer to complete the task than l1 magic.

It is known that the performance of OMP for the sparse signal recovery varies according to the nature of the distribution of the magnitudes of the sparse coefficients. More specifically, OMP performs very well when the coefficients follow independent Gaussian distributions while it does not work as well when the coefficients are drawn from independent Rademacher distributions (1 or -1 with equal probabilities). This phenomonon turns out to be true in the case of GAP as well, and Figure 7 shows the corresponding result when Ω is an orthogonal matrix. Note, however, such an unfavorable distribution as Rademacher may not make sense or may have different effects for redundant Ω's. An interesting phenomenon observed in the plots for overcomplete Ω is that there seems to be some threshold δ * such that if the observation to dimension ratio δ is less than δ * , one could not recover any signal however cosparse it may be. We may explain this heuristically as follows: If m measurements are available, then the amount of information we have for the signal is c 1 m where c 1 is the number of bits each observation represent. In order to recover a cosparse signal, we need first to identify which subspace the signal belongs to out of p ℓ , and then to obtain the d -ℓ coefficients for the signal with respect to a basis of the d -ℓ dimensional subspace. Therefore, roughly speaking, one may hope to recover the signal when

c 1 m ≥ log 2 p ℓ + c 1 (d -ℓ) = log 2 p ℓ + ρc 1 m.
Thus, the recovery is only possible when (1 -ρ)δ ≥ log 2 p d /(c 1 d). Using the relation p = σd and Stirling's approximation, this leads to an asymptotic relation

δ ≥ (1 -ρ)δ ≥ σ log σ -(σ -1) log(σ -1) c 1 ,
which explains the phenomenon. The calculation above and the experimental evidence from the figures confirm the intuition we had in Section 2.3: The combinatorial number of lowdimensional cosparse subspaces arising from analysis operators in general position is not desirable. This strengthens our view on the necessity of designing/learning analysis operators with high linear dependencies.

Analysis-based Compressed Sensing

We observed in Section 6.1 that the cosparse analysis model facilitates effective algorithms to recover partially observed cosparse signals. In this section, we demonstrate the effectiveness of GAP algorithm on a standard toy problem: the Shepp Logan phantom recovery problem.

We consider the following problem that is related to computed tomography (CT): There is an image, say of size n × n, which we are interested in but cannot observe directly. It can only be observed indirectly by means of its 2D Fourier transform coefficients. However, due to high cost of measurements or some physical limitation, the Fourier coefficients can only be observed along a few radial lines. These limited observations or the locations thereof can be modeled by a measurement matrix M, and with the obtained observation we want to recover the original image. As an ideal example, we consider the Shepp Logan phantom. One can easily see that this image is a good example of cosparse signals in Ω DIF which consists of all the vertical and horizontal gradients (or one step differences). This image has been used extensively as an example in the literature in the context of compressed sensing (see, e.g., [START_REF] Blumensath | Accelerated iterative hard thresholding[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]).

Figure 8 is the result obtained using GAP. The number of measurements that corresponds to 12 radial lines is m = 3032. Compared to the number of pixels in the image d = 65536, it is approximately 4.63%. The number of analysis atoms that give non-zero coefficients is p-ℓ = 2546. The size of Ω DIF is roughly twice the image size d = 65536, namely p = 130560. At first glance, this corresponds to very high cosparsity level (ℓ = 130560 -2546), or put differently, given the high cosparsity level ℓ = 128014, we seem to have required too many measurements. However, using the near optimal guarantee for uniqueness [START_REF] Do | The contourlet transform: An efficient directional multiresolution image representation[END_REF], we have a uniqueness guarantee when m ≥ 2552. In view of this, the fact that GAP recovered the signal perfectly for 3032 measurements is encouraging.

We have also ran the GAP algorithm for a larger sized 512 × 512 problem. The results (not shown here) are visually similar to Figure 8. In this case, the number of measurements (m = 7112) represents approximately 2.71% of the image size (d = 262144). The number of non-zero analysis coefficients is p -ℓ = 5104. The sufficient uniqueness condition [START_REF] Do | The contourlet transform: An efficient directional multiresolution image representation[END_REF] gives m ≥ 5110 as a number of measurements for the uniqueness.

While encouraged by the result on the Phantom images, we acknowledge that these images are unrealistic. Hence, further experiments on more realistic images will be desirable. Some progress has been made [START_REF] Nam | Recovery of cosparse signals with greedy analysis pursuit in the presence of noise[END_REF] to address this. conjugate gradient method) to approximate pseudo-inverses. Also, due to high computational cost, we eliminated many rows at each iteration (super greedy) instead of one. Although this was not implemented using a selection factor, this can be interpreted as using varying selection factors 0 < t k < 1 along the iterations.

To conclude this section, we have repeated the 256 × 256 Shepp Logan phantom image recovery problem using several algorithms while varying the number of radial observation lines. Given that we know the minimal theoretical number and a theoretically sufficient number of radial observation lines for the uniqueness guarantee, the experimental result gives us an insight on how various algorithms actually perform in the recovery problem in relation to the amount of observation available. Figure 9 shows the outcome. The algorithms used in the experiment are the GAP, the TV-minimization from l1magic, the AIHT from [START_REF] Blumensath | Accelerated iterative hard thresholding[END_REF], and the back-projection algorithm. 7 The GAP and l1magic can be viewed as analysis-based reconstruction algorithms while the AIHT is a synthesis-based reconstruction algorithm. The AIHT is seen to use Haar wavelets as the synthesis dictionary, hence the algorithm implicitly assumes that the phantom image has sparse representation in that dictionary. We remark that while Figure 9 gives an impression that the AIHT does not have any improvement over the baseline back-projection algorithm, perfect reconstructions were observed for the former when sufficient measurements were available, which is not the case for the back-projection. x 2

x-x 2 " where x is an approximation to the true signal x.

Remark 13. It must be noted that in our experiment, each radial line consists of N pixels for an N × N image; this is in contrast to the fact that the radial lines in the existing codes, e.g. l1magic, have N -1 pixels. We have made appropriate changes for our experiment. The radial lines with N -1 pixels do make the recovery problem more difficult and more observations were required for perfect recovery for the GAP.

Conclusions and Further Work

In this work, we have described the cosparse analysis data model as an alternative to the popular sparse synthesis model. We have shown that the cosparse analysis model is distinctly different from the sparse synthesis one in spite of their apparent similarities. In particular, treating the cosparse model as the synthesis model by assuming that the analysis representations of cosparse signals are sparse was demonstrated to be not very meaningful. Having presented the model, we stated conditions that guarantee the uniqueness of cosparse solutions in the context of linear inverse problems based on the work [START_REF] Lu | A theory for sampling signals from a union of subspaces[END_REF]. We then presented some algorithms for the cosparse recovery problem and provided some theoretical result for the analysis ℓ 1 -minimization and the newly proposed GAP. Lastly, the model and the proposed algorithm were validated via experimental results.

Although our work in this paper shows that the cosparse analysis model together with algorithms based on the model is an interesting subject to study and viable for practical applications, there is much more to be learned about the model. Among possible future avenues for related research, we list the following 8 :

1. The stability of measurement matrices M on the analysis union of subspaces ∪ Λ W Λ ;

2. The effect of noise on the cosparse analysis model and associated algorithms;

3. Adaptation to the cases where the signals of interest are not exactly cosparse;

4. The designing / learning of analysis operators for classes of signals of interest;

5. More concrete and/or optimal theoretical success guarantees for algorithms. As an example, one may seek a similar concept as coherence for the analysis operators;

6. Better understanding of the role of linear dependencies between rows of the analysis operator.

On top of these, one can also develop and study more algorithms for cosparse signal recovery; for example, if we view that the GAP is the 'dual' of the OMP, then we could ask: what is the dual of the CoSaMP [START_REF] Needell | Cosamp: iterative signal recovery from incomplete and inaccurate samples[END_REF] (or subspace pursuit [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF])?

Appendix A. Proof of Theorem 7 and Corollary 8

Let us begin with the simplest case. For a fixed x 0 with cosupport Λ, the analysis ℓ 1 -minimization (23) recovers x 0 as the unique minimizer if and only if

| Ω Λ c z, sign(Ω Λ c x 0 ) | < Ω Λ z 1 , ∀z ∈ Null(M), z = 0.
This follows from two facts: a) the above condition characterizes strict local minima of the optimization problem; b) the optimization problem is convex and can have at most one strict local minimum, which must be the unique global optimum. From this, we derive the following: The analysis ℓ 1 -minimization (23) recovers x 0 as a unique minimizer for any x 0 with cosupport Λ, if and only if

sup xΛ:ΩΛxΛ=0 | Ω Λ c z, sign(Ω Λ c x Λ ) | < Ω Λ z 1 , ∀z ∈ Null(M), z = 0
and the proof of Theorem 7 is complete.

8 Some progress [START_REF] Nam | Recovery of cosparse signals with greedy analysis pursuit in the presence of noise[END_REF] has been made for items 2 and 3 in the listing.

Since Ω Λ x 0 = 0, we can also write

NΩ T Λ Ω Λ u = -NΩ T Λ c Ω Λ c x0 (B.1)
with u = x0x 0 . On the other hand, from Mx 0 = y = Mx 0 , we have Mu = 0. This means that u can be expressed as u =: N T w for some w. Plugging this into (B.1), we have

NΩ T Λ Ω Λ N T w = -NΩ T Λ c Ω Λ c x0 . Hence, w = -NΩ T Λ Ω Λ N T -1 NΩ T Λ c Ω Λ c x0 . This gives us x0 -x 0 = u = -N T NΩ T Λ Ω Λ N T -1 NΩ T Λ c Ω Λ c x0 .
Again, using Ω Λ x 0 = 0, we have

Ω Λ x0 = -Ω Λ N T NΩ T Λ Ω Λ N T -1 NΩ T Λ c Ω Λ c x0 = -(NΩ T Λ ) † NΩ T Λ c Ω Λ c x0 .

Appendix C. Proof of Proposition 6

All the statements in this section are about a 2D regular graph consisting of d = N × N vertices (V ) and the vertical and horizontal edges (E) connecting these vertices. To prove the proposition, we will need three basic lemmas. Lemma 14. For a fixed ℓ, the value

α(ℓ) := min Λ⊂E:|Λ|≥ℓ {|V (Λ)| -J(Λ)}
is achieved for a subgraph (V (Λ), Λ)-we will simply identify Λ with the subgraph from here on-satisfying |Λ| = ℓ and J(Λ) = 1.

Proof. It is not difficult to check that the minimum if achieved for Λ with |Λ| = ℓ. Thus, we will assume |Λ| = ℓ. Now, we need to show that there is also a Λ with J(Λ) = 1. Suppose that Λ with | Λ| = ℓ achieves α(ℓ) and J( Λ) > 1. We will show that we can obtain Λ from Λ that also achieves the value α(ℓ), and |Λ| = ℓ and J(Λ) = 1. For simplicity, we will consider the case J( Λ) = 2 only; one can deal with other cases by the repetition of the same argument.

Let Λ1 and Λ2 be the two connected components of Λ. Note that on a 2D regular graph, we can shift a subgraph horizontally or vertically unless the subgraph has vertices on all four boundaries of V . Since Λ1 and Λ2 are disconnected, not all of them can have vertices on all four boundaries of V . Therefore, one of them, say Λ1 , can be shifted towards the other. Let us consider the first moment when they touched each other. Let t be the number of vertices that coincided. Then, at most t -1 edges must have coincided. Thus, denoting the number of edges coincided by s < t, the resulting subgraph Λ′ has |V ( Λ)| -t vertices and | Λ| -s edges and one connected components. Now let Λ be a subgraph obtained from Λ′ by adding s additional edges that are connected to Λ′ . Then,

|V (Λ)| ≤ |V ( Λ′ )| + s ≤ |V ( Λ)| -t + s, |Λ| = | Λ| = ℓ, and J(Λ) = 1. Hence, |V (Λ)| -J(Λ) ≤ |V ( Λ)| -t + s -1 = |V ( Λ)| -J( Λ) -t + s + 1 ≤ |V ( Λ)| -J( Λ),
which is what we wanted to show.

The next Lemma provides a lower bound for the minimum number of vertices min Λ,|Λ|=l |V (Λ)|.

Lemma 15. For ℓ ≥ 1, min Λ,|Λ|=ℓ |V (Λ)| ≥ ℓ 2 + 1 2 + ℓ 2 + 1 4 1/2 .
Proof. From Lemma 14 we can restrict ourselves to sets Λ that are connected. Let e ↓ (v) denote the edge descending from a vertex v, for which we may need to extend the boundary of the lattice. Similarly let e → (v) denote the edge extending rightwards from v. We can now define the following enlargement of the edge set Λ. The goal of our third Lemma is not just to derive a lower bound on κ ΩDIF but a lower bound that is close to optimal. By Lemma 14, κ ΩDIF (ℓ) is achieved for connected Λ, so, as with Lemma 15, we will consider such Λ's only (J(Λ) = 1). With J(Λ) = 1, the formula (12) tells us to look for the cases when |V (Λ)| is minimal in order to compute κ ΩDIF (ℓ).

What is the shape of the collection of edges Λ yielding the minimum? Recalling Euler's formula for graphs on plane: where F (Λ) is the faces of Λ which includes the 'unbounded one', we see that we are seeking Λ such that |F (Λ)| is maximal, i.e., there is maximum number of faces. By intuition, we conjecture that this happens when Λ consists of all the edges in an almost square, by which we mean V (Λ) is an r × r or r × (r + 1) rectangular grid or the inbetweens (e.g., an r×r grid of pixels to which 1 ≤ j ≤ r pixels have been added on one side). These considerations lead to the following: In fact, we observe that two edges can be added while adding one additional vertex until Λ corresponds to r × (r + 1) rectangle. Summarizing all these, a graph Λ that is constructed as above, is contained r × (r + 1) rectangle (included), and contains r×r square; satisfies either |Λ| = 2(r 2 -r)+2j or |Λ| = 2(r 2 -r)+2j+1, and |V (Λ)| = r 2 + j + 1, for j = 1, . . . , r -1-this holds for j = r as well. (Here, the case |Λ| = 2(r 2 -r) + 1 is not stated.) By a similar observation, we observe that a graph Λ that is constructed similarly as above, is contained in (r + 1) × (r + 1) square (included), and contains r × (r + 1) square; satisfies either |Λ| = 2r 2 -1 + 2j or |Λ| = 2r 2 -1 + 2j + 1, and |V (Λ)| = r 2 + r + j + 1, for j = 1, . . . , r-this holds for j = r + 1 as well. Of course, in all cases, J(Λ) = 1.

The above observation leads to the following inequalities-which we conjecture to be in fact equalities: α(2(r 2 -r) + 2j) ≤ r 2 + j, j = 1, . . . , r, α(2(r 2 -r) + 2j + 1) ≤ r 2 + j, j = 1, . . . , r, α(2r 2 -1 + 2j) ≤ r 2 + r + j, j = 1, . . . , r + 1, α(2r 2 -1 + 2j + 1) ≤ r 2 + r + j, j = 1, . . . , r + 1.

We will now express these in a simpler form in terms of |Λ| = ℓ. In the first case, letting ℓ = 2(r 2 -r) + 2j, we have

r 2 + j = ℓ 2 + r.
Since 2(r 2 -2r + 1) ≤ 2(r 2 -r + 1) ≤ ℓ ≤ 2r 2 , we have r -1 ≤ ℓ 2 ≤ r. Hence, we can write α(ℓ) ≤ ℓ 2 + ℓ 2 + 1. The other sense that we can derive a condition-which is not a guarantee-that reflects empirical results more faithfully. The idea is, instead of obtaining an upper bound of the left hand side of (25) by disregarding (or considering the worst case of) sign patterns, to model the effects of the sign patterns by estimating the size of the left hand side in terms of the maximum ℓ 2 -norm of the rows of NΩ T Λ † NΩ T Λ c (up to some constants). Though further investigation is desirable, we have empirically observed that the condition derived in this way better reflected the success rates of GAP and ℓ 1 -minimization.

Desirable properties for Ω and M. At this point, one may ask a practical tion: what are desirable properties of Ω and M that would help the performance of GAP or ℓ 1 -minimization? Can we gain some insights from our theoretical result? For this, we look for scenarios where the entries of R 0 := Ω Λ N T NΩ T

Λ Ω Λ N T -1 NΩ T Λ c are small (hence, it is likely that condition (26) is satisfied). We start with the inner expression NΩ T

Λ Ω Λ N T -1 . The larger the minimum singular value of NΩ T

Λ Ω Λ N T , the smaller the entries of R 0 . First, assuming that the rows of Ω are normalized, we note that the minimum singular value is larger when the size Λ is larger. Second, the closer the minimum singular value is to the maximum one (this is in some sense an RIP-like condition for Ω), the larger it is. These two observations tell us that Ω should have high linear dependencies (to allow large cosupport Λ) and the rows of Ω should be close to uniformly distributed on S d-1 .

Suppose that Ω has the properties described above. Then, R 0 is well approximated by R 1 := γΩ Λ N T NΩ T Λ c for some γ > 0. Therefore, we ask when the entries of Ω Λ N T NΩ T Λ c are small. Each entry of Ω Λ N T NΩ T Λ c can be guaranteed to be small if a) N satisfies an RIP condition for the space spanned by two rows of Ω and the rows of Ω are incoherent. In summary, it is desirable that:

• The rows of Ω are close to uniformly distributed in S d-1 .

• Ω is highly redundant and has highly linearly dependent structure.

• M is 'independent' from Ω. This has to do with the RIP-like properties.

• The rows of Ω are incoherent.

• The cosparsity ℓ is large.

Remark 17. The 2D finite difference operator Ω DIF may be considered incoherent even though the coherence is relatively large (1/4). This is because the majority of pairs of rows of Ω DIF are in fact uncorrelated. Heuristic comparison of success guarantees for analysis-ℓ 1 and GAP. We point out that one can obtain from (27) a condition for the GAP that is similar to [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF]. For this, we observe from [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming[END_REF] that

Ω Λ x0 ∞ = (NΩ T Λ ) † NΩ T Λ c Ω Λ c x0 ∞ = [Ω Λ c N T (Ω Λ N T ) † ] T Ω Λ c x0 ∞ .

Figure 1 :

 1 Figure 1: Number of subspaces of a given dimension, for n = p = 2d. The solid blue curve shows the log number of subspaces for the synthesis model as the dimension of subspaces vary, while the dashed red curve shows that for the analysis model.

Figure 2 :

 2 Figure 2: Top: a piecewise quadratic signal. Bottom: the support set (white region) for the wavelet coefficients of the signal using a J = 4 level shift invariant Daubechies wavelet transform with s = 3 vanishing moments. Scaling coefficients are not shown. The support set contains 122 coefficients out of a possible 512, yet the analysis subspace has a dimension of only 3.

Figure 3 :

 3 Figure 3: Top left: an example of a piecewise constant image: the 256 × 256 Shepp Logan phantom; Top right: an image with the same cosparsity, ℓ = 128014, but whose cosupport is associated with an empirically maximum subspace dimension. Bottom: zoom on the top of the top right image.

Figure 4 :

 4 Figure 4: A schematic overview of analysis cosparse vs synthesis sparse models in relation with compressed sensing.

Ω Λk x 2 2

 2 |αj|}, -Update Support: Λk = Λk-1 \ Γ k , and -Update Solution: xk = arg min x subject to y = Mx. -Stopping Criterion: If k ≥ pd + m (or k ≥ pℓ), stop. • Output: The proposed solution is xGAP = xk obtained after k iterations.

Figure 5 :

 5 Figure 5: Greedy Analysis Pursuit Algorithm (GAP)

Figure 6 :

 6 Figure 6: Recovery Rate of Analysis Algorithms for d = 200.

Figure 7 :

 7 Figure 7: Recovery Rate of Analysis Algorithms for d = 200 and σ = 1 when Ωx 0 follows the Rademacher distribution. GAP (left) and L1 (right).

Remark 12 .

 12 Due to the large size of these problems, GAP algorithm as described in Section 4 had to be modified: We used numerical optimization (the

Figure 8 :

 8 Figure 8: Recovery of 256 × 256 Shepp Logan phantom image. (a) Original Image. (b) Sampling locations of Fourier coefficients. (c) Locations where one-step difference of the original image is non-zero. Upper half corresponds to the horizontal differences and lower half the vertical differences. (d) Locations that GAP identified/eliminated to be the ones where the differences are likely non-zero. Perfect reconstruction is implied by the fact that this image 'contains' image (c).

Figure 9 :

 9 Figure 9: SNR vs the number of radial observation lines in 256 × 256 Shepp Logan phantom image recovery. The output line for the GAP is clipped due to high SNR value. SNR was computed as 20 log 10 "

+ ℓ 2 + 1 4 1 / 2

 112 Λ = {e ↓ (v), e → (v) : v ∈ V (Λ)}. (C.1)Since each edge can only descend or extend rightwards from a single vertex we have | Λ| = 2|V (Λ)|. We also have Λ ⊂ Λ.We now wish to estimate how much larger | Λ| is to |Λ|. Let us define the width, w, of V (Λ) as the number of columns spanned by V (Λ). Similarly define the height, h, as the number of rows spanned by V (Λ). For every row spanned by V (Λ) the edge extending rightwards from the right-most vertex is in Λ\Λ. Similarly for every column spanned by by V (Λ) the edge descending from the lowest vertex is also in Λ\Λ. Hence we have the following bound2|V (Λ)| = | Λ| ≥ |Λ| + w + h. (C.2)Intuitively to minimize |V (Λ)| we should choose a set of vertices with maximal area to perimeter ratio. Given that V (Λ) lies in a rectangle of size h × w we can bound the number of edges in Λ using a counting argument to obtain:ℓ ≤ h(w -1) + w(h -1) (C.3)Substitiuting this into (C.2) we get:2|V (Λ)| ≥ ℓ + w + ℓ + w 2w -1 (C.4)and similarly we can now bound the minimum possible number of edges by: min right hand side we are minimizing over all w ≥ 1. Although this includes non-integer values of w this does not invalidate the bound. It only makes it less tight.The right hand side of (C.5) is convex over w ≥ 1 with the minimum occurring at w ⋆ = 1 2 for which we also have:w ⋆ = (ℓ + w ⋆ ) (2w ⋆ -1) ,i.e. square cosupports are optimal. Inserting this into (C.5) then gives: min

  |V (Λ)| -|Λ| + |F (Λ)| = 2, (C.7)

for ℓ ≥ 5 .

 5 Proof. For r ≥ 2, we consider a subgraph corresponding to an r ×r square (solid lines) and consider graphs obtained by adding additional edges in the fashion depicted in Figure C.10.

Figure C. 10 :

 10 Figure C.10: Add dashed edges (from longer to shorter dashed) to r × r square subgraph (solid lines).

Put differently, we assume that the spark of the matrix Ω T is full, implying that every set of d rows from Ω are linearly independent.

Note that the notation Ω Λ refers to restricting rows from Ω indexed by Λ, whereas in the synthesis case we have taken the columns. We shall use this convention throughout this paper, where from the context it should be clear whether rows or columns are extracted.

When measuring sparsity with an ℓp norm, 0 < p ≤ 1, rather than with p = 0, it has been shown[START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF] that for so-called localized frames the analysis coefficients Ωx obtained with Ω = D † the canonical dual frame of D are near optimally sparse: Ωx p ≤ Cp min z|Dz=x z p, where the constant Cp does not depend on x.

Note that the same lack of sharpness holds true for our results based on (26), yet we will see that these can actually provide cosparse signal recovery guarantees in simple but nontrivial cases.

One could also construct Ω by simply drawing the rows of it randomly and independently from S d-1 without the tight frame constraint. We have run the experiment for such operators and observed that the result was similar.

The code for l1magic was downloaded from http://www.acm.caltech.edu/l1magic/ and the one for AIHT from http://www.personal.soton.ac.uk/tb1m08/sparsify/AIHT Paper Code.zip. The result for the back-projection was obtained using the code for AIHT.

This work was supported in part by the EU through the project SMALL (Sparse Models, Algorithms and Learning for Large-Scale data), FET-Open programme, under grant number: 225913

To obtain Corollary 8, observe that we can remove the constraint z ∈ Null(M) by writing z = N T α where N T is an d × (d -m) basis matrix for Null(M) and α ∈ R d-m is an appropriate coefficient sequence. Thus, the necessary and sufficient condition becomes

Since the ℓ × (d -m) matrix Ω Λ N T is thin (ℓ ≥ d -m) and full-rank, defining β := Ω Λ N T α, we have α = (Ω Λ N T ) † β. Therefore, a sufficient (but no longer necessary) recovery condition for analysis ℓ 1 -minimization is

Equivalently, for all x Λ with Ω Λ x Λ = 0, sup

that is to say

Condition [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF] follows from the above. To conclude the proof of Corollary 8, we note that since sign(Ω Λ c x Λ ) ∞ ≤ 1, the left hand side of (A.4) is bounded above by [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] implies [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF] and the proof is complete.

Appendix B. Proof of Lemma 9

Since x0 is the solution of arg min x Ωx 2 2 subject to y = Mx, applying the Lagrange multiplier method, we observe that x0 satisfies

three cases can be treated similarly and we obtain

Therefore, for all ℓ ≥ 5, we have α(ℓ) ≤ ℓ 2 + ℓ 2 + 1.

We now put these ingredients together.

Proof of Proposition 6. The proof of the lower bound comes directly from Lemma 16.

To prove the upper bound we note that from Lemma 14, Lemma 15 and Eq. ( 12) we have:

as required.

Appendix D. Discussion on the analysis exact recovery condition

We observe that the analysis ERC condition [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] is not sharp in general, especially for the redundant Ω. In the case of GAP, tracing the arguments of Lemma 9 and Theorem 10, we conclude that in order for [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] to be sharp, there must exist a cosparse signal x 0 such that, with x0 defined as in [START_REF] Dragotti | Wavelet footprints: theory, algorithms, and applications[END_REF], Ω Λ c x0 matches the exact sign pattern of the row of (NΩ T Λ ) † NΩ T Λ c with the largest ℓ 1 -norm and is of constant magnitude in absolute value. We remind that x0 is the initial estimate that appears in the algorithm. Since the collection of Ω Λ c x0 may not span the whole R Λ c , especially when Ω is over-complete, it is unreasonable to expect the existence of such an x 0 . Similarly, in the case of analysis ℓ 1 , we know that ( 26) is obtained from [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF] in a crude way without taking into account the sign patterns of Ω Λ c x Λ , which is not sharp in general for redundant Ω.

Average case performance guarantees? Can we think of a way to obtain a more realistic success guarantee? We have a partial answer for this question in the Since Ω Λ x0 ∞ < Ω Λ c x0 ∞ is the necessary and sufficient condition for the (one-step) success of the GAP, we can derive a necessary and sufficient condition:

where x 0 is varied over all signals with cosupport Λ and x0 is the signal resulting from the first step of GAP. The above condition can be rewritten in a form similar to [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF]:

where x0 is derived as in [START_REF] Dragotti | Wavelet footprints: theory, algorithms, and applications[END_REF], ⊙ denotes the element-wise multiplication of vectors, and v is obtained from Ω Λ c x0 by taking element-wise absolute values and normalizing it to a unit ℓ ∞ -norm (v :

and ( 25) are in a similar form, but there are two differences between the two: First, for (D.1), the signal x0 that apears is not in general a vector with cosupport Λ. It is rather a signal that arises from an approximation. Second, there is a 'weight' vector v in (D.1). One can heuristically deduce that such a v favors condition (D.1) to hold true since the size of most entries of v likely be smaller than 1. Beside these differences, one should keep in mind that condition (D.1) is only for one step.

Appendix E. Derivation of the Solution (19)

By the method of Lagrange multiplier, x0 is the solution of (18

We first assume that M is of full rank, thus MM T is invertible. With this assumption, we can solve for µ and obtain µ = (MM T ) -1 MΩ T Ωx 0 . Substituting this in and combining the two equation, we have (Id -M T (MM T ) -1 M)Ω T Ω M x0 = 0 y . Now we assume that the matrix on the right hand side is of full rank. Multiplying the both sides by (Id -M T (MM T ) -1 M) M T (which is of full rank), we arrive at (M T M + (Id -M T (MM T ) -1 M)Ω T Ω)x 0 = M T y, from which [START_REF] Elad | Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing[END_REF] follows. Above, we have used the fact that Id-M T (MM T ) -1 M is a projection.