Samuel Vaiter
email: samuel.vaiter@ceremade.dauphine.fr

Gabriel Peyré
email: gabriel.peyre@ceremade.dauphine.fr

Jalal Fadili
email: jalal.fadili@greyc.ensicaen.fr

Low Complexity Regularization of Linear Inverse Problems

Inverse problems and regularization theory is a central theme in imaging sciences, statistics and machine learning. The goal is to reconstruct an unknown vector from partial indirect, and possibly noisy, measurements of it. A now standard method for recovering the unknown vector is to solve a convex optimization problem that enforces some prior knowledge about its structure. This chapter delivers a review of recent advances in the field where the regularization prior promotes solutions conforming to some notion of simplicity/low-complexity. These priors encompass as popular examples sparsity and group sparsity (to capture the compressibility of natural signals and images), total variation and analysis sparsity (to promote piecewise regularity), and low-rank (as natural extension of sparsity to matrix-valued data). Our aim is to provide a unified treatment of all these regularizations under a single umbrella, namely the theory of partial smoothness. This framework is very general and accommodates all low-complexity regularizers just mentioned, as well as many others. Partial smoothness turns out to be the canonical way to encode low-dimensional models that can be linear spaces or more general smooth manifolds. This review is intended to serve as a one stop shop toward the understanding of the theoretical properties of the so-regularized solutions. It covers a large spectrum including: (i) recovery guarantees and stability to noise, both in terms of ℓ 2 -stability and model (manifold) identification; (ii) sensitivity analysis to perturbations of the parameters involved (in particular the observations), with applications to unbiased risk estimation ; (iii) convergence properties of the forward-backward proximal splitting scheme, that is particularly well suited to solve the corresponding large-scale regularized optimization problem.

Inverse Problems and Regularization

In this chapter, we deal with finite-dimensional linear inverse problems.

Forward Model

Let x 0 ∈ R N be the unknown vector of interest. Suppose that we observe a vector y ∈ R P of P linear measurements according to y = Φx 0 + w, [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF] where w ∈ R P is a vector of unknown errors contaminating the observations. The forward model ( 1) offers a model for data acquisition that describes a wide range of problems in data processing, including signal and image processing, statistics, and machine learning. The linear operator Φ : R N → R P , assumed to be known, is typically an idealization of the acquisition hardware in imaging science applications, or the design matrix in a parametric statistical regression problem. The noise w can be either deterministic (in this case, one typically assumes to know some bound on its ℓ 2 norm w ), or random (in which case its distribution is assumed to be known). Except in Sections 4.4 and 5.3 where the noise is explicitly assumed random, w is deterministic throughout the rest of the chapter. We refer to [191] and [START_REF] Boncelet | Handbook of Image and Video Processing, chapter Image Noise Models[END_REF] for a comprehensive account on noise models in imaging systems.

Solving an inverse problem amounts to recovering x 0 , to a good approximation, knowing y and Φ according to [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF]. Unfortunately, the number of measurements P can be much smaller than the ambient dimension N of the signal. Even when P = N, the mapping Φ is in general ill-conditioned or even singular. This entails that the inverse problem is in general ill-posed. In signal or image processing, one might for instance think of Φ as a convolution with the camera point-spread function, or a subsampling accounting for low-resolution or damaged sensors. In medical imaging, typical operators represent a (possibly subsampled) Radon transform (for computerized tomography), a partial Fourier transform (for magnetic resonance imaging), a propagation of the voltage/magnetic field from teh dipoles to the sensors (for electro-or magneto-encephalography). In seismic imaging, the action of Φ amounts to a convolution with a wavelet-like impulse response that approximates the solution of a wave propagation equation in media with discontinuities. For regression problems in statistics and machine learning, Φ is the design matrix whose columns are P covariate vectors.

Variational Regularization

As argued above, solving an inverse problem from the observations [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF] is in general ill-posed. In order to reach the land of well-posedness, it is necessary to restrict the inversion process to a well-chosen subset of R N containing the plausible solutions including x 0 ; e.g. a linear space or a union of subspaces. A closely related procedure, that we describe next, amounts to adopting a variational framework where the sought-after solutions are those where a prior penalty/regularization function is the smallest. Though this approach may have a maximum a posteriori Bayesian interpretation, where a random prior is placed on x 0 , this is not the only interpretation. In fact, we put no randomness whatsoever on the class of signals we look for. We will not elaborate more on these differences in this chapter, but the reader may refer to [119] for an insightful discussion.

The foundations of regularization theory can be traced back to the pioneering work of the Russian school, and in particular of Tikhonov in 1943 when he proposed the notion of conditional well-posedness. In 1963, Tikhonov [217, 218] introduced what is now commonly referred to as Tikhonov (or also Tikhonov-Phillips) regularization, see also the book [219]. This corresponds, for λ > 0, to solving an optimization problem of the form

x ⋆ ∈ Argmin x∈R N 1 2λ
Φxy 2 + J(x).

(P y,λ )

Data fidelity

In (P y,λ ), Φxy 2 stands for the data fidelity term. If the noise happens to be random, then using a likelihood argument, an appropriate fidelity term conforming to the noise distribution can be used instead of the quadratic data fidelity. Clearly, it is sufficient then to replace the latter by the negative log-likelihood of the distribution underlying the noise. Think for instance of the Csiszár's I-divergence for Poisson noise. We would also like to stress that many of the results provided in this chapter extend readily when the quadratic loss in the fidelity term, i.e. µ → y -µ 2 , is replaced by any smooth and strongly convex function, see in particular Remark 13. To make our exposition concrete and digestible, we focus in the sequel on the quadratic loss.

Regularization

The function J : R N → R is the regularization term which is intended to promote some prior on the vector to recover. We will consider throughout this chapter that J is a convex finite-valued function. Convexity plays an important role at many locations, both on the recovery guarantees and the algorithmic part. See for instance Section 6 which gives a brief overview of recent algorithms that are able to tackle this class of convex optimization problems. It is however important to realize that non-convex regularizing penalties, as well as non-variational methods (e.g. greedy algorithms), are routinely used for many problems such as sparse or low-rank recovery. They may even outperform in practice their convex counterparts/relaxation. It is however beyond the scope of this chapter to describe these algorithms and the associated theoretical performance guarantees. We refer to Section 2.1 for a brief account on non-convex model selection approaches. The scalar λ > 0 is the regularization parameter. It balances the trade-off between fidelity and regularization. Intuitively, and anticipating on our theoretical results hereafter, this parameter should be adapted to the noise level w and the known properties of the vector x 0 to recover. Selecting optimally and automatically λ for a given problem is however difficult in general. This is at the hear of Section 5, where unbiased risk estimation strategies are shown to offer a versatile solution.

Note that since Φ is generally not injective and J is not coercive, the objective function of (P y,λ ) is neither coercive nor strictly convex. In turn, there might be existence (of minimizers) issues, and even if minimizers exist, there are not unique in general.

Under mild assumptions, problem (P y,λ ) is formally equivalent to the constrained formulations min {J(x) ; y -Φx ε} , (P 1 y,ε ) min { y -Φx ; J(x) γ} , (P 2 y,γ ) in the sense that there exists a bijection between each pair of parameters among (λ , ε, γ) so that the corresponding problems share the same set of solutions. However, this bijection is not explicit and depends on y, so that both from an algorithmic point of view and a theoretical one, each problem may need to be addressed separately. See the recent paper [START_REF] Ciak | Homogeneous penalizers and constraints in convex image restoration[END_REF] and references therein for a detailed discussion, and [154, Theorem 2.3] valid also in the non-convex case. We focus in this chapter on the penalized/Tikhonov formulation (P y,λ ), though most of the results stated can be extended to deal with the constrained ones (P 1 y,ε ) and (P 2 y,γ ) (the former is known as the residual method or Mozorov regularization and the latter as Ivanov regularization in the inverse problems literature).

The value of λ should typically be an increasing function of w . In the special case where there is no noise, i.e. w = 0, the fidelity to data should be perfect, which corresponds to considering the limit of (P y,λ ) as λ → 0 + . Thus, assuming that y ∈ Im(Φ), as is the case when w = 0, it can be proved that the solutions of (P y,λ ) converge to the solutions of the following constrained problem [217,198] x ⋆ ∈ Argmin x∈R N J(x) subject to Φx = y.

(P y,0 )

Notations

For any subspace T of R N , we denote P T the orthogonal projection onto T , x T = P T (x) and Φ T = Φ P T . For a matrix A, we denote A * its transpose, and A + its Moore-Penrose pseudo-inverse. For a convex set E, aff(E) denotes its affine hull (i.e. the smallest affine space containing it), and lin(E) its linear hull (i.e. the linear space parallel to aff(E)). Its relative interior ri(E) is the interior for the topology of aff(E) and rbd(E) is its relative boundary. For a manifold M , we denote T M (x) the tangent space of M at x ∈ M . A good source on smooth manifold theory is [144].

A function J : R N → R ∪ {+∞} is said to be proper if it is not identically +∞. It is said to be finite-valued if J(x) ∈ R for all x ∈ R N . We denote dom(J) the set of points x where

J(x) ∈ R is finite. J is said to be closed if its epigraph {(x, y) ; J(x) y} is closed. For a set C ⊂ R N , the indicator function ι C is defined as ι C (x) = 0 if x ∈ C and ι C (x) = +∞ otherwise.
We recall that the subdifferential at x of a proper and closed convex function

J : R N → R ∪ {+∞} is the set ∂ J(x) = η ∈ R N ; ∀ δ ∈ R N , J(x + δ ) J(x) + η, δ .
Geometrically, when J is finite at x, ∂ J(x) is the set of normals to the hyper-planes supporting the graph of J and tangent to it at x. Thus, ∂ J(x) is a closed convex set. It is moreover bounded, hence compact, if and only if x ∈ int(dom(J)). The size of the subdifferential at x ∈ dom(J) reflects in some sense the degree of nonsmoothness of J at x. The larger the subdifferential at x, the larger the "kink" of the graph of J at x. In particular, if J is differentiable at x, then ∂ J(x) is a singleton and ∂ J(x) = {∇J(x)}.

As an illustrative example, the subdifferential of the absolute value is

∀ x ∈ R, ∂ | • |(x) = sign(x) if x = 0, [-1, 1] otherwise. ( 2
)
The ℓ 1 norm

∀ x ∈ R N , x 1 = N ∑ i=1 |x i |
is a popular low-complexity prior (see Section 2.3.1 for more details). Formula ( 2) is extended by separability to obtain the subdifferential of the ℓ 1 norm

∂ • 1 (x) = η ∈ R N ; η ∞ 1 and ∀ i ∈ I, sign(η i ) = sign(x i ) (3)
where

I = supp(x) = {i ; x i = 0}. Note that at a point x ∈ R N such that x i = 0 for all i, • 1 is differentiable, and ∂ • 1 (x) = {sign(x)}.

Low Complexity Priors

A recent trend in signal and image processing, statistics and machine learning is to make use of large collections of so-called "models" to account for the complicated structures of the data to handle. Generally speaking, these are manifolds M (most of the time linear subspaces), and hopefully of low complexity (to be detailed later), that capture the properties of the sought after signal, image or higher dimensional data. In order to tractably manipulate these collections, the key idea underlying this approach is to encode these manifolds in the non-smooth parts of the regularizer J.

As we detail here, the theory of partial smoothness turns out to be natural to provide a mathematically grounded and unified description of these regularizing functions.

Model Selection

The general idea is thus to describe the data to recover using a large collection of models M = {M } M ∈M , which are manifolds. The "complexity" of elements in such a manifold M is measured through a penalty pen(M ). A typical example is simply the dimensionality of M , and it should reflect the intuitive notion of the number of parameters underlying the description of the vector x 0 ∈ M that one aims at recovering from the noisy measurements of the form [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF]. As popular examples of such low complexity, one thinks of sparsity, piecewise regularity, or low rank. Penalizing in accordance to some notion of complexity is a key idea, whose roots can be traced back to the statistical and information theory literature, see for instance [162,[START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]. Within this setting, the inverse problem associated to the measurements ( 1) is solved by restricting the inversion to an optimal manifold as selected by pen(M ). Formally, this would correspond to solving (P y,λ ) with the combinatorial regularizer

J(x) = inf {pen(M ) ; M ∈ M and x ∈ M } . (4)
A typical example of such a model selection framework is that with sparse signals, where the collection M corresponds to a union of subspaces, each of the form

M = x ∈ R N ; supp(x) ⊆ I .
Here I ⊆ {1, . . ., N} indexes the supports of signals in M , and can be arbitrary. In this case, one uses pen(M ) = dim(M ) = |I|, so that the associated combinatorial penalty is the so-called ℓ 0 pseudo-norm

J(x) = x 0 = | supp(x)| = | {i ∈ {1, . . ., N} ; x i = 0} |. (5)
Thus, solving (P y,λ ) is intended to select a few active variables (corresponding to non-zero coefficients) in the recovered vector.

These sparse models can be extended in many ways. For instance, piecewise regular signals or images can be modeled using manifolds M that are parameterized by the locations of the singularities and some low-order polynomial between these singularities. The dimension of M thus grows with the number of singularities, hence the complexity of the model.

Literature review. The model selection literature [START_REF] Birgé | From model selection to adaptive estimation[END_REF][START_REF] Barron | Risk bounds for model selection via penalization[END_REF][START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] proposes many theoretical results to quantify the performance of these approaches. However, a major bottleneck of this class of methods is that the corresponding J function defined in ( 4) is non-convex, and even not necessarily closed, thus typically leading to highly intractable combinatorial optimization problems. For instance, in the case of the ℓ 0 penalty [START_REF] Asif | Sparse recovery of streaming signals using L1-homotopy[END_REF] and for an arbitrary operator Φ, (P y,λ ) is known to be NP-hard, see e.g. [168].

It then appears crucial to propose alternative strategies which allow to deploy fast computational algorithms. A first line of work consists in finding stationary points of (P y,λ ) using descent-like schemes. For instance, in the case of the ℓ 0 pseudonorm, this can be achieved using iterative hard thresholding [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF]124], or iterative reweighting schemes which consist of solving a sequence of weighted ℓ 1 -or ℓ 2minimization problems where the weights used for the next iteration are computed from the values of the current solution, see for instance [189,[START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF][START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF] and references therein. Another class of approaches is that of greedy algorithms. These are algorithms which explore the set of possible manifolds M by progressively, actually in a greedy fashion, increasing the value of pen(M ). The most popular schemes are matching pursuit [161] and its orthogonal variant [180,[START_REF] Davis | Adaptive time-frequency approximations with matching pursuits[END_REF], see also the comprehensive review [169] and references therein. The last line of research, which is the backbone of this chapter, consists in considering convex regularizers which are built in such away that they promote the same set of low-complexity manifolds M. In some cases, the convex regularizer proves to be the convex hull of the initial (restricted) non-convex combinatorial penalty [START_REF] Amelunxen | Living on the edge: A geometric theory of phase transitions in convex optimization[END_REF]. But these convex penalties can also be designed without being necessarily convexified surrogates of the original non-convex ones.

In the remainder of this section, we describe in detail a general framework that allows model selection through the general class of convex partly smooth functions.

Encoding Models into Partly Smooth Functions

Before giving the precise definition of our class of convex priors, we define formally the subspace T x .

Definition 1 (Model tangent subspace). For any vector x ∈ R N , we define the model tangent subspace of x associated to J

T x = lin(∂ J(x)) ⊥ ,
In fact, the terminology "tangent" originates from the sharpness property of Definition 2(ii) below, when x belongs to the manifold M .

When J is differentiable at x, i.e. ∂ J(x) = {∇J(x)}, one has T x = R N . On the contrary, when J is not smooth at x, the dimension of T x is of a strictly smaller dimension, and J essentially promotes elements living on or close to the affine space x + T x .

We can illustrate this using the ℓ 1 norm J = • 1 defined in [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]. Using formula (3) for the subdifferential, one obtains that

T x = u ∈ R N ; supp(u) ⊆ supp(x) ,
which is the set of vector having the same sparsity pattern as x.

Toward the goal of studying the recovery guarantees of problem ( 4), our central assumption is that J is a partly smooth function relative to some manifold M . Partial smoothness of functions was originally defined [146]. Loosely speaking, a partly smooth function behaves smoothly as we move on the manifold M , and sharply if we move normal to it. Our definition hereafter specializes that of [146] to the case of finite-valued convex functions. Definition 2. Let J be a finite-valued convex function. J is partly smooth at x relative to a set M containing x if

(i) (Smoothness) M is a C 2 -manifold around x and J restricted to M is C 2 around x. (ii) (Sharpness) The tangent space T M (x) is T x . (iii) (Continuity) The set-valued mapping ∂ J is continuous at x relative to M .
J is said to be partly smooth relative to a set M if M is a manifold and J is partly smooth at each point x ∈ M relative to M . J is said to be locally partly smooth at x relative to a set M if M is a manifold and there exists a neighbourhood U of x such that J is partly smooth at each point of M ∩U relative to M .

Remark 1 (Uniqueness of M ). In the previous definition, M needs only to be defined locally around x, and it can be shown to be locally unique, see [130,Corollary 4.2]. In the following we will thus often denote M x any such a manifold for which J is partly smooth at x. Taking once again the example of J = • 1 , one sees that in this case, M x = T x because this function is polyhedral. Section 2.3.6 bellow defines fonctions J for which M x differs in general from T x .

Examples of Partly Smooth Regularizers

We describe below some popular examples of partly smooth regularizers that are widely used in signal and image processing, statistics and machine learning. We first expose basic building blocks (sparsity, group sparsity, anti-sparsity) and then show how the machinery of partial smoothness enables a powerful calculus to create new priors (using pre-and post-composition, spectral lifting, and positive linear combinations).

ℓ 1 Sparsity

One of the most popular non-quadratic convex regularization is the ℓ 1 norm

J(x) = x 1 = N ∑ i=1 |x i |,
which promotes sparsity. Indeed, it is easy to check that J is partly smooth at x relative to the subspace

M x = T x = u ∈ R N ; supp(u) ⊆ supp(x) .
Another equivalent way to interpret this ℓ 1 prior is that it is the convex enveloppe (restricted to the ℓ 2 -ball) of the ℓ 0 pseudo-norm [START_REF] Asif | Sparse recovery of streaming signals using L1-homotopy[END_REF], in the sense that the ℓ 1 -unit ball is the convex hull of the restriction of the unit ball of the ℓ 0 -pseudo norm to the ℓ 2 -unit ball.

Literature review. The use of the ℓ 1 norm as a sparsity-promoting regularizer traces back several decades. An early application was deconvolution in seismology [START_REF] Claerbout | Robust modeling with erratic data[END_REF]197,212]. Rigorous recovery results began to appear in the late 1980's [START_REF] Donoho | Uncertainty principles and signal recovery[END_REF][START_REF] Donoho | Signal recovery and the large sieve[END_REF]. In the mid-1990's, ℓ 1 regularization of least-square problems has been popularized in the signal processing literature under the name basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] and in the statistics literature under the name Lasso [213]. Since then, the applications and understanding of ℓ 1 minimization have continued to increase dramatically.

ℓ 1 -ℓ 2 Group Sparsity

To better capture the sparsity pattern of natural signals and images, it is useful to structure the sparsity into non-overlapping groups B such that b∈B b = {1, . . . , N}. This group structure is enforced by using typically the mixed ℓ 1 -ℓ 2 norm

J(x) = x 1,B = ∑ b∈B x b , (6)
where

x b = (x i ) i∈b ∈ R |b| .
Unlike the ℓ 1 norm, and except the case |b| = 1 for all b ∈ B, the ℓ 1 -ℓ 2 norm is not polyhedral, but is still partly smooth at x relative to the linear manifold

M x = T x = {u ; supp B (u) ⊆ supp B (x)} where supp B (x) = {b ; x b = 0} .
Literature review. The idea of group/block sparsity has been first proposed by [127,126,[START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF] for wavelet block shrinkage, i.e. when Φ = Id. For over-determined regression problems of the form [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF], it has been introduced by [START_REF] Bakin | Adaptive regression and model selection in data mining problems[END_REF]244]. Group sparsity has also been extensively used in machine learning in e.g. [START_REF] Bach | Consistency of the group Lasso and multiple kernel learning[END_REF] (regression and mutiple kernel learning) and [175] (for multi-task learning).The wavelet coefficients of a natural image typical exhibit some group structure, see [160] and references therein on natural image modeling. Indeed, edges and textures induce strong dependencies between coefficients. In audio processing, it has proved useful to structure sparsity in multi-channel data [121]. Group sparsity is also at the heart of the socalled multiple measurements vector (MMV) model, see for instance [START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF][START_REF] Chen | Theoretical results on sparse representations of multiple-measurement vectors[END_REF]. It is possible to replace the ℓ 2 norm with more general functionals, such as ℓ p norms for p > 1, see for instance [225,171,237].

2.3.3 ℓ ∞ Anti-sparsity
In some cases, the vector to be reconstructed is expected to be flat. Such a prior can be captured using the ℓ ∞ norm

J(x) = x ∞ = max i∈{1,...,n} |x i |.
It can be readily checked that this regularizer is partly smooth (in fact polyhedral) relative to the subspace

M x = T x = {u ; u I = ρx I for some ρ ∈ R} , where I = {i ; x i = x ∞ } .
Literature review.

The ℓ ∞ regularization has found applications in computer vision, such as for database image retrieval [137]. For this application, it is indeed useful to have a compact signature of a signal x, ideally with only two values ± x ∞ (thus achieving optimal anti-sparsity since dim(T x ) = 1 in such a case). An approach proposed in [138] for realizing this binary quantification is to compute these vectors as solutions of (P y,λ ) for J = • ∞ and a random Φ. A study of this regularization is done in [108], where an homotopy-like algorithm is provided. The use of this ℓ ∞ regularization is also connected to Kashin's representation [157], which is known to be useful in stabilizing the quantization error for instance. Others applications such as wireless network optimization [211] also rely on the ℓ ∞ prior.

Synthesis Regularizers

Sparsity or more general low-complexity regularizations are often used to model coefficients α ∈ R Q describing the data x = Dα in a dictionary D ∈ R N×Q of Q atoms in R N . Given a partly smooth function J 0 : R Q → R, we define the following synthesis-type prior J : R N → R as the pre-image of J 0 under the linear mapping

D J(x) = min α∈R Q J 0 (α) s.t. Dα = x
Since J 0 is bounded below and convex, J is convex. If D is surjective (as in most cases with redundant dictrionaries), then J is also finite-valued. The initial optimization (P y,λ ) can equivalently been solved directly over the coefficients domain to obtain x ⋆ = Dα ⋆ where

α ⋆ ∈ Argmin α∈R Q 1 2λ y -ΦDα 2 + J 0 (α) (7)
which can be interpreted as a regularized inversion of the operator ΦD using the prior J 0 .

It is possible to study directly the properties of the solutions α ⋆ to [START_REF] Bach | Consistency of the group Lasso and multiple kernel learning[END_REF], which involves directly partial-smoothness of J 0 . A slightly different question is to understand the behavior of the solutions x ⋆ = Dα ⋆ of (P y,λ ), which requires to study partial smoothness of J itself. In the case where D is invertible, both problems are completely equivalent.

Literature review. Sparse synthesis regularization using J 0 = • 1 is popular in signal and image processing to model natural signals and images, see for instance [160,207] for a comprehensive account. The key problem to achieve good performance in these applications is to design a dictionary to capture sparse representations of the data to process. Multiscale dictionaries built from wavelet pyramids are popular to sparsely represent transient signals with isolated singularities and natural images [159]. The curvelet transform is known to provide non-adaptive near-optimal sparse representation of piecewise smooth images away from smooth edges (so-called cartoon images) [START_REF] Candès | Curvelets: A surprisingly effective nonadaptive representation for objects with edges[END_REF]. Gabor dictionaries (made of localized and translated Fourier atoms) are popular to capture locally stationary oscillating signals for audio processing [START_REF] Allen | Short-term spectral analysis, and modification by discrete Fourier transform[END_REF]. To cope with richer and diverse contents, researchers have advocated to concatenate several dictionaries to solve difficult problems in signal and image processing, such as component separation or inpainting, see for instance [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA)[END_REF]. A line of current active research is to learn and optimize the dictionary from exemplars or even from the available data themselves. We refer to [START_REF] Elad | Sparse and redundant representations: from theory to applications in signal and image processing[END_REF]Chapter 12] for a recent overview of the relevant literature.

Analysis Regularizers

Analysis-type regularizers (following the terminology introduced in [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]) are of the form

J(x) = J 0 (D * x) ,
where D ∈ R N×Q is a linear operator. Such a prior controls the low complexity (as measured by J 0 ) of the correlations between the columns of D and the signal x. If J 0 is partly smooth at z = D * x for the manifold M 0 z , then it is shown in [146, Theorem 4.2] that J is partly smooth at x relative to the manifold

M x = u ∈ R N ; D * u ∈ M 0 z
provided that the following transversality condition holds [144,Theorem 6.30(a)]

Ker(D) ∩ T M 0 z (z) ⊥ = {0} ⇐⇒ Im(D * ) + T M 0 z (z) = R N .
Literature review. A popular example is when J 0 = • 1 and D * a finite-difference discretization of the derivative of a 1-D signal or a 2-D image. This defines the anisotropic total variation semi-norm, which promotes piecewise constant signals or images [196]. The 2-D isotropic total variation semi-norm can be interpreted as taking J 0 = • 1,2 with blocks of size two. A comprehensive review of total variation regularization can be found in [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF]. TV regularization has been extended in several ways to model piecewise polynomial functions, see in particular the Total Generalized Variation prior [START_REF] Bredies | Total generalized variation[END_REF].

One can also use a wavelet dictionary D which is shift-invariant, such that the corresponding regularization J can be seen as a kind of multi-scale total variation. This is typically the case of the Haar wavelet dictionary [208]. When using higher order wavelets, the corresponding priors favors models M composed of discrete piecewise polynomials.

The Fused Lasso [214] corresponds to J 0 being the ℓ 1 -norm and D is the concatenation of the identity and the adjoint of a finite-difference operator. The corresponding models M are composed of disjoint blocks over which the signals are constant.

Defining a block extracting operator D * x = (x b ) b∈B allows to re-write the group ℓ 1 -ℓ 2 norm ( 6), even with overlapping blocks [139,182,245,[START_REF] Cai | Incorporating information on neighbouring coefficients into wavelet estimation[END_REF]. To cope with correlated covariates in linear regression, analysis-type sparsity-enforcing prior were proposed in [118,193] using J 0 = • * the nuclear norm (as defined in Section 2.3.6).

(i.e. ∃(b, b ′ ) ∈ B 2 with b ∩ b ′ = / 0), as J = J 0 • D * where J 0 = • 1,2 without overlap, see
For unitary D, the solutions of (P y,λ ) with synthesis and analysis regularizations are obviously the same. In the general case (e.g. D overcomplete), however, these two regularizations are different. Some authors have reported results comparing these two priors for the case where J 0 is the ℓ 1 norm [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]199]. A first discussion on the relation and distinction between analysis and synthesis ℓ 1 -sparse regularizations can be found in [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. But only very recently, some theoretical recovery results and algorithmic developments on ℓ 1 -analysis regularization (so-called cosparse model) have began to be developed, see e.g. [167,230].

Spectral Functions

The natural extension of low-complexity priors to matrix-valued data x ∈ R N 0 ×N 0 (where N = N 2 0 ) is to impose the low-complexity on the singular values of the matrix. We denote

x = U x diag(Λ x )V * x an SVD decomposition of x, where Λ x ∈ R N 0 + . If j : R N 0 → R is a permutation-invariant closed convex function, then one can con- sider the function J(x) = j(Λ x )
which can be shown to be a convex function as well [147]. When restricted to the linear space of symmetric matrices, j is partly smooth at Λ x for a manifold m Λ x , if and only if J is partly smooth at x relative to the manifold

M x = U diag(Λ )U * ; Λ ∈ m Λ x ,U ∈ O N 0 ,
where O N 0 ⊂ R N 0 ×N 0 is the orthogonal group. The proof of this assertion can be found in [START_REF] Daniilidis | Orthogonal invariance and identifiability[END_REF]Theorem 3.19], which builds upon the work of [START_REF] Daniilidis | Spectral (isotropic) manifolds and their dimension[END_REF] on manifold smoothness transfer under spectral lifting. This result can be extended to nonsymmetric matrices by requiring that j is an absolutely permutation-invariant closed convex function, see [START_REF] Daniilidis | Orthogonal invariance and identifiability[END_REF]Theorem 5.3].

Literature review. The most popular spectral prior is obtained for j = • 1 . This defines the nuclear norm, or 1-Schatten norm, as

J(x) = x * = Λ x 1 . (8)
It can be shown that the nuclear norm is the convex hull of the rank function with respect to the spectral norm ball, see [102,133]. It then corresponds to promoting a low-rank prior. Moreover, the nuclear norm can be shown to be partly smooth at x relative to the set [148, Example 2]

M x = {u ; rank(u) = rank(x)}
which is a manifold around x.

The nuclear norm has been used in signal and image processing, statistics and machine learning for various applications, including low rank matrix completion [205,190,[START_REF] Candès | Exact matrix completion via convex optimization[END_REF], principal component pursuit [START_REF] Candès | Robust principal component analysis?[END_REF], model reduction [103], and phase retrieval [START_REF] Candès | Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming[END_REF]. It is also used for some imaging applications, see for instance [152].

Mixed Regularizations

Starting from a collection of convex functions {J ℓ } ℓ∈L , L = {1, . . . , L}, it is possible to design a convex function as

J ℓ (x) = ∑ ℓ∈L ρ ℓ J ℓ (x),
where ρ ℓ > 0 are weights. If each J ℓ is partly smooth at x relative to a manifold M ℓ x , then it is shown in [146, Corollary 4.8] that J is also partly smooth at x for

M x = ℓ∈L M ℓ x ,
with the proviso that the manifolds M ℓ x intersect transversally [144, Theorem 6.30(b)], i.e. the sum of their respective tangent spaces T M ℓ x (x) spans the whole ambient space R N . Literature review. A popular example is to impose both sparsity and low rank of a matrix, when using J 1 = • 1 and J 2 = • * , see for instance [114,177].

Separable Regularization

Let {J ℓ } ℓ∈L , L = {1, . . . , L}, be a family of convex functions. If J ℓ is partly smooth at x ℓ relative to a manifold M ℓ

x ℓ , then the separable function Proposition 4.5]. Literature review. One fundamental problem that has attracted a lot of interest in the recent years in data processing involves decomposing an observed object into a linear combination of components/constituents x ℓ , ℓ ∈ L = {1, . . . , L}. One instance of such a problem is image decomposition into texture and piece-wise-smooth (cartoon) parts. The corresponding forward model can be cast in the form [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF], where

J ({x ℓ } ℓ∈L ) = ∑ ℓ∈L J ℓ (x ℓ ) is partly smooth at (x 1 , . . . , x L ) relative to M 1 x 1 × • • • × M L x L [146,
x 0 = x 1 x 2
, x 1 and x 2 are the texture and cartoon components, and Φ = [Id Id]. The decomposition is then achieved by solving the variational problem (P y,λ ), where J 1 is designed to promote the discontinuities in the image, and J 2 to favor textures; see [206,[START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF]184] and references therein. Another example of decomposition is principal component pursuit, proposed in [START_REF] Candès | Robust principal component analysis?[END_REF], to decompose a matrix which is the superposition of a low-rank component and a sparse component. In this case J 1 = • 1 and

J 2 = • * .

ℓ 2 Stability

In this section, we assume that J is a finite-valued convex function, but it is not assumed to be partly smooth.

The observations y are in general contaminated by noise, as described by the forward model [START_REF] Agarwal | Exact recovery of sparsely used overcomplete dictionaries[END_REF]. It is thus important to study the ability of (P y,λ ) to recover x 0 to a good approximation in presence of such a noise w, and to assess how the reconstruction error decays as a function of the noise level. In this section, we present a generic result ensuring a so-called "linear convergence rate" in terms of ℓ 2 -error between a recovered vector and x 0 (see Theorem 1), which encompasses a large body of literature from the inverse problems community.

Dual Certificates

It is intuitively expected that if (P y,λ ) is good at recovering an approximation of x 0 in presence of noise, then (P y,0 ) should be able to identify x 0 uniquely when the noise vanishes, i.e. y = Φx 0 . For this to happen, the solution to (P y,0 ) has to satisfy some non-degeneracy condition. To formalize this, we first introduce the notion of dual certificate.

Definition 3 (Dual certificates). For any vector x ∈ R N , the set of dual certificates at x is defined as

D(x) = Im(Φ * ) ∩ ∂ J(x) .
The terminology "dual certificate" was introduced in [START_REF] Candès | Exact matrix completion via convex optimization[END_REF]. One can show that the image by Φ * of the set of solutions of the Fenchel-Rockafellar dual to (P y,0 ) is precisely D(x).

It is also worth noting that x 0 being a solution of (P y,0 ) for y = Φx 0 is equivalent to D(x 0 ) = / 0. Indeed, this is simply a convenient re-writing of the first order optimality condition for (P y,0 ).

To ensure stability of the set of minimizers (P y,λ ) to noise perturbing the observations Φx 0 , one needs to introduce the additional requirement that the dual certificates should be strictly inside the subdifferential of J at x 0 . This is precisely the non-degeneracy condition mentioned previously.

Definition 4 (Non-degenerate dual certificates). For any vector x ∈ R N , we define the set of non-degenerate dual certificates of x

D(x) = Im(Φ * ) ∩ ri(∂ J(x)) .

Stability in ℓ 2 Norm

The following theorem, proved in [101], establishes a linear convergence rate valid for any regularizer J, without any particular assumption beside being a proper closed convex function. In particular, its does not assume partial smoothness of J. This generic result encompasses many previous works, as discussed in Section 3.3.

Theorem 1. Assume that

Ker(Φ) ∩ T x 0 ∩ = {0} and D(x 0 ) = / 0 (9)
and consider the choice λ = c w , for some c > 0. Then we have for all minimizers x ⋆ of (P y,λ )

x ⋆ -x 0 2 C w , ( 10
)
where C > 0 is a constant (see Remark 4 for details).

In plain words, this bound tells us that the distance of x 0 to the set of minimizers of (P y,λ ) is within a factor of the noise level, which justifies the terminology "linear convergence rate".

Remark 2 (The role of non-smoothness). The injectivity of Φ when restricted to T x 0 is intimately related to the fact that J is non-smooth at x 0 . The higher the degree of non-smoothness, the lower the dimension of the subspace T x 0 , and hence the more likely the restricted injectivity. If J is smooth around x 0 (e.g. quadratic regularizers), however, the restricted injectivity condition cannot be fulfilled, unless Φ is itself injective. The reason is that T x 0 is the whole R N at the smoothness points. For smooth regularizations, it can be shown that the convergence rate is slower than linear, we refer to [198] for more details.

Remark 3 (Uniqueness). One can show that condition [START_REF] Bakin | Adaptive regression and model selection in data mining problems[END_REF] implies that x 0 is the unique solution of (P y,0 ) for y = Φx 0 . This condition however does not imply in general that (P y,λ ) has a unique minimizer for λ > 0.

Remark 4 (Stability constant). The result [START_REF] Bandeira | Certifying the restricted isometry property is hard[END_REF] ensures that the mapping y → x ⋆ (that might be set-valued) is C-Lipschitz-continuous at y = Φx 0 . Condition D(x 0 ) = / 0 is equivalent to the existence of some η ∈ D(x 0 ). The value of C (in fact an upperbound) can be found in [101]. It depends on Φ, T x 0 , c and the chosen non-degenerate dual certificate η. In particular, the constant degrades critically as η gets closer to the relative boundary of D(x 0 ), which reflects the intuition of how far is η from being a non-degenerate certificate.

Remark 5 (Source condition). The condition D(x 0 ) = / 0 is often called "source condition" or "range condition" in the literature of inverse problems. We refer to the monograph [198] for a general overview of this condition and its implications. It is an abstract condition, which is not easy to check in practice, since exhibiting a valid non-degenerate certificate is not trivial. We give in Section 4.1 further insights about this in the context of compressed sensing. Section 4.1 describes a particular construction of a good candidate (the so-called linearized pre-certificate) for being such an η ∈ D(x 0 ), and it is shown to govern stability of the manifold M x 0 for partly smooth regularizers.

Remark 6 (Infinite dimension).

It is important to remind that, in its full general form, Theorem 1 only holds in finite dimension. The constant C indeed may depend on the ambient dimension N, in which case the constant can blow-up as the discretization grid of the underlying continuous problem is made finer (i.e. as N grows). We detail below some relevant literature where similar results are shown in infinite dimension.

Related Works

Convergence Rates

For quadratic regularizations of the form J = D * • 2 for some linear operator D * , the ℓ 2 -error decay can be proved to be O( w ), which is not linear, see [198,Chapter 3] for more details and extensions to infinite dimensional Hilbert spaces. For non-smooth priors, in [START_REF] Burger | Convergence rates of convex variational regularization[END_REF], the authors show the Bregman distance between x ⋆ and x 0 exhibits a linear convergence rate for both the Lagrangian (P y,λ ) and the constrained (P 1 y,ε ) problems under the source condition D(x 0 ) = 0. These results hold more generally over infinite dimensional Banach spaces. They have been subsequently generalized to ill-posed non-linear inverse problems by [192] and [134]. It is important to observe that in order to prove convergence rates in terms of ℓ 2error, as done in [START_REF] Bandeira | Certifying the restricted isometry property is hard[END_REF], it is necessary to strengthen the source condition to its nondegenerate version, i.e. D(x 0 ) = 0.

In [155], the authors consider the case where J is a ℓ p norm with 1 p 2 and establish convergence rates of Φx 0 -Φx ⋆ in O( w ) and of x ⋆x 0 in O( w ).

[117] prove Theorem 1 for J = • 1 . They show that the non-degeneracy condition is also necessary for linear convergence, and draw some connections with the restricted isometry property (RIP), see below. Under a condition that bears similarities with ( 9), linear convergence with respect to J, i.e. J(x ⋆x 0 ) = O( w ), is proved in [116] for positively homogeneous regularizers. This result is equivalent to Theorem 1 but only when J is coercive, which precludes many important regularizers, such as for instance analysis-type regularizers including total variation.

RIP-based Compressed Sensing

The recovery performance of compressed sensing (i.e. when Φ is drawn from suitable random ensembles) for J = • 1 has been widely analyzed under the so-called restricted isometry property (RIP) introduced in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies? Information Theory[END_REF]. For any integer k 0, the k-th order restricted isometry constant of a matrix Φ is defined as the smallest

δ k 0 such that (1 -δ k ) x 2 Φx 2 (1 + δ k ) x 2 ,
for all vectors x such that x 0 k. It is shown [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] that if δ 2k + δ 3k < 1, then for every vector x 0 with x 0 0 k, there exists a non-degenerate certificate [START_REF] Candès | Decoding by linear programming[END_REF]Lemma 2.2], see also the discussion in [117]. In turn, this implies linear convergence rate, and is applied in [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF] to show ℓ 2 -stability to noise of compressed sensing. This was generalized in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] to analysis sparsity J = D * • 1 , where D is assumed to be a tight frame, structured sparsity in [START_REF] Candès | Compressed sensing with coherent and redundant dictionaries[END_REF] and matrix completion in [190,[START_REF] Candès | Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements[END_REF] using J = • * . The goal is then to design RIP matrices Φ with constants such that δ 2k + δ 3k (or a related quantity) is small enough. This is possible if Φ is drawn from an appropriate random ensemble for some (hopefully optimal) scaling of (N, P, k). For instance, if Φ is drawn from the standard Gaussian ensemble (i.e. with i.i.d. zero-mean standard Gaussian entries), there exists a constant C such that the RIP constants of Φ/ √ P obey δ 2k + δ 3k < 1 with overwhelming probability provided that P Ck log(N/k) , [START_REF] Barron | Risk bounds for model selection via penalization[END_REF] see for instance [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies? Information Theory[END_REF]. This result remains true when the entries of Φ are drawn independently from a subgaussian distribution. When Φ is a structured random matrix, e.g. random partial Fourier matrix, the RIP constants of Φ/ √ P can also satisfy the desired bound, but at the expense of polylog terms in the scaling [START_REF] Barron | Risk bounds for model selection via penalization[END_REF], see [105] for a comprehensive treatment. Note that in general, computing the RIP constants for a given matrix is an NP-hard problem [START_REF] Bandeira | Certifying the restricted isometry property is hard[END_REF]220].

RIP-less Compressed Sensing

RIP-based guarantees are uniform, in the sense that the recovery holds with high probability for all sparse signals. There is a recent wave of work in RIP-less analysis of the recovery guarantees for compressed sensing. The claims are non-uniform, meaning that they hold for a fixed signal with high probability on the random matrix Φ. This line of approaches improves on RIP-based bounds providing typically sharper constants. When Φ is drawn from the Gaussian ensemble, it is proved in [195] for J = • 1 that if the number of measurements P obeys P Ck log(N/k) for some constant C > 0, where k = x 0 0 , then condition ( 9) holds with high probability on Φ. This result is based on Gordon's comparison principle for Gaussian processes and depends on a summary parameter for convex cones called the Gaussian width. Equivalent lower bounds on the number of measurements for matrix completion from random measurements by minimizing the nuclear norm were provided in [START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF] to ensure that ( 9) holds with high probability. This was used to prove ℓ 2 -stable matrix completion in [START_REF] Candès | Matrix completion with noise[END_REF].

The authors in [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF] have recently showed that the Gaussian-width based approach leads to sharp lower bounds on P required to solve regularized inverse problems from Gaussian random measurements. For instance, they showed for J = • 1 that P > 2k log(N/k) [START_REF] Bauschke | A dykstra-like algorithm for two monotone operators[END_REF] guarantees exact recovery from noiseless measurements by solving (P y,0 ). An overhead in the number of measurements is necessary to get linear convergence of the ℓ 2 -error in presence of noise by solving (P 1 y,ε ) with ε = w , i.e. x 0 is feasible. Their results handle for instance the case of group sparsity [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF] and the nuclear norm [START_REF] Bach | Consistency of trace norm minimization[END_REF]. In the polyhedral case, it can be shown that [START_REF] Bauschke | A dykstra-like algorithm for two monotone operators[END_REF] implies the existence of a nondegenerate dual certificate, i.e. [START_REF] Bakin | Adaptive regression and model selection in data mining problems[END_REF], with overwhelming probability. The Gaussian width is closely related to another geometric quantity called the statistical dimension in conic integral geometry. The statistical dimension canonically extends the linear dimension to convex cones, and has been proposed in [START_REF] Amelunxen | Living on the edge: A geometric theory of phase transitions in convex optimization[END_REF] to deliver reliable predictions about the quantitative aspects of the phase transition for exact noiseless recovery from Gaussian measurements.

To deal with non-Gaussian matrix measurements (such as for instance partial Fourier matrices), [123] introduced the "golfing scheme" for noiseless low-rank matrix recovery guarantees using J = • * . The golfing scheme is an iterative procedure to construct an (approximate) non-degenerate certificate. This construction is also studied in [START_REF] Candès | A probabilistic and RIPless theory of compressed sensing. Information Theory[END_REF] for noiseless and noisy sparse recovery with J = • 1 . In another chapter of this volume [221], the author develops a technique, called the "bowling scheme", which is able to deliver bounds on the number of measurements that are similar with the Gaussian width-based bounds for standard Gaussian measurements, but the argument applies to a much wider class of measurement ensembles.

Model Stability

In the remainder of this chapter, we assume that J is finite-valued convex and locally partly smooth around x 0 , as defined in Section 2.2. This means in particular that the prior J promotes locally solution which belong to the manifold M = M x 0 . In the previous section, we were only concerned with ℓ 2 stability guarantees and partial smoothness was not necessary then. Owing to the additional structure conveyed by partial smoothness, we will be able to provide guarantees on the identification of the correct M = M x 0 by solving (P y,λ ), i.e. whether the (unique) solution x ⋆ of (P y,λ ) satisfies x ⋆ ∈ M . Such guarantees are of paramount importance for many applications. For instance, consider the case where ℓ 1 regularization is used to localize some (sparse) sources. Then x ⋆ ∈ M means that one perfectly identifies the correct source locations. Another example is that of the nuclear norm for low-rank matrix recovery. The correct model identification implies that x ⋆ has the correct rank, and consequently that the eigenspaces of x ⋆ have the correct dimensions and are close to those of x 0 .

Linearized Pre-certificate

We saw in Section 3.2 that ℓ 2 -stability of the solutions to (P y,λ ) is governed by the existence of a non-degenerate dual certificate p ∈ D(x 0 ). It turns out that not all dual certificates are equally good for stable model identification, and toward the latter, one actually needs to focus on a particular dual certificate, that we call "minimal norm" certificate.

Definition 5 (Minimal norm certificate). Assume that x 0 is a solution of (P y,0 ). We define the "minimal-norm certificate" as

η 0 = Φ * argmin Φ * p∈∂ J(x 0 ) p . (13)
A remarkable property, stated in Proposition 1 below, is that, as long as one is concerned with checking whether η 0 is non-degenerate, i.e. η 0 ∈ ri(∂ J(x 0 )), one can instead use the vector η F defined below, which can be computed in closed form.

Definition 6 (Linearized pre-certificate). Assume that

Ker(Φ) ∩ T x 0 = {0}. ( 14
)
We define the "linearized pre-certificate" as

η F = Φ * argmin Φ * p∈aff(∂ J(x 0 )) p . ( 15
)
Remark 7 (Well-posedness of the definitions). Note that the hypothesis that x 0 is a solution of (P y,0 ) is equivalent to saying that D(x 0 ) is a non-empty convex compact set. Hence in [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF], the optimal p is the orthogonal projection of 0 on a non-empty closed convex set, and thus η 0 is uniquely defined. Similarly, the hypothesis (14) implies that the constraint set involved in ( 15) is a non-empty affine space, and thus η F is also uniquely defined.

Remark 8 (Certificate vs. pre-certificate). Note that the only difference between ( 13) and ( 15) is that the convex constraint set ∂ J(x 0 ) is replaced by a simpler affine constraint. This means that η F does not always qualify as a valid certificate, i.e. η F ∈ ∂ J(x 0 ), hence the terminology "pre-certificate" used. This condition is actually at the heart of the model identification result exposed in Theorem 2.

For now on, let us remark that η F is actually simple to compute, since it amounts to solving a linear system in the least-squares sense.

Proposition 1. Under condition [START_REF] Bauschke | Dykstras algorithm with bregman projections: A convergence proof[END_REF], one has

η F = Φ * Φ +, * T x 0 e x 0 where e x 0 = P T x 0 (∂ J(x 0 )) ∈ R N . ( 16
)
Remark 9 (Computating e x ). The vector e x appearing in ( 16) can be computed in closed form for most of the regularizers discussed in Section 2.2. For instance, for J = • 1 , e x = sign(x). For J = • 1,B , it reads e x = (e b ) b∈B , where e b = x b / x b if x b = 0, and e b = 0 otherwise. For J = • * and a SVD decomposition

x = U x diag(Λ x )V * x , one has e x = U x V * x .
The following proposition, whose proof can be found in [231], exhibits a precise relationship between η 0 and η F . In particular, it implies that η F can be used in place of η 0 to check whether η 0 is non-degenerate, i.e. η 0 ∈ ri(∂ J(x 0 )).

Proposition 2. Under condition [START_REF] Bauschke | Dykstras algorithm with bregman projections: A convergence proof[END_REF], one has

η F ∈ ri(∂ J(x 0 )) =⇒ η F = η 0 , η 0 ∈ ri(∂ J(x 0 )) =⇒ η F = η 0 .

Model Identification

The following theorem provides a sharp sufficient condition to establish model selection. It is proved in [231]. It encompasses as special cases many previous works in the signal processing, statistics ans machine learning literatures, as we discuss in Section 4.5.1.

Theorem 2. Let J be locally partly smooth at x 0 relative to M = M x 0 . Assume that

Ker(Φ) ∩ T x 0 = {0} and η F ∈ ri(∂ J(x 0 )). ( 17
)
Then there exists C such that if

max(λ , w /λ ) C, ( 18
)
the solution x ⋆ of (P y,λ ) from the measurements (1) is unique and satisfies

x ⋆ ∈ M and x 0 -x ⋆ = O(max(λ , w )). ( 19
)
Remark 10 (Linear convergence rate vs. model identification). Obviously, the assumptions ( 17) of Theorem 2 imply those of Theorem 1. They are of course stronger, but imply a stronger result, since uniqueness of x ⋆ and model identification (i.e.

x ⋆ ∈ M ) are not guaranteed by Theorem 1 (which does not even need J to be partly smooth). A chief advantage of Theorem 2 is that its hypotheses can be easily checked and analyzed for a particular operator Φ. Indeed, computing η F only requires solving a linear system, as clearly seen from formula [START_REF] Beck | Gradient-based algorithms with applications to signal recovery[END_REF].

Remark 11 (Minimal signal-to-noise ratio). Another important distinction between

Theorems 1 and 2 is the second assumption [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF]. In plain words, it requires that the noise level is small enough and that the regularization parameter is wisely chosen. Such an assumption is not needed in Theorem 2 to ensure linear convergence of the ℓ 2 -error. In fact, this condition is quite natural. To see this, consider for instance the case of sparse recovery where J = • 1 . If the minimal signal-to-noise ratio is low, the noise will clearly dominate the amplitude of the smallest entries, so that one cannot hope to recover the exact support, but it is still possible to achieve a low ℓ 2 -error by forcing those small entries to zero.

Remark 12 (Identification of the manifold). For all the regularizations considered in Section 2.3, the conclusion of Theorem 2 is even stronger as it guarantees that M x ⋆ = M . The reason is that for any x and nearby points x ′ with x ′ ∈ M x , one has

M x ′ = M x .
Remark 13 (General loss/data fidelity). It is possible to extend Theorem 2 to account for general loss/data fidelity terms beyond the quadratic one, i.e. 1 2 y -Φx 2 . More precisely, this result holds true for loss functions of the form F(Φx, y), where F : R P × R P → R is a C 2 strictly convex function in its first argument, ∇F is C 1 in the second argument, with ∇F(y, y) = 0, where ∇F is the gradient with respect to the first variable. In this case, the expression ( 16) of η F becomes simply η F = Γ (P T Γ P T ) + e x 0 where

T = T x 0 Γ = Φ * ∂ 2 F(Φx 0 , Φx 0 )Φ ,
and where ∂ 2 F is the Hessian with respect to the first variable (which is a positive definite operator). We refer to [231] for more details.

Sharpness of the Model Identification Criterion

The following proposition, proved in [231], shows that Theorem 2 is in some sense sharp, since the hypothesis η F ∈ ri(∂ J(x 0 )) (almost) characterizes the stability of M .

Proposition 3. We suppose that x 0 is the unique solution of (P y,0 ) for y = Φx 0 and that Ker(Φ) ∩ T x 0 = {0}, and η

F / ∈ ∂ J(x 0 ). ( 20
)
Then there exists C > 0 such that if [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] holds, then any solution x ⋆ of (P y,λ ) for λ > 0 obeys x ⋆ / ∈ M .

In the particular case where w = 0 (no noise), this result shows that the manifold M is not correctly identified when solving (P y,λ ) for y = Φx 0 and for any λ > 0 small enough.

Remark 14 (Critical case). The only case not covered by neither Theorem 2 nor Proposition 3 is when η F ∈ rbd(∂ J(x 0 )), where rbd stands for the boundary relative to the affine hull. In this case, one cannot conclude, since depending on the noise w, one can have either stability or non-stability of M . We refer to [230] where an example illustrates this situation for the 1-D total variation J = D * DIF • 1 , where D * DIF is a finite-difference discretization of the 1-D derivative operator.

Probabilistic Model Consistency

Theorem 2 assumes a deterministic noise w, and the operator Φ is fixed. For applications in statistics and machine learning, it makes sense to rather assume a random model for both Φ and w. The natural question is then to assert that the estimator defined by solving (P y,λ ) is consistent in the sense that it correctly estimates x 0 and possibly the model M x 0 as the number of observations P → +∞. This requires to handle operators Φ with an increasing number of rows, and thus to also assess sensitivity of the optimization problem (P y,λ ) to perturbations of Φ (and not only to (w, λ ) as done previously).

To be more concrete, in this section, we work under the classical setting where N an x 0 are fixed as the number of observations P → +∞. The data (ϕ i , w i ) are assumed to be random vectors in R N × R, where ϕ i is the i-th row of Φ for i = 1, . . . , P. These vectors are supposed independent and identically distributed (i.i.d.) samples from a joint probability distribution such that E (w i |ϕ i ) = 0, finite fourthorder moments, i.e. E w 4 i < +∞ and E ϕ i 4 < +∞. Note that in general, w i and ϕ i are not necessarily independent. It is possible to consider other distribution models by weakening some of the assumptions and strengthening others, see e.g. [143,246,[START_REF] Bach | Consistency of the group Lasso and multiple kernel learning[END_REF]. Let us denote Γ = E(ϕ * i ϕ i ) ∈ R N×N , where ϕ i is any row of Φ. We do not make any assumption on the invertibility of Γ .

In this setting, a natural extension of η F defined by ( 16) in the deterministic case is ηF

= Γ Γ + T x 0 e x 0
where Γ T x 0 = P T x 0 Γ P T x 0 , and we use the fact that Γ T x 0 is symmetric and Im(Γ

+ T x 0 ) ⊂ T x 0 .
It is also implicitly assumed that Ker(Γ ) ∩ T x 0 = {0} which is the equivalent adaptation of the restricted injectivity condition in [START_REF] Birgé | From model selection to adaptive estimation[END_REF] to this setting.

To make the discussion clearer, the parameters (λ = λ P , Φ = Φ P , w = w P ) are now indexed by P. The estimator x ⋆ P obtained by solving (P λ P ,y P ) for y P = Φ P x 0 + w P is said to be consistent for x 0 if, lim P→+∞ Pr (x ⋆ P is unique) = 1 and x ⋆ P → x 0 in probability. The estimator is said to be model consistent if

lim P→+∞ Pr (x ⋆ P ∈ M ) = 1,
where M = M x 0 is the manifold associated to x 0 .

The following result, whose proof can be found in [231], guarantees model consistency for an appropriate scaling of µ P . It generalizes several previous works in the statistical and machine learning literature as we review in Section 4.5.1.

Theorem 3. If Ker(Γ ) ∩ T x 0 = {0} and ηF ∈ ri(∂ J(x 0 )), (21)
and λ P = o(P) and λ -1

P = o(P -1/2 ). ( 22
)
Then the estimator x ⋆ P of x 0 is model consistent.

Related Works

Model Consistency

Theorem 2 is a generalization of a large body of results in the literature. For the Lasso, i.e. J = • 1 , to the best of our knowledge, this result was initially stated in [107]. In this setting, the result ( 19) corresponds to the correct identification of the support, i.e. supp(x ⋆ ) = supp(x 0 ). Condition ( 21) for J = • 1 is known in the statistics literature under the name "irrepresentable condition" (generally stated in a non-geometrical form), see e.g. [246].

[143] have shown estimation consistency for Lasso for fixed N and x 0 and asymptotic normality of the estimates. The authors in [246] prove Theorem 3 for J = • 1 , though under slightly different assumptions on the covariance and noise distribution. A similar result is established in [141] for the elastic net, i.e. J = • 1 + ρ • 2 2 for ρ > 0. In [START_REF] Bach | Consistency of the group Lasso and multiple kernel learning[END_REF] and [START_REF] Bach | Consistency of trace norm minimization[END_REF], the author proves Theorem 3 for two special cases, namely the group Lasso and nuclear norm minimization. Note that these previous works assume that the asymptotic covariance Γ is invertible. We do not impose such an assumption, and only require the weaker restricted injectivity condition Ker(Γ ) ∩ T = {0}. In a previous work [230], we have proved an instance of Theorem 2 when J(x) = D * x 1 , where D ∈ R N×Q is an arbitrary linear operator. This covers as special cases the discrete anisotropic total variation or the fused Lasso. This result was further generalized in [229] when J belongs to the class of partly smooth functions relative to linear manifolds M , i.e. M = T x . Typical instances encompassed in this class are the ℓ 1 -ℓ 2 norm, or its analysis version, as well as polyhedral gauges including the ℓ ∞ norm. Note that the nuclear norm (and composition of it with linear operators as proposed for instance in [118,193]), whose manifold is not linear, does not fit into the framework of [229], while it is covered by Theorem 2. Lastly, a similar result is proved in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] for a continuous (infinite dimensional) sparse recovery problem over the space of Radon measures normed by J the total variation of a measure (not to be confused with the total variation of functions). In this continuous setting, an interesting finding is that, when η 0 ∈ ri(∂ J(x 0 )), η 0 is not equal to η F but to a different certificate (called "vanishing derivative" certificate in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]) that can also be computed by solving a linear system.

Stronger Criteria for ℓ 1

Many sufficient conditions have been proposed in the literature to ensure that η F is a non-degenerate certificate, and hence to guarantee stable identification of the support (i.e. model). We illustrate this here for J = • 1 , but similar reasoning can be carried out for • 1,B or • * .

The strongest criterion makes use of mutual coherence, first considered in [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF] µ

(Φ) = max i = j | ϕ i , ϕ j |
where each column ϕ i of Φ are assumed normalized to a unit ℓ 2 norm. Mutual coherence measures the degree of ill-conditioning of Φ through the correlation of its columns (ϕ i ) 1 i N . Mutual coherence is always lower-bounded by N-P P(N-1) , and equality holds if and only if (ϕ i ) 1 i N is an equiangular tight frame, see [210]. Finer variants based on cumulative coherences have been proposed in [120,[START_REF] Borup | Beyond coherence : recovering structured timefrequency representations[END_REF]. To take into account the influence of the support I = supp(x 0 ) of the vector x 0 to recover, Tropp introduced in [222] the Exact Recovery Condition (ERC), defined as

ERC(I) = Φ * I c Φ +, * I ∞,∞ = max j / ∈I Φ + I ϕ j 1
where • ∞,∞ is the matrix operator norm induced by the ℓ ∞ vector norm, Φ I = (ϕ i ) i∈I and I c is the complement of the set I. Φ I is assumed injective which, in view of Section 2.3.1, is nothing but a specialization to ℓ 1 of the restricted injectivity condition in [START_REF] Birgé | From model selection to adaptive estimation[END_REF]. A weak ERC criterion, which does not involve matrix inversion, is derived in [START_REF] Dossal | Sparse spike deconvolution with minimum scale[END_REF] wERC(I) =

max j∈I c ∑ i∈I | ϕ i , ϕ j | 1 -max j∈I ∑ i = j∈I | ϕ i , ϕ j | .
Given the structure of the subdifferential of the ℓ 1 norm, it is easy to check that

η F ∈ ri(∂ J(x 0 )) ⇐⇒ IC(x 0 ) = Φ * I c Φ +, * I sign(x 0,I ) ∞ < 1.
The right hand side in the equivalence is precisely what is called the irrepresentable condition in statistics and machine learning. Clearly, IC(x 0 ) involves both the sign vector and the support of x 0 . The following proposition gives ordered upper bounds of IC(x 0 ) in terms of the cruder criteria ERC, wERC and mutual coherence. A more elaborate discussion of them can be found in [160].

Proposition 4. Assume that Φ I is injective and denote k = |I| = x 0 0 . Then,

IC(x 0 ) ERC(I) wERC(I) kµ(Φ) 1 -(k -1)µ(Φ)
.

Linearized Pre-certificate for Compressed Sensing Recovery

Stable support identification has been established in [240,[START_REF] Dossal | Sharp support recovery from noisy random measurements by ℓ 1 -minimization[END_REF] for the Lasso problem when Φ is drawn from the Gaussian ensemble. These works show that for k = x 0 0 , if P > 2k log(N) then indeed η F ∈ ri(∂ J(x 0 )), and this scaling can be shown to be sharp. This scaling should be compared with [START_REF] Bauschke | A dykstra-like algorithm for two monotone operators[END_REF] ensuring that there exists a non-degenerate certificate.

The gap in the log term indicates that there exists vectors that can be stably recovered by ℓ 1 minimization in ℓ 2 -error sense, but whose support cannot be stably identified. Equivalently, for these vectors, there exists a non-degenerate certificate but it is not η F . The pre-certificate η F is also used to ensure exact recovery of a low-rank matrix from incomplete noiseless measurements by minimizing the nuclear norm [START_REF] Candès | Exact matrix completion via convex optimization[END_REF][START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF]. This idea is further generalized by [START_REF] Candès | Simple bounds for recovering low-complexity models[END_REF] for a family of decomposable norms (including in particular ℓ 1 -ℓ 2 norm and the nuclear norm), which turns to be a subset of partly smooth regularizers. In these works, lower bounds on the number of random measurements needed for η F to be a non-degenerate certificate are developed. In fact, these measurement lower bounds combined with Theorem 2 allow to conclude that matrix completion by solving (P y,λ ) with J = • * identifies the correct rank at high signal-to-noise levels.

Sensitivity Analysis

Sensitivity analysis is a central theme in variational analysis. Comprehensive monographs on the subject are [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]166]. The function to be analyzed underlying problems (P y,λ ) and (P y,0 ) is

f (x, θ ) = 1 2λ y -Φx 2 + J(x) if λ > 0, ι H y (x) + J(x) if λ = 0, , (23)
where H y = {y ; Φx = y} and where the parameters are θ = (λ , y, Φ) for λ 0. Theorems 2 and 3 can be understood as a sensitivity analysis of the minimizers of f at a point (x = x 0 , θ = θ 0 = (0, Φx 0 , Φ)).

Classical sensitivity analysis of non-smooth optimization problems seeks conditions to ensure smoothness of the mapping θ → x θ where x θ is a minimizer of f (•, θ ), see for instance [194,[START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]. This is usually guaranteed by the non-degenerate source condition and restricted injectivity condition [START_REF] Bakin | Adaptive regression and model selection in data mining problems[END_REF], which, as already exposed in Section 3.2, ensure linear convergence rate, and hence Lipschitz behavior of this mapping. The analysis proposed by Theorem 2 goes one step further, by assessing that M x 0 is a stable manifold (in the sense of [241]), since the minimizer x θ is unique and remains in M x 0 for θ close to θ 0 . Our starting point for establishing Theorem 2 is the inspiring work of Lewis [146] who first introduced the notion of partial smoothness and showed that this broad class of functions enjoys a powerful calculus and sensitivity theory. For convex functions (which is the setting considered in our work), partial smoothness is closely related to U -V -decompositions developed in [145]. In fact, the behavior of a partly smooth function and of its minimizers (or critical points) depend essentially on its restriction to the manifold, hence offering a powerful framework for sensitivity analysis theory. In particular, critical points of partly smooth functions move stably on the manifold as the function undergoes small perturbations [149]. A important and distinctive feature of Theorem 2 is that, partial smoothness of J at x 0 relative to M transfers to f (•, θ ) for λ > 0, but not when λ = 0 in general. In particular, [146,Theorem 5.7] does not apply to prove our claim.

Sensitivity Analysis and Parameter Selection

In this section, we study local variations of the solutions of (P y,λ ) considered as functions of the observations y. In a variational-analytic language, this corresponds to analyzing the sensitivity of the optimal values of (P y,λ ) to small perturbations of y seen as a parameter. This analysis will have important implications, and we exemplify one of them by constructing unbiased estimators of the quadratic risk, which in turn will allow to have an objectively-guided way to select the optimal value of the regularization parameter λ .

As argued in Section 4.5.4, assessing the recovery performance by solving (P y,λ ) for w and λ small amounts to a sensitivity analysis of the minimizers of f in ( 23) at (x = x 0 , θ = θ 0 = (0, Φx 0 , Φ)). This section involves again sensitivity analysis of ( 23) to perturbations of y but for λ > 0. Though we focus our attention on sensitivity to y, our arguments extend to any parameters, for instance λ or Φ.

Similarly to the previous section, we suppose here that J is a finite-valued convex and partly smooth function. For technical reasons, we furthermore assume that the partial smoothness manifold is linear, i.e. M x = T x . We additionally suppose that the set of all possible models T = {T x } x∈R N is finite. All these assumptions hold true for the regularizers considered in Section 2.3, with the notable exception of the nuclear norm, whose manifolds of partial smoothness are non-linear.

Differentiability of Minimizers

Let us denote x ⋆ (y) a minimizer of (P y,λ ) for a fixed value of λ > 0. Our main goal is to study differentiability of x ⋆ (y) and find a closed-form formula of the derivative of x ⋆ (y) with respect to the observations y. Since x ⋆ (y) is not necessarily a unique minimizer, such a result means actually that we have to single out one solution x ⋆ (y), which hopefully should be a locally smooth function of y. However, as J is non-smooth, one cannot hope for such a result to hold for any observation y ∈ R P . For applications to risk estimation (see Section 5.3), it is important to characterize precisely the smallest set H outside of which x ⋆ (y) is indeed locally smooth. It turns out that one can actually write down an analytical expression of such a set H , containing points where one cannot find locally a smooth parameterization of the minimizers. This motivates our definition of what we coin a "transition space". Definition 7 (Transition space). We define the transition space H as

H = T ∈T bd(H T ),
where bd(C) is the boundary of a set C, and

H T = y ∈ R P ; ∃x ∈ T , λ -1 Φ * T (Φx -y) ∈ rbd(∂ J(x)) .
where

T = x ∈ R N ; T x = T .
The set H contains the observations y ∈ R P such that the model subspace T x(y) associated to a well chosen solution x(y) of (P y,λ ) is not stable with respect to small perturbations of y. In particular, when J = • 1 , it can be checked that H is a finite union of hyperplanes and when J = • 1,2 it is a semi-algebraic set (see Definition 8). This stability is not only crucial to prove smoothness of x(y), it is also important to be able to write down an explicit formula for the derivative, as detailed in the following theorem whose proof is given in [227].

Theorem 4. Let y ∈ H and x ⋆ a solution of (P y,λ ) such that

Ker Φ T ∩ Ker D 2 J T (x ⋆ ) = {0} (I x ⋆ )
where T = T x ⋆ . Then, there exists an open neighborhood V ⊂ R N of y, and a mapping x : V → T such that 1. for every ȳ ∈ V , x( ȳ) is a solution of (P λ , ȳ), and x(y) = x ⋆ ; 2. the mapping x is C 1 (V ) and

∀ ȳ ∈ V , D x( ȳ) = (Φ * T Φ T + λ D 2 J T (x ⋆ )) -1 Φ T .
Here D 2 J T is the Hessian (second order derivative) of J restricted to T . This Hessian is surely well-defined owing to partial smoothness, see Definition 2(i).

Semi-algebraic Geometry

Our goal now is to show that the set H is in some sense "small" (in particular to show that it has zero Lebesgue measure), which will entail differentiability of y → x ⋆ Lebesgue almost everywhere. For this, additional geometrical structure on J is needed. Such a rich class of functions is provided by the notion of a semialgebraic subset of R N to be defined shortly. Semi-algebraic sets and functions have been broadly applied to various areas of optimization. The wide applicability of semi-algebraic functions follows largely from their stability under many mathematical operations. In particular, the celebrated Tarski-Seidenberg theorem states, loosely, that the projection of a semi-algebraic set is semi-algebraic. These stability properties are crucial to obtain the following result, proved in [227].

Definition 8 (Semi-algebraic set and function).

A set E is semi-algebraic if it is a finite union of sets defined by polynomial equations and (possibly strict) inequalities. A function f :

E → F is semi-algebraic if E and its graph {(u, f (u)) ; u ∈ E} are semi-algebraic sets.
Remark 15 (From semi-algebraic to o-minimal geometry). The class of semi-algebraic functions is large, and subsumes, for instance, all the regularizers J described in Section 2.3. The qualitative properties of semi-algebraic functions are shared by a much bigger class called functions definable in an o-minimal structure over R, or simply definable functions. O-minimal structures over R correspond in some sense to an axiomatization of some of the prominent geometrical properties of semialgebraic geometry [START_REF] Coste | An introduction to o-minimal geometry[END_REF] and particularly of the stability under projection. For example, the function J(x) = ∑ i |x i | s , for an arbitrary s 0, is semi-algebraic only for rational s ∈ Q, while it is always definable in an o-minimal structure [236]. Due to the variety of regularizations J that can be formulated within the framework of ominimal structures, all our results stated in this section apply to definable functions, see [227] for a detailed treatment.

Semi-algebraic functions are stable for instance under (sub)differentiation and projection. These stability properties are crucial to obtain the following result, proved in [227].

Proposition 5. If J is semi-algebraic, the transition space H is semi-algebraic and has zero Lebesgue measure.

Unbiased Risk Estimation

A problem of fundamental practical importance is to automatically adjust the parameter λ to reach the best recovery performance when solving (P y,λ ). Parameter selection is a central theme in statistics, and is intimately related to the question of model selection, as introduced in Section 2.1.

We then adopt a statistical framework in which the observation model ( 1) becomes

Y = Φx 0 + W ( 24
)
where W is random noise having an everywhere strictly positive probability density function, assumed to be known. Though the forthcoming results can be stated for a large family of distributions, for the sake of concreteness, we only consider the white Gaussian model where W ∼ N (0, σ 2 Id P×P ), with known variance σ 2 .

Under the observation model [START_REF] Borup | Beyond coherence : recovering structured timefrequency representations[END_REF], the ideal choice of λ should be the one which minimizes the quadratic estimation risk E W ( x ⋆ (Y )x 0

2 ). This is obviously not realistic as x 0 is not available, and in practice, only one realization of Y is observed. To overcome these obstacles, the traditional approach is to replace the quadratic risk with an some estimator that solely depends on Y . The risk estimator is also expected to enjoy nice statistical properties among which unbiasedness is highly desirable.

However, it can be shown, see e.g. [100, Section IV], that the quadratic risk

E W ( x ⋆ (Y ) -x 0 2
) cannot be reliably estimated on Ker(Φ). Nonetheless, we may still obtain a reliable assessment of the part that lies in Im(Φ * ) = Ker(Φ) ⊥ or any linear image of it. For instance, the most straightforward surrogate of the above risk is the so-called prediction risk E W ( µ(Y ) -µ 0 2 ), where

µ 0 = Φx 0 and µ(y) = Φx ⋆ (y),
where x ⋆ (y) is any solution of (P y,λ ). One can easily show that µ(y) ∈ R P is well-defined as a single-valued mapping and thus does not depend on the particular choice of x ⋆ (y), see [227]. Consequently, Theorem 4 shows that y → µ(y) is a C 1 mapping on R P \ H .

Degrees of Freedom

The degrees of freedom (DOF) quantifies the model "complexity" of a statistical modeling procedure [START_REF] Efron | How biased is the apparent error rate of a prediction rule[END_REF]. It is at the heart of several risk estimation procedures. Therefore, in order to design estimators of the prediction risk, an important step is to get an estimator of the corresponding DOF.

Definition 9 (Empirical DOF). Suppose that y → µ(y) is differentiable Lebesgue almost everywhere, as is the case when it is Lipschitz-continuous (Rademacher's theorem). The empirical number of degrees of freedom is defined as

df(y) = div(µ)(y) = tr(Dµ(y)),
where the derivative is to be understood in the weak sense, i.e. to hold Lebesgue almost everywhere (a.e.).

An instructive example to get the gist of this formula is the case where µ is the orthogonal projection onto some linear subspace V . We then get easily that df(y) = dim(V ), which is in agreement with the intuitive notion of the number of DOF.

The following result delivers the closed-form expression of df(y), valid on a full Lebesgue measure set, for µ(y) = Φx ⋆ (y) and x ⋆ (y) an appropriate solution of (P y,λ ). At this stage, it is important to realize that the main difficulty does not lie in showing almost everywhere differentiability of µ(y); this mapping is in fact Lipschitz-continuous by classical arguments of sensitivity analysis applied to (P y,λ ). Rather, it is the existence of such a formula and its validity Lebesgue a. e. that requires more subtle arguments obtained owing to partial smoothness of J. For this, we need also to rule out the points y where (I x ⋆ ) does not hold. This is the rationale behind the following set.

Definition 10 (Non-injectivity set). We define the Non-injectivity set G as G = y / ∈ H ; (I x ⋆ ) does not hold for any minimizer x ⋆ of (P y,λ ) .

Theorem 5. For every y / ∈ H ∪ G , there is x ⋆ such that (P y,λ ) holds and

df(y) = tr(∆ x ⋆ (y)) where ∆ x ⋆ (y) = Φ T • (Φ T * Φ T + λ D 2 J T (x ⋆ )) -1 • Φ T * , ( 25
)
were T = T x ⋆ .

Remark 16 (Non-injectivity set). It turns out that G is in fact empty for many regularizers. This is typically the case for J = • 1 [START_REF] Dossal | The degrees of freedom of the Lasso for general design matrix[END_REF], J = D * • 1 [228], and the underlying reasoning can be more generally extended to polyhedral regularizers. The same result was also shown for J = • 1,2 in [232]. More precisely, in all these works, it was shown that for each y / ∈ H , there exists a solution x ⋆ of (P y,λ ) that fulfills (I x ⋆ ). The proof is moreover constructive allowing to build such a solution starting from any other one.

Stein Unbiased Risk Estimator (SURE)

We now have all necessary ingredients at hand to design an estimator of the prediction risk.

Definition 11. Suppose that y → µ(y) is differentiable Lebesgue almost everywhere, as is the case when it is Lipschitz-continuous. The SURE is defined as

SURE(y) = y -µ(y) 2 + 2σ 2 df(y) -Pσ 2 . ( 26
)
In this definition, we have anticipated on unbiasedness of this estimator. In fact, this turns out to be a fundamental property owing to the celebrated lemma of Stein [209], which indeed asserts that the SURE ( 26) is an unbiased estimator of the prediction risk. Therefore, putting together Theorem 5, Proposition 5 and Stein's lemma, we get the following. Theorem 6. Suppose that J is semi-algebraic and G is of zero Lebesgue measure. Then,

E W (SURE(Y )) = E W ( µ(Y ) -µ 0 2 )
where ( 25) is plugged into [START_REF] Bredies | Total generalized variation[END_REF], and µ(Y ) = Φx ⋆ (Y ).

Remark 17 (Parameter selection).

A practical usefuleness of the SURE is its ability to provide an objectively guided way to select a good λ from a single observation y by minimizing SURE(y). While unbiasedeness of the SURE is guaranteed, it is hard to control its variance and hence its consistency. This is an open problem in general, and thus little can be said about the actual theoretical efficiency of such an empirical parameter selection method. It works however remarkably well in practice, see the discussion in Section 5.6.5 and references therein.

Remark 18 (Projection risk). The SURE can be extended to unbiasedly estimate other risks that the prediction one. For instance, as argued in Section 5.3, one can estimate the so-called projection risk defined as E W ( P Ker(Φ) ⊥ (x ⋆ (Y )x 0 ) 2 ). This is obviously better that the prediction risk as a surrogate for the estimation risk.

Related Works

Sensitivity Analysis

In Section 4.5.4, we reviewed the relevant literature pertaining to sensitivity analysis for partly smooth functions, which is obviously very connected to Theorem 4. See also [START_REF] Bolte | Generic optimality conditions for semialgebraic convex programs[END_REF] for the case of linear optimization over a convex semi-algebraic partly smooth feasible set, where the authors prove a sensitivity result with a zero-measure transition space. A distinctive feature of our analysis toward proving unbiasedness of the SURE is the need to ensure that sensitivity analysis can be carried out on a full Lebesgue measure set. In particular, it necessitates local stability of the manifold M x ⋆ associated to an appropriate solution x ⋆ , and this has to hold Lebesgue almost everywhere. Thus the combination of partial smoothness and semi-algebraicity is the key.

Risk Estimators

In this section, we put emphasis on the SURE as an unbiased estimator of the prediction risk. There are other alternatives in the literature which similarly rely on estimator of the DOF. One can think for instance of the generalized cross-validation (GCV) [115]. Thus our results apply equally well to such risk estimators. Extensions of the SURE to independent variables from a continuous exponential family are considered in [135].

[100] generalizes the SURE principle to continuous multivariate exponential families, see also [181,228] for the multivariate Gaussian case.

The results described here can be extended to these setting as well, see [227].

Applications of SURE in Statistics and Imaging

Applications of SURE emerged for choosing the parameters of linear estimators such ridge regression or smoothing splines [150]. After its introduction in the wavelet community through the SURE-Shrink estimator [START_REF] Donoho | Adapting to Unknown Smoothness Via Wavelet Shrinkage[END_REF], it has been extensively used for various image restoration problems, e.g. with sparse regularization [START_REF] Blu | The SURE-LET approach to image denoising[END_REF]238,186,[START_REF] Chaux | A nonlinear stein-based estimator for multichannel image denoising[END_REF]181,[START_REF] Cai | A data-driven block thresholding approach to wavelet estimation[END_REF]156,187,188] or with non-local means [234,[START_REF] Duval | A bias-variance approach for the non-local means[END_REF][START_REF] Deledalle | Non-local Methods with Shape-Adaptive Patches (NLM-SAP)[END_REF]235].

Closed-form Expressions for SURE

For the Lasso problem, i.e. J = • 1 , the divergence formula [START_REF] Boyd | Convex optimization[END_REF] reads

df(y) = | supp(x ⋆ )|,
where x ⋆ is a solution of (P y,λ ) such that (I x ⋆ ) holds, i.e. Φ supp(x ⋆ ) has full rank. This result is proved in [247] for injective Φ and in [START_REF] Dossal | The degrees of freedom of the Lasso for general design matrix[END_REF] for arbitrary Φ. This result is extended to analysis ℓ 1 -sparsity, i.e. J = D * • 1 , in [216,228]. A formula for the DOF in the case where x ⋆ (y) is the orthogonal projection onto a partly smooth convex set C is proved in [142]. This work extends that of [164] which treats the case where C is a convex polyhedral cone. These two works allow one to compute the degrees of freedom of estimators defined by solving (P 2 y,γ ) in the case where Φ is injective. [128] studied the DOF of the metric projection onto a closed set (non-necessarily convex), and gave a precise representation of the bias when the projection is not sufficiently differentiable.

A formula of an estimate of the DOF for the group Lasso, i.e. J = • 1,2 when Φ is orthogonal within each group was conjectured in [244]. An estimate is also given by [202] using heuristic derivations that are valid only when Φ is injective, though its unbiasedness is not proved. [226] derived an estimator of the DOF of the group Lasso and proved its unbiasedness when Φ is injective. Closed-form expression of the DOF estimate for denoising with the nuclear norm, i.e. Φ = Id and J = • * , were concurrently provided in [START_REF] Deledalle | Risk estimation for matrix recovery with spectral regularization[END_REF][START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF].

Numerical Methods for SURE

Deriving the closed-form expression of the DOF is in general challenging and has to be addressed on a case by case basis. The implementation of the divergence formula such as ( 25) can be computationally expensive in high dimension. But since only the trace of the Jacobian is needed, it is possible to speed up these computations through Monte-Carlo sampling, but at the price of mild approximations. If the Jacobian is not known in closed-form or prohibitive to compute, one may appeal to finite-difference approximations along Monte Carlo sampled directions [243,201], see [111,186] for applications to imaging problems.

In practice, the analytical formula ( 25) might be subject to serious numerical instabilities, and thus cannot always be applied safely when the solution x ⋆ is only known approximately. Think for instance of the case where x ⋆ is approximated by an an iterate computed after finitely-many iterations of an algorithm as detailed in Section 6. A better practice is then to directly compute the DOF, hence the SURE, recursively from the iterates themeselves, as proposed by [238,112,[START_REF] Deledalle | Proximal Splitting Derivatives for Risk Estimation[END_REF].

Proximal Splitting for Structured Optimization

Though problems (P y,λ ), (P y,0 ), (P 1 y,ε ) or (P 2 y,γ ) are non-smooth, they enjoy enough structure to be solved by efficient algorithms. The type of algorithm to be used depends in particular on the properties of J. We first briefly mention some popular non-smooth optimization schemes in Section 6.1, and focus our attention on proximal splitting schemes afterwards.

Convex Optimization for Regularized Inverse Problems

(Sub)-gradient Descent

Consider for example problem (P y,λ ). This is a convex composite optimization problem where one of the functions is smooth with a Lipschitz-continuous gradient. If J were smooth enough, then a simple gradient (or possibly (quasi-)Newton) descent method could be used. However, as detailed in Section 2.2, low-complexity regularizers J are intended to be non-smooth in order to promote models M of low intrinsic dimension, and J is precisely non-smooth transverse to M . One can think of replacing gradients by subgradients (elements of the subdifferential), since J is assumed finite-valued (hence closed) convex, which are bounded. This results in a subgradient descent algorithm which is guaranteed to converge but under stringent assumptions on the descent step-sizes, which in turn makes their global convergence rate quite slow, see [172].

Interior Point Methods

Clearly, the key to getting efficient algorithms is to exploit the structure of the optimization problems at hand while handling non-smoothness properly. For a large class of regularizers J, such as those introduced in Section 2.3, the corresponding optimization problems can be cast as conic programs. The cone constraint can be enforced using a self-concordant barrier function, and the optimization problem can hence be solved using interior point methods, as pioneered by [174], see also the monograph [START_REF] Boyd | Convex optimization[END_REF]. This class of methods enjoys fast convergence rate. Each iteration however is typically quite costly and can become prohibitive as the dimension increases.

Conditional Gradient

This algorithm is historically one of the first method for smooth constrained convex optimization (a typical example being (P 2 y,γ )), and was extensively studied in the 70's. It is also known as Frank-Wolfe algorithm, since it was introduced by [106] for quadratic programming and extended in [START_REF] Dunn | Conditional gradient algorithms with open loop step size rules[END_REF]. The conditional gradient algorithm is premised on being able to easily solve (at each iteration) linear optimization problems over the feasible region of interest. This is in contrast to other first-order methods, such as forward-backward splitting and its variants (see Section 6.3), which are premised on being able to easily solve (at each iteration) a projection problem. Moreover, in many applications the solutions to the linear optimization subproblem are highly structured and exhibit particular sparsity and/or low-rank properties. These properties have renewed interest in the conditional gradient method to solve sparse recovery (ℓ 1 and total variation), low-rank matrix recovery (nuclear norm minimization), anti-sparsity recovery, and various other problems in signal processing and machine learning; see e.g. [START_REF] Clarkson | Coresets, sparse greedy approximation, and the frank-wolfe algorithm[END_REF]136,200,[START_REF] Dudík | Lifted coordinate descent for learning with tracenorm regularization[END_REF]129].

Homotopy/Path-following

Homotopy and path-following-type methods have been introduced in the case of ℓ 1minimization to solve (P y,λ ) by [176]. They were then adapted to analysis ℓ 1 , i.e. J = D * • 1 , in [215], and ℓ ∞ regularization, • ∞ , in [108]. One can in fact show that these methods can be applied to any polyhedral regularization (see [233]), because these methods only rely on the crucial fact that the solution path λ → x ⋆ λ , where x ⋆ λ is a solution of (P y,λ ), is piecewise affine. The LARS algorithm [START_REF] Efron | Least angle regression[END_REF] is an accelerated version of homotopy which computes an approximate homotopy path for J = • 1 along which the support increases monotonically along the course of iterations. In the noiseless compressed sensing case, with Φ drawn from the Gaussian ensemble, it is shown in [START_REF] Donoho | Fast solution of ℓ 1 -norm minimization problems when the solution may be sparse[END_REF] that if x 0 is k-sparse with P > 2k log(N), the homotopy method reaches x 0 in only k iterations. This k-solution property was empirically observed for other random matrix ensembles, but at different thresholds for P. In [158], the authors proved that in the worst case, the number of segments in the solution path is exponential in the number of variables, and thus the homotopy method can then take as many iterations to converge.

As for interior points, the cost per iteration of homotopy-like methods, without particular ad hoc optimization, scales badly with the dimension, thus preventing them to be used for large-scale problems such as those encountered in imaging. This class of solvers is thus a wise choice for problems of medium size, and when high accuracy (or even exact computation up to machine precision for the homotopy algorithm) is needed. Extensions of these homotopy methods can deal with progressive changes in the operator Φ or the observations y, and are thus efficient for these settings, see [START_REF] Asif | Sparse recovery of streaming signals using L1-homotopy[END_REF].

Approximate Message Passing

In the last five years, ideas from graphical models and message passing and approximate message passing algorithms have been proposed to solve large-scale problems of the form (P y,λ ) for various regularizers J, in particular ℓ 1 , ℓ 1 -ℓ 2 and the nuclear norm. A comprehensive review is given in [165]. However, rigorous convergence results have been proved so far only in the case in which Φ is standard Gaussian, though numerical results show that the same behavior should apply for broader random matrix ensembles.

Proximal Splitting Algorithms

Proximal splitting methods are first-order iterative algorithms that are tailored to solve structured non-smooth (essentially convex) optimization problems. The first operator splitting method has been developed from the 70's. Since then, the class of splitting methods have been regularly enriched with increasingly sophisticated algorithms as the structure of problems to handle becomes more complex.

To make our discussion more concrete, consider the general problem of minimizing the proper closed convex function

f = h + K ∑ k=1 g k • A k
where h : R N → R is convex and smooth, the A k : R N → R N k are linear operators and g k : R N k → R are proper closed convex functions for which the so-called proximity operator (to be defined shortly) can be computed easily (typically in closed form). We call such a function g k "simple". Definition 12. The proximity operator of a proper closed convex function g is defined as, for γ > 0

prox γg (x) = argmin u∈R N 1 2
xu 2 + γg(u).

The proximal operator generalizes the notion of orthogonal projection onto a nonempty closed convex set C that one recovers by taking g = ι C . Proximal splitting algorithms may evaluate (possibly approximately) the individual operators (e.g. gradient of h), the proximity operators of the g k 's, the linear operators A k , all separately at various points in the course of iteration, but never those of sums of functions nor composition by a linear operator. Therefore, each iteration is cheap to compute for large-scale problems. They also enjoy rigorous convergence guarantees, stability to errors, with possibly quantified convergence rates and iteration complexity bounds on various quantities. This justifies their popularity in contemporary signal and image processing or machine learning, despite that their convergence is either sublinear or at best linear.

It is beyond the scope of this Chapter to describe thoroughly the huge literature on proximal spliting schemes, as it is a large and extremely active research field in optimization theory. Good resources and reviews on the subject are [START_REF] Beck | Gradient-based algorithms with applications to signal recovery[END_REF][START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]178]. We instead give a brief classification of the most popular algorithms according to the class of structured objective functions they are able to handle:

• Forward-Backward (FB) algorithm [163,179,[START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. It is designed to minimize (6.2) when h has a Lipschitz-continuous gradient, K = 1, A 1 = Id, and g 1 is simple. There are accelerated (optimal) variants of FB, such as the popular Nesterov [173] or Fista [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], but the convergence of the iterates is not longer guaran-teed for these schemes. FB and its variants are good candidates to solve (P y,λ ). We will further elaborate on FB in Section 6.3. • Douglas-Rachford (DR) algorithm [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF]153]. It is designed to minimize (6.2) for h = 0, K = 2, A k = Id and g k is simple for k = 1, 2. It can be easily extended to the case of K > 2 by either lifting to a product space, see e.g. [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], or through projective splitting [START_REF] Eckstein | General projective splitting methods for sums of maximal monotone operators[END_REF]. DR can be used to solve (P y,0 ), (P 1 y,ε ) or (P 2 y,γ ) for certain operators Φ.

• Generalized Forward-Backward (GFB) algorithm [185]. It can handle the case of an arbitrary K with A k = Id, g k simple and h has a Lipschitz-continuous gradient.

It can be interpreted as hybridization of FB scheme and the DR scheme on a product space.

• Alternate Direction Method of Multipliers (ADMM) algorithm [104,109,110,113]. It is adapted to minimize (6.2) for h = 0, K = 2 with A 1 = Id and A 2 is injective. It can be shown [110,[START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF] that ADMM is equivalent to DR applied to the Fenchel-Rockafellar dual problem min u g * 1 • -A * 2 (u) + g * 2 (u), where g * k is the Legendre-Fenchel conjugate of g k . While DR applies when g 1 and g 2 • A 2 are simple, ADMM is a better alternative whereas both g 1 • -A * 2 and g * 2 are simple. Extension to the case K > 2 was proposed for instance in [START_REF] Eckstein | Parallel alternating direction multiplier decomposition of convex programs[END_REF].

• Dykstra algorithm [START_REF] Dykstra | An algorithm for restricted least squares regression[END_REF]. It is able to solve the case where h(x) = xy 2 , A k = Id and the g k are simple functions. It was initially introduced by [START_REF] Dykstra | An algorithm for restricted least squares regression[END_REF]in the case where the g k are indicator functions of closed convex sets, and is generalized in [START_REF] Bauschke | A dykstra-like algorithm for two monotone operators[END_REF] to arbitrary convex functions. It is also extended in [START_REF] Censor | The dykstra algorithm with bregman projections[END_REF][START_REF] Bauschke | Dykstras algorithm with bregman projections: A convergence proof[END_REF] to the case where h is a Bregman divergence. • Primal-Dual schemes. Recently, primal-dual splitting algorithms have been proposed to minimize (6.2) in its full generality, and even more complex objectives, see for instance [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF]223,203,[START_REF] Briceño Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, lipschitzian, and parallel-sum type monotone operators[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms[END_REF]239]. Primal-dual schemes can be used to solve (P y,λ ), (P y,0 ), (P 1 y,ε ) or (P 2 y,γ ).

Finite Model Identification with Forward Backward

The FB algorithm is a good candidate to solve (P y,λ ) when J is simple. Starting from some x (0) ∈ R N , the FB iteration applied to (P y,λ ) reads

x (n+1) = Prox τ n λ J x (n) + τ n Φ * (y -Φx (n) ) ,
where the step-size sequence should satisfy 0 < τ τ n τ < 2/ Φ 2 to ensure convergence of of the sequence x (n) to a minimizer of (P y,λ ).

In fact, owing to partial smoothness of J, much more can be said about the iterates of the FB algorithm. More precisely, after a finite number of iterations, Forward-Backward algorithm correctly identifies the manifold M . This is made formal in the following theorem whose proof can be found in [151].

Theorem 7. Under the assumptions of Theorem 2, x (n) ∈ M for n large enough.

This result sheds some light on the convergence behavior of this algorithm in the favorable case where condition [START_REF] Bauschke | Dykstras algorithm with bregman projections: A convergence proof[END_REF] holds and ( w /λ , λ ) are sufficiently small. In fact, it is shown in [151] that FB identifies in finite time the manifold of any non-degenerate minimizer x ⋆ . As a corollary, if condition [START_REF] Bauschke | Dykstras algorithm with bregman projections: A convergence proof[END_REF] holds at x 0 and ( w /λ , λ ) are sufficiently small, then we recover Theorem 7. These results shed light on the typical convergence behavior of FB observed in such circumstances (e.g. in compressed sensing problems).

Remark 19 (Local linear convergence). The FB generally exhibits a global sublinear O(1/n) convergence rate in terms of the objective function. However, under partial smoothness of J, it is shown in [151] that once the active manifold is identified, the FB algorithm enters a local linear convergence regime (Q-linear in general and Rlinear if M is a linear manifold), whose rate can be characterized precisely in terms of the condition number of Φ T x 0 .

Related Works

Finite support identification and local R-linear convergence of FB to solve (P y,λ ) is established in [START_REF] Bredies | Linear convergence of iterative soft-thresholding[END_REF] under either a very restrictive injectivity assumption, or a nondegeneracy assumption that is a specialization of ours to the ℓ 1 norm. A similar result is proved in [125]. The ℓ 1 norm is a partly smooth function and is therefore covered by Theorem 7. [170] proved Q-linear convergence of FB to solve (P y,λ ) with a data fidelity satisfying restricted smoothness and strong convexity assumptions, and J a so-called convex decomposable regularizer. Again, the latter falls within the class of partly smooth functions, and their result is then subsumed by our analysis.

For general programs, a variety of algorithms, such as proximal and projectedgradient schemes were observed to have the finite identification property of the active manifold. In [132,130], the authors have shown finite identification of manifolds associated to partly smooth functions via the (sub)gradient projection method, Newton-like methods, and the proximal point algorithm. Their work extends that of e.g. [242] on identifiable surfaces from the smooth constrained convex case to a general non-smooth setting. Using these results, [131] considered the algorithm [224] to solve (6.2) when h is C 2 , K = 1, A 1 = Id, and g 1 is simple and partly smooth, but not necessarily convex, and proved finite identification of the active manifold. However, the convergence rates remain an open problem in all these works.

Summary and Perspectives

In this chapter, we have reviewed work covering a large body of literature on the regularization of linear inverse problems. We also showed how these previous works can be all seen as particular instances of a unified framework, namely sensitivity analysis for minimization of convex partly smooth functions. We believe this general framework is the one that should be adopted as long as one is interested in studying fine properties and guarantees of these regularizers, and in particular when the stability of the low-complexity manifold associated to the data to recover is at stake. This analysis is however only the tip of the iceberg, and there is actually a flurry of open problems to go beyond the theoretical results presented in this chapter. We list here a few ones that we believe are important avenues for future works:

• Non-convexity and/or non-finiteness: in this chapter, for the sake of simplicity, we focused on smooth convex fidelity terms and finite-valued convex regularizers.

All the results stated in this chapter extend readily to proper lower semicontinuous convex regularizers, since any such a function is subdifferentially regular. Generalizations of some of the results to non-convex regularizers is possible as well, though some regularity assumptions are needed. This is of practical importance to deal with settings where Φ is not a linear operator, or to impose more agressive regularization (for instance when using ℓ p functional with 0 p < 1 instead of the ℓ 1 norm). There are however many difficulties to tackle in this case. For instance, regularity properties that hold automatically for the convex case have to be either imposed or proved. Another major bottleneck is that some of the results presented here, if extended verbatim, will only assess the recovery of a stationary/critical point. The latter is not a local minimum in general, and even less global. • Dictionary learning: a related non-convex sensitivity analysis problem is to understand the recovery of the dictionary D in synthesis regularization (as defined in Section 2.3.4) when solving problems of the form min

{α k } k ,D∈D ∑ k 1 2 y -ΦDα k 2 + λ J 0 (α k )
where the (y k ) k are a set of input exemplars and D stands for the set of constraints imposed on the dictionary to avoid trivial solutions. Such a non-convex variational problem is popular to compute adapted dictionaries, in particular when J 0 = • 1 , see [START_REF] Elad | Sparse and redundant representations: from theory to applications in signal and image processing[END_REF] and references therein. Although the dictionary learning problem has been extensively studied when J 0 = • 1 , most of the methods lack of theoretical guarantees. The theory of dictionary learning is only beginning to develop, see e.g. [122, 140, 204, 1]. Tackling other regularizers, including analysis ℓ 1 of the form J = J 0 • D * is even more difficult, see e.g. [183,[START_REF] Chen | Learning ℓ 1 -based analysis and synthesis sparsity priors using bi-level optimization[END_REF] for some computational schemes. • Infinite dimensional problems: we dealt in this chapter with finite-dimensional vector spaces. It is not straightforward to extend these results to infinite-dimensional cases. As far as ℓ 2 -stability is concerned, the constants involved in the upper bounds depend on the dimension N, and the scaling might diverge as N → +∞. We refer to Section 3.3 for previous works on convergence rates of Tikhonov regularization in infinite-dimensional Hilbert or Banach spaces. Extending Theorem 2 for possibly non-reflexive Banach spaces is however still out of reach (non-reflexivity is a typical degeneracy when considering lowcomplexity regularization). There exists however some extensions of classical stability results over spaces of measures, such as weak convergence [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], exact recovery [START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF][START_REF] Castro | Exact reconstruction using beurling minimal extrapolation[END_REF] and stable support recovery [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]. • Compressed sensing: as highlighted in Sections 3.3.3 and 4.5.3, the general machinery of partly smooth regularizers (and the associated dual certificates) is well adapted to derive optimal recovery bounds for compressed sensing. Unfortunately, this analysis has been for now only applied to norms ( • 1 , • 1,B ,

• * and • ∞ ). Extending this framework for synthesis and analysis regularizers (see Sections 2.3.4 and 2.3.5) is a difficult open problem.

• Convergence and acceleration of the optimization schemes: Section 6.3 showed how partial smoothness can be used to achieve exact manifold identification after a finite number of iterations using the FB algorithm. This in turn implies a local linear convergence of the iterates, and raises the hope of acceleration using either first-order or second-order information fo the function along the identified manifold (in which we recall it is C 2 ). Studying such accelerations and their guarantees as well as extending this idea other proximal splitting schemes is thus of practical importance to tackle more complicated problems such as e.g. (P y,0 ), (P 1 y,ε ) or (P 2 y,γ ).

Acknowledgements

This work has been supported by the European Research Council (ERC project SIGMA-Vision). We would like to thank our collaborators Charles Deledalle, Charles Dossal, Mohammad Golbabaee and Vincent Duval who have helped to build this unified view of the field.