Charles-Alban Deledalle

Samuel Vaiter

Jalal Fadili

Gabriel Peyré

Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection *

Keywords: Inverse problem, risk estimation, SURE, parameter selection, proximal splitting, sparsity

. Under mild assumptions on the estimation mapping, we show that this approximation is a weakly differentiable function of the parameters and its weak gradient, coined the Stein Unbiased GrAdient estimator of the Risk (SUGAR), provides an asymptotically (with respect to the data dimension) unbiased estimate of the gradient of the risk. Moreover, in the particular case of softthresholding, it is proved to be also a consistent estimator. This gradient estimate can then be used as a basis to perform a quasi-Newton optimization. The computation of the SUGAR relies on the closed-form (weak) differentiation of the non-smooth function. We provide its expression for a large class of iterative methods including proximal splitting ones and apply our strategy to regularizations involving non-smooth convex structured penalties. Illustrations on various image restoration and matrix completion problems are given.

Introduction

In this paper, we consider the recovery problem of a signal x 0 ∈ X (where X = R N or is a suitable finite-dimensional Hilbert space that can be identified to R N ) from a realization y ∈ Y = R P of the normal random vector

Y = µ 0 + W with µ 0 = Φx 0 (1)
where W ∼ N (0, σ 2 Id P ), and the linear imaging operator Φ : X → Y entails some loss of information. Typically, P = dim(Y) is smaller than N = dim(X ), or Φ is rank-deficient, and the recovery problem is ill-posed. Let (y, θ) → x(y, θ) be some recovery mapping, possibly multivalued, which attempts to approach x 0 from a given realization y ∈ Y of Y and is parametrized by a collection of continuous parameters θ ∈ Θ. Throughout, Θ is considered as a subset of a linear subspace of dimension dim(Θ). We also denote µ(y, θ) = Φx(y, θ) ∈ Y and assume in the rest of the paper that it is always a single-valued mapping though x(y, θ) may not.

Depending on the smoothness of the mapping y → µ(y, θ), the recovered estimate enjoys different regularity properties. For instance, µ(y, θ) can be built by solving a variational problem with some regularizing penalty parametrized by θ (see the example in [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], as well as Section 4 and 5). This regularization is generally chosen so as to preserve/promote the interesting structure underlying x 0 , e.g. singularities, textures, etc. Also, depending on its choice and that of the data fidelity, the resulting mapping y → µ(y, θ) may be smooth or not. To cover most of these situations, throughout the paper, we will assume that (y, θ) → µ(y, θ) is weakly differentiable with respect to both the observation y, and the collection of parameters θ.

Recall that for a locally integrable function f : a ∈ Ω → R, Ω is an open subset of R N , its weak partial derivative with respect to a i in Ω is the locally integrable function g i on Ω such that

Ω g i (a)ϕ(a)da = - Ω f (a)
∂ϕ(a) ∂a i da holds for all functions ϕ ∈ C 1 c (Ω), i.e. the space of continuously differentiable functions of compact support. The weak partial derivative, if it exists, is uniquely defined Lebesguealmost everywhere (a.e.). Thus we write

g i = ∂f ∂a i
and all such pointwise relations involving weak derivatives will be accordingly understood to hold Lebesgue-a.e. A function is said to be weakly differentiable if all its weak partial derivatives exist. Similarly, a vector-valued function h : a ∈ R N → h(a) = (h 1 (a), . . . , h P (a)) ∈ R P is weakly differentiable if h k (a) is weakly differentiable ∀k ∈ {1, . . . , P }, and we will denote ∂h(a) its weak Jacobian, and ∇g(a) = ∂h(a) * its adjoint. Remark that weak differentiation concepts boil down to the classical ones when the considered function is C 1 . A comprehensive account on weak differentiability can be found in e.g. [START_REF] Evans | Measure theory and fine properties of functions[END_REF][START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Orde[END_REF].

Getting back to the estimator x(y, θ), we now discuss some typical examples covered in this paper.

• Given (y, θ), consider a minimizer of a convex variational problem of the form x(y, θ) = Argmin x∈X {E(x, y, θ) = H(y, Φx) + R(x, θ)} [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] where x(y, θ) is the set of minimizers of x → E(x, y, θ) which is considered nonempty (the minimizer may not be unique but is assumed to exist). The data fidelity term x → H(y, Φx) is defined using a strongly convex map µ → H(y, µ). The regularization term x → R(x, θ) is assumed to be a closed proper and convex function, that accounts for the prior structure of x 0 . Typical priors correspond to non-smooth regularizers such as sparsity in a suitable domain, e.g. Fourier, wavelet [START_REF] Mallat | A wavelet tour of signal processing[END_REF], or gradient [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Such regularizers are usually parametrized with a collection of parameters θ. A typical example is R(x, θ) = θR 0 (x) where θ ∈ R + is a scaling which controls the strength of the regularization. Of course, more complicated (multi-parameters) regularizations, are often considered in the applications, and our methodology aims at dealing with these higher dimensional sets of parameters.

An important observation is that even though x(y, θ) may not be a singleton (minimizer of E(x, y, θ) may not be unique), strict convexity of H(y, •) implies that all minimizers share the same image under Φ, see e.g. [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF]. Hence (y, θ) → µ(y, θ) is defined without ambiguity as a single-valued mapping. Moreover, strong convexity of H(y, •) implies that y → µ(y, θ) is non-expansive (i.e. uniformly 1-Lipschitz) [START_REF] Vaiter | The degrees of freedom of partly smooth regularizers[END_REF], hence weakly differentiable [28, Theorem 5, Section 4.2.3].

• Consider now the ℓ-th iterate, denoted by x (ℓ) (y, θ), of an iterative algorithm converging to a fixed point of an operator acting on X . In this case, θ can include the parameters of the fixed point operator, as well as other continuous parameters inherent to the fixed point iteration (such as, e.g., step sizes). Section 4 is completely dedicated to this setting, and appropriate sufficient conditions will be exhibited to ensure weak differentiability of x (ℓ) (y, θ) with respect to both its arguments.

This general setting encompasses the case of proximal splitting methods that have become popular to solve large-scale optimization problems of the form (2), especially with convex non-smooth terms, e.g. those encountered in sparsity regularization. The precise splitting algorithm to be used depends on the structure of the optimization problem at hand. See for instance [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF]Fixed-Point Algorithms for Inverse Problems in Science and Engineering[END_REF] for an overview. Some of these algorithms are considered in detail in Section 4.

The choice of θ is generally a challenging task, especially as the dimension of Θ gets large. Ideally, one would like to choose the parameters θ ⋆ that makes µ(y, θ ⋆ ) (or some appropriate image of it) as faithful as possible to µ 0 (or some appropriate image of it). Formally, this can be cast as selecting θ ⋆ that minimizes the expected reconstruction error (a.k.a., mean-squared error or quadratic risk), i.e.

θ ⋆ ∈ Argmin θ∈Θ {R A {µ}(µ 0 , θ) = E W A(µ(Y, θ) -µ 0 ) 2 } ( 3
)
where the matrix A ∈ R M ×P is typically chosen to counterbalance the effect of Φ, see Section 2.1 for a precise discussion. If θ → R A {µ}(µ 0 , θ) were sufficiently smooth, at least locally (e.g. Lipschitz), one could expect to solve (3) using a (sub)gradient-descent scheme relying on the (weak) gradient of the risk ∇ 2 {R A {µ}}(µ 0 , θ), where the subscript 2 specifies that the (weak) gradient is with respect to the second argument θ. However, this would only apply if µ 0 were available. In the context of our observation model [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF], µ 0 is however considered to be unknown. Our motivation is then to build an estimator of ∇ 2 {R A {µ}}(µ 0 , θ) that depends solely on y, without prior knowledge of µ 0 .

Toward this goal, we adopt the framework of the (generalized) Stein Unbiased Risk Estimator (SURE) [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF][START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF][START_REF]Least squares estimation without priors or supervision[END_REF][START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF][START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF]. For a fixed θ, the celebrated Stein's lemma [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] allows to unbiasedly estimate R A {µ}(µ 0 , θ) through the weak Jacobian ∂ 1 µ(y, θ), where the subscript 1 specifies that the (weak) Jacobian is with respect to the first argument y. This idea have been exploited for years in several statistical and signal processing applications, typically for selecting thresholds in wavelet based reconstruction algorithms, see e.g. [START_REF] Benazza-Benyahia | Building robust wavelet estimators for multicomponent images using stein's principle[END_REF]5,[START_REF] Chaux | A nonlinear stein-based estimator for multichannel image denoising[END_REF][START_REF] Deledalle | Stein consistent risk estimator (SCORE) for hard thresholding[END_REF][START_REF] Donoho | Adapting to Unknown Smoothness Via Wavelet Shrinkage[END_REF][START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF][START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and MRI reconstruction using GCV and SURE-based methods[END_REF][START_REF] Raphan | Optimal denoising in redundant representations[END_REF]. Given such an estimator R A {µ}(y, θ) (see Section 2.1), the idea is so to replace the optimization problem [START_REF] Benazza-Benyahia | Building robust wavelet estimators for multicomponent images using stein's principle[END_REF] with

θ ⋆ ∈ Argmin θ∈Θ R A {µ}(y, θ) . (4)
Providing that the variance of 1 P R A {µ}(Y, θ) can be made arbitrarily small or even asymptotically vanishing as P increases, so that it becomes a consistent estimator of 1 P R A {µ}(µ 0 , θ), one can expect that the minimizers of (4) becomes close to those of [START_REF] Benazza-Benyahia | Building robust wavelet estimators for multicomponent images using stein's principle[END_REF].

It remains to find an efficient way to solve the optimization (4). Again, a (sub)gradientdescent algorithm can qualify as a good candidate if θ → R A {µ}(y, θ) were sufficiently smooth. To our knowledge, only [START_REF] Deledalle | Poisson NL means: Unsupervised non local means for Poisson noise[END_REF] have performed such an optimization with Newton's method where (y, θ) → R A {µ}(y, θ) was C ∞ . Unfortunately, being a function of ∂ 1 µ(y, θ), θ → R A {µ}(y, θ) is in general not differentiable, not even continuous (think of a simple soft-thresholding). This then precludes the use of standard descent schemes.

The common practice has been to apply an exhaustive search by evaluating the risk estimate R A {µ}(y, θ) at different values of θ. Even if in some particular cases this can be done efficiently (see for instance [START_REF] Donoho | Adapting to Unknown Smoothness Via Wavelet Shrinkage[END_REF]), the computational expense can become prohibitive in general especially as dim(Θ) increases.

Derivative-free optimization algorithms have also been investigated (see for instance [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF] for the case of 2 parameters). But such approaches typically do not scale up to problems where Θ has a linear vector space structure with dimension larger than 2. Their performance are known to degrade exponentially with problem size, and they require to compute a lower and an upper bound on the optimal value over a given region.

Contributions

In this paper, we address the challenging problem of solving efficiently (4): a main subject of interest for applications that has been barely investigated. Our main contribution (Section 3) is an effective strategy to optimize automatically a collection of parameters θ independently of their dimension. While classical unbiased risk estimates entail optimizing a non-continuous function of the parameters, we show that the biased risk estimator introduced in [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF] is differentiable in the weak sense. This allows us, whenever the derivatives exist, to perform a quasi-Newton optimization driven by a biased estimator of the gradient of the risk based on the evaluation of ∂ 2 µ(y, θ). Such optimization technique can be provably faster thanks to first-order information compared to derivative-free approaches. We prove that, under mild assumptions, this estimator is asymptotically (with respect to P ) unbiased, hence the name: Stein Unbiased GrAdient estimator of the Risk (SUGAR). Moreover, in the particular case of soft-thresholding, we go a step further and show that it is actually a consistent estimator of the gradient of the risk.

As a second contribution (Section 4), we propose a versatile approach to compute the derivatives ∂ 1 µ(y, θ) and ∂ 2 µ(y, θ), involved respectively in the computation of the SURE and SUGAR, when µ(y, θ) is computed through an iterative algorithm, typically proximal splitting methods. We illustrate the versatility of our method by applying it to both primal (forwardbackward [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], Douglas-Rachford [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] and generalized foward-backward [START_REF] Raguet | A generalized forward-backward splitting[END_REF]) and primal-dual [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] algorithms (see [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF] for a recent review). The proposed methodology can however be adapted to any other proximal splitting method and more generally to any algorithm whose iteration operator is weakly differentiable.

Numerical simulations involving multi-parameter selection for image restoration and matrix completion problems are reported in Section 5. The proofs of our results are collected in the appendix.

Overview on Risk Estimation

This section gives an overview of the literature to estimate the risk via the SURE and its variants for ill-posed inverse problems contaminated by additive white Gaussian noise.

Stein Unbiased Risk Estimator

Degrees of freedom (DOF) is often used to quantify the complexity of a statistical modeling procedure, see for instance, GCV (generalized cross-validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]). From [START_REF] Efron | How biased is the apparent error rate of a prediction rule?[END_REF][START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF], the degrees of freedom of a function y → µ(y, θ) relatively to a matrix A ∈ R M ×P is given by

df A {µ}(µ 0 , θ) = P i=1 cov(AY i , (Aµ(Y, θ)) i ) σ 2 , (5)
such that df A {µ}(µ 0 , θ) is maximal when Aµ(Y, θ) is highly correlated with the random vector AY . Taking A = Id, leads to the standard definition of the DOF defined in the seminal work of Efron [START_REF] Efron | How biased is the apparent error rate of a prediction rule?[END_REF]. But other choices of A allow to counterbalance the undesirable effect of the linear operator Φ (recall that µ 0 = Φx 0 ). For instance, setting A = (Φ * Φ) -1 Φ * when Φ has full-rank, or A = Φ * (ΦΦ * ) + when Φ is rank deficient * , provides a measure of the DOF relatively to the least-squares estimate of x 0 , i.e. x LS (y) = Ay [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF][START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF][START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF].

With the proviso that y → µ(y, θ) is weakly differentiable with essentially bounded weak partial derivatives, an unbiased estimate of the DOF can be used to unbiasedly estimate the risk in [START_REF] Benazza-Benyahia | Building robust wavelet estimators for multicomponent images using stein's principle[END_REF]. This leads to the (generalized) SURE (also known as weighted SURE [START_REF] Ramani | Iterative weighted risk estimation for nonlinear image restoration with analysis priors[END_REF]) given as

SURE A {µ}(y, θ) = A(µ(y, θ) -y) 2 -σ 2 tr(A * A) + 2σ 2 df A {µ}(y, θ) (6)
with df A {µ}(y, θ) = tr (A∂ 1 µ(y, θ)A * )

where we recall that ∂ 1 µ(y, θ) is the weak Jacobian of µ(y, θ) with respect to the first argument y. It can be shown that (see e.g. [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF][START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF][START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF])

E W [ df A {µ}(Y, θ)] = df A {µ}(µ 0 , θ) and E W [SURE A {µ}(Y, θ)] = R A {µ}(µ 0 , θ) .
Expression ( 6) is general enough to encompass unbiased estimates of the prediction risk

E W µ(Y, θ) -µ 0 2 (i.e. A = Id), the projection risk E W Π(x(Y, θ) -x 0 ) 2 , where Π
is the orthogonal projector on ker(Φ) ⊥ (i.e. A = Φ * (ΦΦ * ) + ), and the estimation risk

E W x(Y, θ) -x 0 2 when Φ has full rank (i.e. A = (Φ * Φ) -1 Φ * )
. This can prove useful when Φ is rank deficient, since in this case, the minimizers of the prediction risk can be far away from the minimizers of the estimation risk [START_REF] Rice | Choice of smoothing parameter in deconvolution problems[END_REF]. The projection risk restricts the estimate to the subspace where there is a signal beside noise, and in this sense, is a good approximation of the estimation risk [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF].

Note that generalization of the SURE have been developed for other noise models, typically within the multivariate canonical exponential family (see e.g. [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF][START_REF] Hudson | A natural identity for exponential families with applications in multiparameter estimation[END_REF][START_REF] Hudson | Maximum likelihood restoration and choice of smoothing parameter in deconvolution of image data subject to poisson noise[END_REF][START_REF]Least squares estimation without priors or supervision[END_REF]).

Applications of SURE emerged for choosing the smoothing parameters in families of linear estimates [START_REF] Li | From Stein's unbiased risk estimates to the method of generalized cross validation[END_REF] such as for model selection, ridge regression, smoothing splines, etc. After its introduction in the wavelet community with the SURE-Shrink algorithm [START_REF] Donoho | Adapting to Unknown Smoothness Via Wavelet Shrinkage[END_REF], it has been widely used for various image restoration problems, e.g. with sparse regularizations [START_REF] Benazza-Benyahia | Building robust wavelet estimators for multicomponent images using stein's principle[END_REF]5,[START_REF] Cai | A data-driven block thresholding approach to wavelet estimation[END_REF][START_REF] Chaux | A nonlinear stein-based estimator for multichannel image denoising[END_REF][START_REF] Deledalle | Stein consistent risk estimator (SCORE) for hard thresholding[END_REF][START_REF] Luisier | SURE-LET for orthonormal wavelet-domain video denoising[END_REF][START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF][START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF][START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and MRI reconstruction using GCV and SURE-based methods[END_REF][START_REF] Ramani | Iterative weighted risk estimation for nonlinear image restoration with analysis priors[END_REF][START_REF] Raphan | Optimal denoising in redundant representations[END_REF][START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF] or with non-local filters [START_REF] Deledalle | Non-local Methods with Shape-Adaptive Patches (NLM-SAP)[END_REF][START_REF] Duval | A bias-variance approach for the non-local means[END_REF][START_REF] Van De Ville | SURE-based Non-Local Means[END_REF][START_REF]Non-local means with dimensionality reduction and SURE-based parameter selection[END_REF].

However, a major practical difficulty when using the SURE lies in the numerical computation of the DOF estimate, i.e. the quantity df A {µ}(y, θ) for a given realization y. We now give a brief overview of some previous works to deal with this computation.

Closed-form SURE

The SURE is based on a DOF estimate df A {µ}(Y, θ) that can be sampled from the observation y ∈ R P by evaluating the Jacobian ∂ 1 µ(y, θ) ∈ R N ×P . A natural way, to evaluate ∂ 1 µ(y, θ) would be to derive its closed-form expression. This has been studied for some classes of variational problems.

In quadratic regularization (e.g. ridge regression), where solutions are of the form x(y, θ) = K(θ)y, where K(θ) is known as the hat or influence matrix, the Jacobian has a closed-form ∂ 1 µ(y, θ) = ΦK(θ). In ℓ 1 -synthesis regularization (a.k.a. the lasso), the Jacobian matrix depends on the support (set of non-zero coefficients) of any lasso solution x(y, θ). An estimator of the DOF can then be retrieved from the number of non-zero entries of this solution [START_REF] Dossal | The degrees of freedom of the lasso for general design matrix[END_REF][START_REF]Degrees of freedom in lasso problems[END_REF][START_REF] Zou | On the "degrees of freedom" of the lasso[END_REF]. These results have in turn been extended to more general sparsity promoting regularizations [START_REF]Unbiased risk estimation for sparse analysis regularization[END_REF][START_REF] Kato | On the degrees of freedom in shrinkage estimation[END_REF][START_REF] Solo | Threshold selection for group sparsity[END_REF][START_REF] Tibshirani | The solution path of the generalized lasso[END_REF][START_REF]Degrees of freedom in lasso problems[END_REF][START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF][START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], and spectral regularizations (e.g. nuclear norm) [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF][START_REF]Risk estimation for matrix recovery with spectral regularization[END_REF].

This approach however has three major bottlenecks. First, deriving the closed-form expression of the Jacobian is in general challenging and has to be addressed on a case by case basis. Second, in large-dimensional problems, evaluating numerically this Jacobian is barely possible. Even if it were possible, it might be subject to serious numerical instabilities. Indeed, solutions of variational problems are achieved via iterative schemes providing iterates x (ℓ) (y, θ) that eventually converge to the set of solutions as ℓ → +∞. And yet, for instance, substituting the support of the true solution by the support of x (ℓ) (y, θ), obtained at a prescribed convergence accuracy, might be imprecise (all the more since the problem is ill-conditioned).

The three next sections review previous work to address one or some of these three points.

Monte-Carlo SURE

To deal with the large dimension of the Jacobian, the standard approach is to exploit the fact that the DOF only depends on the trace of A∂ 1 µ(y, θ)A * . In denoising applications where A = Id and Φ = Id, this trace can generally be obtained by closed-form computations of the P diagonal elements of ∂ 1 µ(y, θ) (see e.g. [START_REF] Donoho | Adapting to Unknown Smoothness Via Wavelet Shrinkage[END_REF][START_REF] Van De Ville | SURE-based Non-Local Means[END_REF]). This can also be done for some particular inverse problems. For instance, the authors of [START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF] provide an expression of this trace for the wavelet-vaguelette estimator when Φ is a convolution matrix and A = Φ + . However, in more general settings, the complexity of the closed-form computation of the trace is non-linear, typically the number of operations is in O(P × P ) (think of Φ a mixing operator or µ an iterative estimator). To avoid such a costly procedure, the authors of [START_REF] Girard | A fast Monte-Carlo cross-validation procedure for large least squares problems with noisy data[END_REF][START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF] suggest making use of the following trace equality

df A {µ}(y, θ) = tr (A∂ 1 µ(y, θ)A * ) = E ∆ ∂ 1 µ(y, θ)[∆], A * A∆ (7)
where ∆ ∼ N (0, Id P ) and ∂ 1 µ(y, θ)[δ] ∈ R P denotes the directional derivative of y → µ(y, θ) at y in direction δ. Remark that ∆ does not necessary have to be Gaussian and higher precisions can be reached in some specific cases, see for instance [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF][START_REF] Dong | Stochastic estimation with Z 2 noise[END_REF][START_REF] Hutchinson | A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines[END_REF][START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF]. As shown in [START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF], the performance of this trace estimator is governed by the distribution of the singular values of the operator A∂ 1 µ(y, θ)A * . More specifically, the slower the decay, the better the performance. While it is difficult to make a general claim, we observed numerically that for the recovery problems we consider, it provide a very accurate estimator of the trace. Hence, following [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF][START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF], an estimate of SURE A {µ}(y, θ) can be obtained by Monte-Carlo simulations using

SURE A MC {µ}(y, θ, δ) = A(µ(y, θ) -y) 2 -σ 2 tr(A * A) + 2σ 2 df A MC {µ}(y, θ, δ) with df A MC {µ}(y, θ, δ) = ∂ 1 µ(y, θ)[δ], A * Aδ . (8)
The evaluation of (8) necessitates only computing the P entries of ∂ 1 µ(y, θ) [δ].

It remains to find a stable and efficient way to evaluate for any vector δ ∈ R P the directional derivative ∂ 1 µ(y, θ)[δ] ∈ R P .

Iterative Differentiation for Monte-Carlo SURE

When considering solutions x(y, θ) of a variational problem, the DOF cannot be robustly estimated if one knows only the iterates µ (ℓ) (y, θ) that eventually converge to some µ(y, θ) as ℓ → +∞. It appears then natural to estimate the DOF of µ (ℓ) (Y, θ) directly and make the assumption that it will converge to that of µ(y, θ). For a realization y ∈ R P , one can sample an estimate of the DOF of the iterate µ (ℓ) (Y, θ) by evaluating its directional derivative ∂ 1 µ (ℓ) (y, θ) [δ]. A practical way, initiated by [START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF], to compute this quantity, consists in recursively differentiating the sequence of iterates. The authors of [START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF] have derived the closedform expression of the directional derivative for the Forward-Backward (FB) algorithm. The directional derivative at iteration ℓ + 1, denoting by D

(ℓ+1) µ = ∂ 1 µ (ℓ+1) (y, θ)[δ], is obtained iteratively as a function of µ (ℓ) (y, θ) and D (ℓ) µ = ∂ 1 µ (ℓ) (y, θ)[δ].
The Monte-Carlo DOF and the Monte-Carlo SURE can in turn be iteratively estimated by plugging

∂ 1 µ (ℓ) (y, θ)[δ] in (8) leading to SURE A MC {µ (ℓ) }(y, θ, δ) = A(µ (ℓ) (y, θ) -y) 2 -σ 2 tr(A * A) + 2σ 2 df A MC {µ (ℓ) }(y, θ, δ) with df A MC {µ (ℓ) }(y, θ, δ) = D (ℓ) µ , A * Aδ . (9)
A similar approach is described in [START_REF] Giryes | The projected GSURE for automatic parameter tuning in iterative shrinkage methods[END_REF]. Pursuing this idea, the authors of [START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and MRI reconstruction using GCV and SURE-based methods[END_REF][START_REF] Ramani | Iterative weighted risk estimation for nonlinear image restoration with analysis priors[END_REF] have recently provided such closed-form expressions in the case of the split Bregman method. Concurrently, in an early short version of this paper [START_REF] Deledalle | Proximal splitting derivatives for risk estimation[END_REF], we have also considered this approach for general proximal splitting algorithms, an approach that we extend in Section 4.

Finite-Difference SURE

An alternative initiated in [START_REF] Shen | Adaptive model selection[END_REF][START_REF] Ye | On measuring and correcting the effects of data mining and model selection[END_REF] and rediscovered in [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF] consists in estimating tr (A∂ 1 µ(y, θ)A * ) via finite differences given, for ε > 0, by

tr (A∂ 1 µ(y, θ)A * ) ≈ P i=1 [A * A(µ(y + εe i , θ) -µ(y, θ))] i ε , (10)
where (e i ) 1 i P is the canonical basis of R P . Plugging this expression in (8) yields the Finite-Difference (FD) SURE A given by

SURE A FD {µ}(y, θ, ε) = A(µ(y, θ) -y) 2 -σ 2 tr(A * A) + 2σ 2 df A FD {µ}(y, θ, ε) with df A FD {µ}(y, θ, ε) = 1 ε P i=1 (A * A(µ(y + εe i , θ) -µ(y, θ))) i . (11)
The main advantage of this method is that (y, θ) → µ(y, θ) can be used as a black-box, i.e., without knowledge on the underlying algorithm that provides µ(y, θ), while, for ε small enough, it performs as well as the approach described in Section 2.4 that requires the knowledge of the derivatives in closed-form. In fact, if y → µ(y, θ) is Lipschitz-continuous, then it is differentiable Lebesgue a.e. (Rademacher's theorem), and its derivative equals its weak derivative Lebesgue a.e. [28, Theorem 1-2, Section 6.2], which in turn implies

lim ε→0 SURE A FD {µ}(y, θ, ε) = SURE A {µ}(y, θ) Lebesgue a.e. ( 12
)
The value of ε can so be chosen as small as possible as soon as it does not rise to numerical instabilities due to limited machine precision. To avoid numerical instabilities, ε should be chosen in a reasonable range of values with respect to the amplitudes of the data. In practice, we observe that precised results can be reached compared to the closed-form derivation, hence yielding to a quasi unbiased risk estimator (i.e., with a negligible bias). It remains that when the data dimension P is large, the evaluation of P finite differences along each axis might be numerically intractable. In that case, the Monte-Carlo approach (see Section 2.3) can also be used in conjunction with finite differences leading to the Finite-Difference Monte-Carlo (FDMC) SURE A given by

SURE A FDMC {µ}(y, θ, δ, ε) = A(µ(y, θ) -y) 2 -σ 2 tr(A * A) + 2σ 2 df A FDMC {µ}(y, θ, δ, ε) with df A FDMC {µ}(y, θ, δ, ε) = 1 ε µ(y + εδ, θ) -µ(y, θ), A * Aδ . ( 13
)
The originality of our approach described in the next section is to devise a grounded choice of ε > 0. This introduces a bias in the estimation of the risk. Nevertheless, as we will see, using ε > 0 plays an important role in risk optimization since, unlike SURE A {µ}, SURE A FD {µ} is a smooth function of θ in the weak sense. This is the key point to optimize the risk. By choosing ε > 0 carefully, a smoother objective function can be used as a basis to perform a quasi-Newton-like optimization at the expense of a controlled bias.

Risk Estimate Minimization

In this section, we investigate how risk estimates can be used for optimizing a collection of continuous parameters.

Stein's Unbiased GrAdient Risk (SUGAR) Estimator

The difficulty is that even if θ → R A {µ}(µ 0 , θ) is differentiable in the weak sense, the function θ → SURE A {µ}(y, θ) might contain discontinuities. Typically, df

A {µ}(y, θ) has discontinuities where (y, θ) → µ(y, θ) is not differentiable. We start with a simple result showing that unlike SURE A {µ}(y, θ), the finite-difference based mapping θ → SURE A FD {µ}(y, θ, ε), for ε > 0, is weakly differentiable. Proposition 1. Assume µ(y, θ) is weakly differentiable with respect to y and θ. Given ε > 0, df A FD {µ}(y, θ, ε) and SURE A FD {µ}(y, θ, ε) are also weakly differentiable with respect to y and θ, and their (weak) gradients with respect to θ are given, for almost all θ ∈ Θ, as

SUGAR A FD {µ}(y, θ, ε) = ∇ 2 {SURE A FD {µ}}(y, θ, ε) = 2∂ 2 µ(y, θ) * A * A(µ(y, θ) -y) + 2σ 2 ∇ 2 { df A FD {µ}}(y, θ, ε) where ∇ 2 { df A FD {µ}}(y, θ, ε) = 1 ε P i=1 (∂ 2 µ(y + εe i , θ) -∂ 2 µ(y, θ)) * A * Ae i .
Thanks to Proposition 1, a quasi-Newton-like method can now be used to optimize SURE A FD {µ}(y, θ, ε) for the vector of continuous parameters θ by implementing the iteration

θ n+1 = θ n -B n SUGAR A FD {µ}(y, θ n , ε)
where ) is a sequence of definite-positive matrices. Typically, if θ → SURE A FD {µ}(y, θ, δ, ε) behaves locally as a C 2 function, B n should approach the inverse of the corresponding Hessian at θ n . Remark that in general there is no guarantee that the risk has a unique global minimizer, though in the 1-D case it is generally the case. When several parameters are involved, such objective can have several local minima and specific quasi-Newton-like methods might be developed to avoid being stuck in one of them.

B n ∈ R dim(Θ)×dim(Θ
In practice, the calculation of SUGAR A FD depends on the computation of the Jacobian matrices with respect to the parameters θ. We will see in Section 4 how this quantity can be efficiently computed when µ results from an iterative algorithm.

We now turn to the asymptotic unbiasedness of the proposed gradient estimator of the risk as ε approaches 0. Toward this goal we need the following assumptions.

(A.1) The mapping y → µ(y, θ) is uniformly Lipschitz continuous with Lipschitz constant L 1 .

(A.2) The mapping y → µ(y, θ) is such that µ(0, θ) = 0 for any θ.

(A.3) The mapping θ → µ(y, θ) is uniformly Lipschitz continuous with Lipschitz constant L 2 independently of y.
Remark 1 (Discussion of the assumptions).

1. Assumption (A.1) is mild and is fulfilled in many situations of interest. In particular, this is the case when y → µ(y, θ) is the proximal operator of proper closed and convex function, as considered in Section 4 (see also Section 3.2 for soft-thresholding). Standard convex analysis arguments [START_REF] Hiriart-Urruty | Convex Analysis And Minimization Algorithms[END_REF] show that the proximity operator is indeed a uniformly Lipschitz of its argument y with constant L 1 = 1, independently of θ.

Assumption (A.2) is very natural and does not entail any loss of generality. It basically

states that, when the observations are zero, so is the estimator.

3. As far as Assumption (A.3) is concerned, it is verified under certain circumstances. This is for instance the case when µ(y, θ) = Prox θG (y), θ > 0, where G is the gauge (see Definition 2 in Appendix B) of any compact convex set containing the origin as an interior point † ; see Proposition 5 in Appendix B. By induction, this also holds when

µ(y, θ) = Prox θ 1 G 1 • • • • • Prox θmGm (y), θ ∈]0, +∞[ m ,
and for any i = 1, . . . , m, G i is the gauge of any compact convex set containing the origin as an interior point, see Corollary 5. Typical instances of these gauges are norms, e.g. ℓ 1 , ℓ 1ℓ 2 or nuclear norms very popular now in the signal and image processing community.

4. Assumption (A.3) can be relaxed to cover the case where L 2 depends on y. In such a situation, additional assumptions on the function y → L 2 (y) are needed for steps 3) and 4) in the proof of Theorem 1 to go through. We omit this case for the sake of clarity and to avoid further technicalities.

We are now ready to state our theorem.

Theorem 1 (Asymptotic unbiasedness of SUGAR). Assume that (A.1)-(A.3) hold. Then, R A {µ}(µ 0 , •) and df A {µ}(µ 0 , •) are weakly differentiable, and for any Lebesgue point θ,

lim ε→0 E W SUGAR A FD {µ}(Y, θ, ε) = ∇ 2 {R A {µ}}(µ 0 , θ) and lim ε→0 E W ∇ 2 { df A FD {µ}}(Y, θ, ε) = ∇ 2 {df A {µ}}(µ 0 , θ) .
Theorem 1 can be given the following interpretation. As ε gets close to 0, e.g. a decreasing function of the dimension P ‡ , the gradient of SURE A FD {µ}(y, •, ε) (normalized by P ) can be used to estimate the gradient of the risk (also normalized by P ) provided that P is large enough.

However, even if ε should decrease towards 0, it should not decrease too fast. In particular, for a fixed dimension P , the step ε cannot be chosen arbitrarily small. This would not be an issue if µ(y, •) were differentiable, but in general, there might be singularities. In fact, for a finite dimension P , the limit when ε → 0 of the sample SUGAR A FD {µ}(y, θ, ε) may not even exist, though that of its expectation does exist Lebesgue a.e. as shown in the proof of Theorem 1. As a consequence, the quantity 1 P SUGAR A FD {µ}(Y, θ, ε) can become very unstable when ε decreases too fast with the dimension P . The underlying statistical question is whether one can control the variance of 1 P SUGAR A FD {µ}(Y, θ, ε) as P increases, and make arbitrarily small or even asymptotically vanishing, so that 1 P SUGAR A FD {µ}(Y, θ, ε) becomes a consistent estimator. Unfortunately, consistency of our gradient estimator of the risk is very intricate to get in the general case, as it is the case for the consistency of the SURE. However, when µ specializes to soft-thresholding, such a result can be achieved.

SUGAR for Soft-tresholding

In this section, we show that the proposed gradient estimator of the risk can be consistent in the case where µ is the soft-thresholding (ST) function and A = Id P . The ST is the proximal operator of the ℓ 1 -norm. Understanding the ST is of chief interest since it is at the heart of any proximal splitting algorithm solving a regularized inverse problem involving terms of the form D * x 1 where D is a linear operator.

Let first recall the definition of soft-thresholding.

Definition 1 (Soft-Thresholding). The soft-thresholding (ST) is defined, for λ > 0, and for all 1 i P , as

ST(y, λ) i =    y i + λ if y i -λ 0 if -λ < y i < λ y i -λ otherwise . ( 14
)
Observe that as a proximity operator of a norm, soft-thresholding satisfies Assumptions (A.1) through (A.3) of Theorem 1, see the corresponding discussion. Hence, we already anticipate from Theorem 1 that our gradient estimator of the soft-thresholding risk is asymptotically unbiased.

We start with following lemma which collects the statistics of the gradient of the finite difference DOF estimator. ‡ as we will see, the higher the dimension P , the smaller ε could be.

Lemma 1 (Statistics of the gradient of the finite difference DOF estimator).

Let 0 < ε < 2λ. The weak gradient of λ → df FD {ST}(Y, λ, ε) is such that E W ∇ 2 { df FD {ST}}(Y, λ, ε) = -1 2 P i=1 ϕ[(µ 0 ) i , λ, ε] ε ,
and

V W ∇ 2 { df FD {ST}}(Y, λ, ε) = 1 2ε P i=1 ϕ[(µ 0 ) i , λ, ε] ε - 1 4 P i=1 ϕ[(µ 0 ) i , λ, ε] ε 2 .
where for

a ∈ R, ϕ[a, λ, ε] = erf a+λ+ε √ 2σ -erf a+λ √ 2σ + erf a-λ+ε √ 2σ -erf a-λ √ 2σ .
We now turn to the asymptotic behavior of the proposed gradient estimator of the risk for large P , at a single realization of Y , i.e., our observation y. To this end, we first have to define how the observation model evolves with the dimension P . Given z 0 ∈ R N , we consider the sequence {Ψ P } P 1 where Ψ P ∈ R P ×N and such that, for all P > 1, Ψ P is the sub-matrix obtained by cutting down one line of Ψ P +1 . We can then define a sequence of observation models as a the sequence of random vector {Y P } P 1 defined as

Y P = Ψ P z 0 + W P where W P ∼ N (0, σ 2 Id P ). ( 15
)
We also define the sequence {(µ 0 ) P } P 1 where (µ 0 ) P = Ψ P z 0 . In the following, for the sake of clarity, we omit the dependency of Y P , W P and (µ 0 ) P on P .

We can now state our consistency result for soft-thresholding.

Theorem 2 (Consistency of SUGAR). Take ε(P ) such that lim P →∞ ε(P ) = 0 and lim P →∞ P -1 ε(P ) -1 = 0. Then for any Lebesgue point λ > 0 (i.e. such that ∀(i, P ), λ = |Y i | and λ = |Y i + ε(P )e i |)

plim P →∞ 1 P (SUGAR FD {ST}(Y, λ, ε(P )) -∇ 2 {R{ST}}(µ 0 , λ)) = 0
and plim

P →∞ 1 P ∇ 2 { df FD {ST}}(Y, λ, ε(P )) -∇ 2 {df {ST}}(µ 0 , λ) = 0 .
In plain words, Theorem 2 asserts that for our gradient estimator of the soft-thresholding risk to be consistent, ε(P ) should not decrease faster than the inverse of the dimension P . With the proviso that ε(P ) fulfills the requirement, for P large enough, 1 P SUGAR FD {ST}(y, λ, ε(P )) is guaranteed to come close to 1 P ∇ 2 {R{ST}}(µ 0 , λ) with high probability. Unfortunately, Theorem 2 does not dictate an explicit choice of ε(P ), and the practitioner may wonder how to choose this value for a given P . It turns out that studying the mean squared error (MSE) of the gradient of the finite difference DOF estimator helps unveiling the link between P and ε through a bias-variance trade-off.

Proposition 2 (MSE of the gradient of the finite difference DOF estimator). The weak gradient of λ → df {ST}(Y, λ, ε) is such that where the statistics of ∇ 2 { df {ST}}(Y, λ, ε) are given in Lemma 1 and

E W 1 P ∇ 2 { df {ST}}(Y, λ, ε) -∇ 2 {df {ST}}(µ 0 , λ) 2 = 1 P 2 E W ∇ 2 { df {ST}}(Y, λ, ε) -∇ 2 {df {ST}}(µ 0 , λ) 2 Bias 2 + 1 P 2 V W ∇ 2 { df {ST}}(Y, λ, ε) . V ariance
∇ 2 {df {ST}}(µ 0 , λ) = -1 √ 2πσ P i=0 exp - ((µ 0 ) i + λ) 2 2σ 2 + exp - ((µ 0 ) i -λ) 2 2σ 2 .
Thus, if µ 0 were given, the quantities in Proposition 2 could be computed in closed-form. The MSE can then be evaluated to select the optimal value of ε for a fixed dimension P and a given threshold λ. See the following numerical experiments which illustrate this relationship. When µ 0 is unknown, an a priori model can be imposed, such as for instance belonging to some ball promoting sparsity, e.g. a weak ℓ γ -ball for γ > 0. For γ sufficiently small, this ball corresponds to compressible or nearly sparse vectors µ 0 whose entries |µ i | sorted in descending order of magnitude behave as O(i -1/γ ). With such a model at hand, the MSE in Proposition 2 can be optimized for ε given P , σ, λ and γ. This however entails a highly non-linear equation that cannot be solved in closed form. We defer such a development to a future work.

Figure 1.(a) shows the evolution of the bias and the variance as a function of the ratio ε/σ for fixed values of σ, λ and a compressible vector µ 0 , i.e. |(µ 0 ) i | = O(i -1/γ ), chosen as illustrated on Figure 1.(d). When ε → 0, for fixed P , the bias vanishes while the variance, and in turn the MSE, increases. However, for a step ε > 0, the MSE is finite and seems to be optimal around the value 0.1σ. Figure 1.(c) shows the evolution of the MSE as a function of the dimension P and the ratio ε/σ for the same fixed values as before. The optimal step, minimizing the MSE, seems to evolve as a power decay function (the scale is log-log) of the form ε ⋆ (P ) = Cσ/P α with C > 0 and 0 < α < 1. Of course the optimal constants C and α depend on the choice of µ 0 , σ and λ. However, whatever C > 0 and 0 < α < 1, or more generally for any admissible choice of ε such that lim P →∞ ε(P ) = 0 and lim P →∞ P -1 ε(P ) -1 = 0, the MSE vanishes with respect to P . Figure 1.(b) shows indeed the evolution of the bias, the variance and the MSE as a function of the dimension P when ε is chosen as a power decay function. For α = 0 or α = 1, the MSE remains constant while, for α > 1, the MSE diverges which suggests the necessity of lim P →∞ P -1 ε(P ) -1 = 0.

Differentiation of an Iterative Scheme

We now turn to iterative algorithms that involve linear and soft-thresholding operators. We observed empirically that for all the inverse problems exposed in Section 5, setting ε ⋆ (P ) = Cσ/P α , as suggested by our study on the soft thresholding, resulted in a reliable way to parametrize our estimator. The effectiveness of this heuristic might be explained by the fact that the singularities encountered in most imaging problems are similar to absolute values, in order to encourage some sort of sparsity in the solution.

In this section, we focus on iterates, defined unambiguously as single-valued mappings (y, θ) → x (ℓ) (y, θ), where ℓ is the iteration counter of the iterative algorithm. In this context, we propose to compute in closed-form the derivatives of x (ℓ) (y, θ) with respect to either y (in a direction δ) or θ. This proves useful to respectively estimate the risk via SURE A MC (see Section 2) and estimate its gradient via SUGAR A FDMC (see Section 3). The iterative schemes we consider can be cast in the same framework, which subsumes proximal splitting algorithms designed to minimize a proper, closed and convex objective function x → E(x, y, θ), whose set of minimizers is supposed non-empty. All these algorithms can be unified as an iterative scheme of the form

x (ℓ) = γ(a (ℓ) ) a (ℓ+1) = ψ(a (ℓ) , y, θ), (16)
where a (ℓ) ∈ A is a sequence of auxiliary variables. ψ : A × Y × Θ → A is a fixed point operator in such a way that a (ℓ) converge to a fixed point a ⋆ , and γ : A → X is non-expansive (i.e, γ(a 1 )γ(a 2 ) a 1a 2 for any a 1 , a 2 ∈ A) entailing that x (ℓ) will converge to x ⋆ = γ(a ⋆ ). Note that for the sake of clarity, we have dropped the dependencies of a ⋆ and x ⋆ to y and θ.

To make our ideas clear, consider the instructive example where x → E(x, y, θ) is convex and C 1 (X ) with L-Lipschitz gradient, in which case A = X , a = x and ψ(x, y, θ) = xτ ∇ 1 E(x, y, θ) where 0 < τ < 2/L.

Iterative Weak Differentiability

A practical way to get the weak directional derivative ∂ 1 x(y, θ)[δ] and the weak Jacobian ∂ 2 x(y, θ), is to compute them iteratively from the sequences ( 16) by relying on the chain rule.

Algorithm Risk estimation of an iterative scheme

Inputs:

observation y ∈ Y = R P , collection of parameters θ ∈ Θ Parameters:

noise variance σ 2 > 0, linear operator Φ ∈ R P ×N , matrix A ∈ R M ×P , number L of iterations Output: solution x(y, θ) ∈ X and its risk estimate R A {x}(y, θ) However, two major issues have to be taken care of. First, one has to ensure weak differentiability of the iterates (16) so that ∂ 1 x (ℓ) (y, θ)[δ] (or resp. to ∂ 2 x (ℓ) (y, θ)) exist Lebesgue a.e. Second, one may legitimately ask whether the sequence of weak derivatives converges, and the properties of its cluster point, if any, with respect to the weak derivatives at a minimizer x ⋆ . Regarding weak differentiability of the iterates, it relies essentially on regularity conditions to apply the chain rule, e.g. [28, Section 4.2.2], i.e. regularity properties of the iteration mappings γ and ψ and of the initialization. For instance, for proximal splitting algorithms, it turns out that γ is the composition of one or several non-expansive operators, hence 1-Lipschitz, operators. In turn, γ is 1-Lipschitz. Furthermore, in all examples we consider, ψ is also 1-Lipschitz with respect to its second and third arguments. Therefore, if one starts at a Lipschitz continuous initialization, by induction, y → x (ℓ) (y, θ) and θ → x (ℓ) (y, θ) are also Lipschitz. Using the chain rule for Lipschitz mappings [28, Theorem 4 and Remark, Section 4.2.2], weak differentiability of x (ℓ) follows with respect to both arguments.

Sample a vector δ from N (0, Id P ) Initialize a (0) ← 0 * Initialize D (0) a ← 0 for ℓ from 0 to L -1 do * a (ℓ+1) ← ψ(a (ℓ) , y, θ) * D (ℓ+1) a ← Ψ (ℓ) a (D (ℓ) a ) + Ψ (ℓ) y (δ) end for * x (L) ← γ(a (L) ) * D (L) x ← Γ (L) a (D (L) a ) df A MC ← ΦD (L) x , A * Aδ SURE A MC ← A(y -Φx (L) ) 2 -σ 2 tr(A * A) + 2σ 2 df A MC return x(y, θ) ← x (L) and R A {x}(y, θ) ← SURE A MC
As far as convergence of the sequence of weak Jacobians is concerned, this remains an open question in the general case, and we believe this would necessitate intricate arguments from non-smooth and variational analysis. This is left to future research.

From now on, we suppose that the Lipschitzian assumptions on γ, ψ and the initial points hold. The next two sections detail the computation of ∂ 1 x (ℓ) (y, θ)[δ] and ∂ 2 x (ℓ) (y, θ) in order to get the estimates SURE A MC and SUGAR A FDMC .

Algorithm Risk and gradient risk estimation of an iterative scheme Inputs: observation y ∈ Y = R P , collection of parameters θ ∈ Θ Parameters:

noise variance σ 2 > 0, linear operator Φ ∈ R P ×N , matrix A ∈ R M ×P , number L of iterations decay parameters C > 0 and 0 < α < 1 Output: solution x(y, θ) ∈ X , its risk estimate R A {x}(y, θ), and its gradient risk estimate ∇2 R A {x}(y, θ)

Sample a vector δ from N (0, Id P ) * Choose ε = Cσ/P α * for y ′ = y and y

′ = y + εδ do * Initialize a (0) ← 0 * Initialize J (0) a ← 0 for ℓ from 0 to L -1 do * a (ℓ+1) ← ψ(a (ℓ) , y ′ , θ) * J (ℓ+1) a ← Ψ (ℓ) a (J (ℓ) a ) + Ψ (ℓ) θ end for * x (ℓ) (y ′ ) ← γ(a (ℓ) ) * J (ℓ) x (y ′ ) ← Γ (ℓ) a (J (ℓ) a ) end for * df FDMC ← 1 ε Φ(x (L) (y + εδ) -x (L) (y)), A * Aδ * SURE A FDMC ← A(y -Φx (L) ) 2 -σ 2 tr(A * A) + 2σ 2 df FDMC * SUGAR A FDMC ← 2J (L) x (y) * Φ * A * A(Φx (L) -y) + 2σ 2 ε J (L) x (y+εδ)-J (L) x (y) * Φ * A * Aδ return x(y, θ) ← x (L) (y), R A {x}(y, θ) ← SURE A FDMC and ∇2 R A {x}(y, θ) ← SUGAR A FDMC Figure 3
: Pseudo-algorithm for risk and gradient risk estimation of an iterative scheme. The symbols * indicate the lines corresponding to the computation of x and its estimated risk R A using approximated Monte Carlo simulation, i.e., as described in [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF]. The others are dedicated to the computation of the estimated gradient of the risk ∇R A . Even if computing the gradient of the risk requires more operations, the global complexity of the algorithm is unchanged.

Computation of SURE A MC for Risk Optimization

We describe here the iterative computation of the directional derivative ∂ 1 x (ℓ) (y, θ)[δ] following the idea introduced in [START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF] (see Section 2.4). Note that we focus on the directional derivative since, on the one hand, ∂ 1 x (ℓ) (y, θ) ∈ R N ×P is never used explicitly but only its trace, not to mention its storage cost, and, on the other hand, the risk can be estimated by applying only the weak directional derivatives on random directions δ (see Section 2.3 for more details).

The next proposition summarizes a recursive scheme to compute the weak derivatives

∂ 1 x (ℓ) (y, θ)[δ].
Proposition 3. For any vector δ ∈ X , the weak directional derivative D (ℓ)

x = ∂ 1 x (ℓ) (y, θ)[δ]
is given by

D (ℓ) x = Γ (ℓ) a (D (ℓ) a ) with D (ℓ+1) a = Ψ (ℓ) a (D (ℓ) a ) + Ψ (ℓ) y (δ),
where

D (ℓ) a = ∂ 1 a (ℓ) (y, θ)[δ]
and we have defined the following linear mappings [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF], and in turn in [START_REF] Cai | A data-driven block thresholding approach to wavelet estimation[END_REF], gives iteratively an unbiased § estimate of the risk at the current iterate x (ℓ) (y, θ). The whole procedure is summarized in Fig. 2. It is worth point out that although estimating the risk entails additional operations, the global complexity is the same as for the original iterative algorithm without risk estimation.

Γ (ℓ) a (•) = ∂ 1 γ(a (ℓ) )[•], Ψ (ℓ) a (•) = ∂ 1 ψ(a (ℓ) , y, θ)[•], Ψ (ℓ) y (•) = ∂ 2 ψ(a (ℓ) , y, θ)[•]. Plugging ∂ 2 x (ℓ) (y, θ)[δ] in

Computation of SUGAR A

FDMC for Risk Optimization

We now focus on the computation of the weak Jacobian ∂ 2 x (ℓ) (y, θ). Unlike for risk estimation that required only weak directional derivatives, for risk optimization we need the full weak Jacobian matrix ∂ 2 x(y, θ) ∈ R dim(Θ)×N . The proposed strategy, known as the forward accumulation, is one of the possible strategies to iteratively evaluate the derivatives by the use of the chain rule. The reverse accumulation is another strategy that does not require computing the full Jacobian matrix at the expense of a large memory load with respect to the number of iterations. Between these two extreme approaches, they are several hybrid strategies that can also be considered, knowing, that finding the optimal Jacobian accumulation strategy is an NP-complete problem. Such strategies have been studied in the field of "automatic differentiation" and the reader is invited to refer to [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Naumann | Optimal jacobian accumulation is NP-complete[END_REF] for a comprehensive account of these approaches.

In our case, we consider that, unlike for ∂ 1 x (ℓ) (y, θ), the matrix ∂ 2 x(y, θ) is in practice quite small since dim(Θ) ≪ P , hence implying only a memory load overhead of small fraction of P . Hence following the forward accumulation strategy, we propose a practical way to compute iteratively the full weak Jacobian matrix ∂ 2 x(y, θ).

The next result describes an iterative scheme to compute ∂ 2 x (ℓ) (y, θ).

Proposition 4. The weak Jacobian J

(ℓ) x = ∂ 2 x (ℓ) (y, θ) is given by J (ℓ) x = Γ (ℓ) a (J (ℓ) a ) with J (ℓ+1) a = Ψ (ℓ) a (J (ℓ) a ) + Ψ (ℓ) θ ,
where J

(ℓ) a = ∂ 2 a (ℓ) (y, θ) and we have defined

Γ (ℓ) a (•) = ∂ 1 γ(a (ℓ) )[•], Ψ (ℓ) a (•) = ∂ 1 ψ(a (ℓ) , y, θ)[•], Ψ (ℓ)
θ = ∂ 3 ψ(a (ℓ) , y, θ). § Expectation is to be taken here with respect to both the Gaussian measure of the noise W and the direction ∆.

Plugging ∂ 2 x (ℓ) (y, θ) in the expression of SUGAR A FDMC given by Proposition 1 provides iteratively an asymptotically (see Theorem 1) unbiased estimate of the gradient of the risk at the current iterate x (ℓ) (y, θ). The main steps of the procedure are summarized in Fig. 3. The estimation of the gradient of the risk entails only a small computational overhead compared to the risk estimation approach of [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF]. Their respective complexity remains however the same. Note finally that in both schemes, another initialization than a (0) = 0 can be chosen, for instance depending on y and θ, in which case the respective derivatives require to be initialized accordingly.

The following sections are devoted to instantiate this approach to more specific iterative algorithms that are able to handle non-smooth convex objective functions E.

Application to Generalized Forward Backward Splitting

The Generalized Forward Backward (GFB) splitting [START_REF] Raguet | A generalized forward-backward splitting[END_REF] allows one to find one element belonging to the set x(y, θ) solution of the structured convex optimization problem

x(y, θ) = Argmin x∈X E(x, y, θ) = F (x, y, θ) + Q k=1 G k (x, y, θ) (17)
under the assumptions that all functions are proper, closed and convex, F is C 1 (X ) with L-Lipschitz continuous gradient, and the G k functions are simple, in the sense their proximity operator can be computed in closed form (e.g. the ℓ 1 norm is simple since its proximal operator is explicitly the soft-thresholding). Recall that the proximal mapping of a proper closed convex function G is defined as

Prox G : x ∈ X → argmin z∈X 1 2 z -x 2 + G(z) .
It is uniquely valued and non-expansive (i.e., Prox G (x 1 ) -Prox G (x 2 )

x 1x 2 for any x 1 , x 2 ∈ X ), in fact even firmly so.

The GFB implements iteration [START_REF] Deledalle | Stein consistent risk estimator (SCORE) for hard thresholding[END_REF] with a (ℓ) = (ξ (ℓ) , z

(ℓ) 1 , . . . , z (ℓ) Q ) ∈ A = X 1+Q , x (ℓ) = γ(a (ℓ) ) = ξ (ℓ)
and a (ℓ+1) = ψ(a (ℓ) , y, θ) chosen such that for all k = 1, . . . , Q,

x (ℓ+1) = 1 Q Q k=1 z (ℓ+1) k and z (ℓ+1) k = z (ℓ) k -x (ℓ) + Prox νQG k (Z (ℓ) k , y, θ) with Z (ℓ) k = 2x (ℓ) -z (ℓ) k -ν∇ 1 F (x (ℓ) , y, θ) .
With the parameter ν ∈]0, 2/L[, the sequence of iterates x (ℓ) is provably guaranteed to converge to a minimizer x(y, θ) of [START_REF] Deledalle | Poisson NL means: Unsupervised non local means for Poisson noise[END_REF]. One recovers as a special cases the Forward-Backward splitting [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] when Q = 1 and the Douglas-Rachford splitting [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] when F = 0. Corollary 1. For any vector δ ∈ X , the GFB weak directional derivatives D

D (ℓ+1) x = 1 Q Q k=1 D (ℓ+1) z k and D (ℓ+1) z k = D (ℓ) z k -D (ℓ) x + G (ℓ) k,x (D (ℓ) Z k ) + G (ℓ) k,y (δ) with D (ℓ) Z k = 2D (ℓ) x -D (ℓ) z k -ν(F (ℓ) x (D (ℓ) x ) + F (ℓ) y (δ)) ,
where we have defined the following linear mappings

G (ℓ) k,x (•) = ∂ 1 {Prox νQG k }(Z (ℓ) k , y, θ)[•], G (ℓ) k,y (•) = ∂ 2 {Prox νQG k }(Z (ℓ) k , y, θ)[•], F (ℓ) x (•) = ∂ 1 {∇ 1 F }(x (ℓ) , y, θ)[•] and F (ℓ) y (•) = ∂ 2 {∇ 1 F }(x (ℓ) , y, θ)[•] .
Corollary 2. The GFB weak Jacobian J (ℓ)

x = Γ (ℓ) a (J (ℓ) a ), where J (ℓ+1) a = Ψ (ℓ) a (J (ℓ) a ) + Ψ (ℓ)
θ , is computed by evaluating iteratively

J (ℓ+1) x = 1 Q Q k=1 J (ℓ+1) z k and J (ℓ+1) z k = J (ℓ) z k -J (ℓ) x + G (ℓ) k,x (J (ℓ) Z k ) + G (ℓ) k,θ with J (ℓ) Z k = 2J (ℓ) x -J (ℓ) z k -ν(F (ℓ) x (J (ℓ) x ) + F (ℓ) θ ) ,
where we have defined

G (ℓ) k,x (•) = ∂ 1 {Prox νQG k }(Z (ℓ) k , y, θ)[•], G (ℓ) k,θ = ∂ 3 {Prox νQG k }(Z (ℓ) k , y, θ), F (ℓ) x (•) = ∂ 1 {∇ 1 F }(x (ℓ) , y, θ)[•] and F (ℓ) θ = ∂ 3 {∇ 1 F }(x (ℓ) , y, θ) .

Application to Primal-dual Splitting

Proximal splitting schemes can be used to find an element of the set x(y, θ) defined as the solution of the large class of variational problems

x(y, θ) = Argmin x∈X {E(x, y, θ) = H(x, y, θ) + G(K(x), y, θ)} , (18)
where both x → H(x, y, θ) and u → G(u, y, θ) are proper closed convex and simple functions, and K : X → U is a bounded linear operator. The primal-dual ¶ relaxed Arrow-Hurwicz algorithm as revitalized recently in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] (that we coin CP) to solve (18) implements ( 16) with a (ℓ) = (ξ (ℓ) , x(ℓ) , u (ℓ) ) ∈ A = X 2 × U , x (ℓ) = γ(a (ℓ) ) = ξ (ℓ) and a (ℓ+1) = ψ(a (ℓ) , y, θ) such that

u (ℓ+1) = Prox τ G * (U (ℓ) , y, θ) where U (ℓ) = u (ℓ) + τ K(x (ℓ) ), x (ℓ+1) = Prox ξH (X (ℓ) , y, θ) where X (ℓ) = x (ℓ) -ξK * (u (ℓ+1) ), x(ℓ+1) = x (ℓ+1) + ζ(x (ℓ+1) -x (ℓ) ). ( 19
)
where the Legendre-Fenchel conjugate of G is defined as G * (u, y, τ ) = max z z, u -G(z, y, τ ), and its proximity operator is given by Moreau's identity as

Prox τ G * (u, y) = u -τ Prox G/τ (u/τ, y) .
The parameters τ > 0, ξ > 0 are chosen such that τ ξ K 2 < 1, and ζ ∈ [0, 1] to ensure provable convergence of x (ℓ) toward an element in the set x(y, θ) of [START_REF] Deledalle | Proximal splitting derivatives for risk estimation[END_REF]. ζ = 0 corresponds to the Arrow-Hurwitz algorithm, and for ζ = 1, a sublinear O(1/ℓ) convergence rate on the partial duality gap was established in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. ¶ We invite the interested reader to consult [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF] for a detailed review on primal-dual algorithms.

Corollary 3. For any vector δ ∈ X , the CP weak directional derivatives D

D (ℓ+1) u = G (ℓ) u (D (ℓ) U ) + G (ℓ) y (δ) where D (ℓ) U = D (ℓ) u + τ K(D (ℓ) x ), D (ℓ+1) x = H (ℓ) x (D (ℓ) X ) + H (ℓ) y (δ) where D (ℓ) X = D (ℓ) x -ξK * (D (ℓ+1) u ), and
D (ℓ+1) x = D (ℓ+1) x + ζ(D (ℓ+1) x -D (ℓ) x )
where we have defined the following linear mappings

H (ℓ) x (•) = ∂ 1 {Prox ξH }(X (ℓ) , y, θ)[•], H (ℓ) y (•) = ∂ 2 {Prox ξH }(X (ℓ) , y, θ)[•], G (ℓ) u (•) = ∂ 1 {Prox τ G * }(U (ℓ) , y, θ)[•] and G (ℓ) y (•) = ∂ 2 {Prox τ G * }(U (ℓ) , y, θ)[•] .
Corollary 4. Similarly to Corollary 3, the CP weak Jacobians J (ℓ)

x = Γ (ℓ) a (J (ℓ) a ) and J (ℓ+1) a = Ψ (ℓ) a (J (ℓ) a ) + Ψ (ℓ)
θ are computed by evaluating iteratively

J (ℓ+1) u = G (ℓ) u (J (ℓ) U ) + G (ℓ) θ where J (ℓ) U = J (ℓ) u + τ K(J (ℓ) x ), J (ℓ+1) x = H (ℓ) x (J (ℓ) X ) + H (ℓ) θ where J (ℓ) X = J (ℓ) x -ξK * (J (ℓ+1) u
), and J

(ℓ+1) x = J (ℓ+1) x + ζ(J (ℓ+1) x -J (ℓ) x )
where we have defined

H (ℓ) x (•) = ∂ 1 {Prox ξH }(X (ℓ) , y, θ)[•], H (ℓ) θ = ∂ 3 {Prox ξH }(X (ℓ) , y, θ), G (ℓ) u (•) = ∂ 1 {Prox τ G * }(U (ℓ) , y, θ)[•] and G (ℓ) θ = ∂ 3 {Prox τ G * }(U (ℓ) , y, θ).
Note that the two proximal splitting schemes described here were chosen for their flexibility and the richness of the class of problems they can handle. Obviously, the methodology and discussion extend easily to the reader's favorite proximal splitting algorithm.

Examples and Numerical Results

In this section, we exemplify the use of the formal differentiation of iterative proximal splitting algorithms for three popular variational problems: nuclear norm regularization, totalvariation regularization and multi-scale wavelet ℓ 1 -analysis sparsity prior. For each of them, the expressions of all quantities including the proximal operators and their derivatives are given in closed-form. On each problem, we illustrate the usefulness of our gradient risk estimators for (multi) continuous parameter optimization.

Implementation Details

All experiments reported below are based on the algorithms detailed in Figure 2 and 3 in conjunction with proximal splitting algorithms presented in the previous section. The step of the finite difference is chosen as ε = 2σ/P 0.3 . Iterative proximal splitting algorithms will be used with L = 100 iterations. For quasi-Newton optimization, we used the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method with the implementation of [START_REF] Lewis | Nonsmooth optimization via quasi-newton methods[END_REF]. The use of BFGS is just given as an example to make the most of the proposed gradient estimator of the risk and seems to be good enough for the cases considered in all the following experiments examples. Of course, other first-order optimization methods can be considered for just as well, essentially if the risk presents several local minima.

An important issue in using quasi-Newton optimization is the choice of the initialization, the initial step and the stopping criteria. For a variation regularization problem expressed as Argmin

x 1 2 Φx -y 2 + K k=1 λ k R k (x) , (20)
where λ k > 0, ∀k ∈ N, the initialization λ k 0 is chosen empirically as

λ k 0 = P σ 2 4 K k=1 R k (x LS (y)) , (21)
where x LS (y) is the least-square estimator. At the first iteration, the approximate inverse Hessian B 1 should be chosen such that, for all k > 0, λ k 1 is of the same order as λ k 0 . To this end, we suggest initializing B 1 as a diagonal matrix with diagonal entries

B k 1 = αλ k 0 SUGAR A FDMC {x}(y, λ 0 , δ, ε) k (22)
such that, for all k, λ k 1 = (1 ± α)λ k 0 , where, in practice, we have chosen α = 0.9. Finally, the BFGS method stops after the following criterion is reached

SUGAR A FDMC {x}(y, λ n , δ, ε) ∞ SUGAR A FDMC {x}(y, λ 0 , δ, ε) ∞ τ ( 23
)
where we have chosen τ = 0.02 meaning that the algorithm stops if all (weak) partial derivatives are at least 50 times lower than the maximal one at initialization. For the sake of reproducibility, the Matlab scripts implementing the SURE and SUGAR for the different problems details hereafter are available online at http://www.math. u-bordeaux1.fr/ ~cdeledal/sugar.php.

Nuclear Norm Regularization

We consider the recovery of a low-rank matrix x 0 ∈ R n 1 ×n 2 from an observation y ∈ R P of Y = Φx 0 + W , W ∼ N (0, σ 2 Id P ), where we have identified the matrix space R n 1 ×n 2 to the vector space R N with N = n 1 n 2 . To this end, we consider the following spectral regularization problem

x ⋆ (y, λ) ∈ Argmin x 1 2 Φx -y 2 + λ x * , (24)
where λ > 0 and • * is the nuclear norm (a.k.a., trace for the symmetric semi-definite positive case or Schatten 1-norm). This is a spectral function defined as the ℓ 1 norm of the singular values Λ x ∈ R n=min(n 1 ,n 2 ) , i.e.

x * = Λ x 1 .

The nuclear norm is a particular case of spectral regularization that accounts for prior knowledge on the spectrum of x, typically low-rank (see, e.g., [START_REF] Fazel | Matrix Rank Minimization with Applications[END_REF]). It is the convex hull of the rank function restricted to the unit spectral ball [START_REF] Candès | Exact matrix completion via convex optimization[END_REF]. The parameter λ balances the sparsity of the spectrum of the recovered matrix, and the tolerated amount of noise. However, except in the random measurements setting, there is no direct relation between λ and the rank of x(y, λ). The optimal value of λ depends indeed on x 0 , Φ and σ confirming the importance of automatic selection procedures. Problem ( 24) is a special instance of ( 17) with the parameter λ = θ ∈ Θ = R + , Q = 1, and

F (x, y, λ) = 1 2 Φx -y 2
and G 1 (x, y, λ) = λ x * .

Hence the GFB algorithm can be used to solve [START_REF] Duval | A bias-variance approach for the non-local means[END_REF] by setting

∇ 1 F (x, y, λ) = Φ * (Φx -y) and Prox τ G 1 (x, y, λ) = V x diag(ST(Λ x , τ λ))U * x
where diag : R n → R n 1 ×n 2 maps the entries of a vector in R n to the main diagonal of a rectangular matrix in R n 1 ×n 2 filled with 0 elsewhere, (V

x , U x , Λ x ) ∈ R n 1 ×n 1 × R n 2 ×n 2 × R n is the singular value decomposition (SVD) of x such that x = V x diag(Λ x )U *
x and ST is the soft-thresholding operator [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]. Corollary 1 and 2 can then be applied using, for any δ x ∈ X and δ y ∈ Y, the relations

∂ 1 {∇ 1 F }(x, y, λ)[δ x ] = Φ * Φδ x , ∂ 2 {∇ 1 F }(x, y, λ)[δ y ] = -Φ * δ y , ∂ 3 {∇ 1 F }(x, y, λ) = 0 and ∂ 1 {Prox τ G 1 }(x, y, λ)[δ x ] = V x (H(Λ x )[ δx ] + Γ S (Λ x )[ δx ] + Γ A (Λ x )[ δx ])U * x , ∂ 2 {Prox τ G 1 }(x, y, λ)[δ y ] = 0, ∂ 3 {Prox τ G 1 }(x, y, λ) = V x diag(∂ 2 ST(Λ x , τ λ))U * x where δx = V * X δ x U X ∈ R n 1 ×n 2 , H(Λ x ) is defined as H(Λ x )[ δx ] = diag(∂ 1 ST(Λ x , ρλ)[diag( δx )])
and Γ S (Λ x ) and Γ A (Λ x ) are defined, for all 1 i n 1 and 1 j n 2 , as

Γ S (Λ x )[ δx ] i,j = ( δx ) i,j + ( δx ) j,i 2 ×      0 if i = j ST(Λx,ρλ) i -ST(Λx,ρλ) j (Λx) i -(Λx) j if (Λ x ) i = (Λ x ) j ∂ 1 ST(Λ x , ρλ) i,i otherwise, Γ A (Λ x )[ δx ] i,j = ( δx ) i,j -( δx ) j,i 2 ×      0 if i = j ST(Λx,ρλ) i +ST(Λx,ρλ) j (Λx) i +(Λx) j if (Λ x ) i > 0 or (Λ x ) j > 0 ∂ 1 ST(Λ x , ρλ) i,i otherwise,
where for i > n we have extended Λ x and ST(Λ x , ρλ) as (Λ x ) i = 0 and ST(Λ x , ρλ) i = 0, and for j > n 1 or i > n 2 , δx as ( δx ) j,i = 0. Recall from (32) that the weak derivatives of the soft-thresholding are defined, for t ∈ R N , ρ > 0, δ t ∈ R N , 1 i N , by

∂ 1 ST(t, ρ) i,i = 0 if |t i | ρ 1 otherwise, ( 25
)
∂ 1 ST(t, ρ)[δ t ] i = ∂ 1 ST(t, ρ) i,i × (δ t ) i and ∂ 2 ST(t, ρ) i = 0 if |t i | ρ -sign(t i ) otherwise.
which corresponds in this case where Q = 1 to the forward-backward algorithm.

The closed-form expression we derived for ∂ 1 {Prox τ G 1 }(x, y, λ)[δ x ] is far from trivial. It is essentially due to [START_REF] Edelman | Matrix Jacobians with wedge products[END_REF][START_REF] Lewis | Twice differentiable spectral functions[END_REF][START_REF] Sun | Nonsmooth matrix valued functions defined by singular values[END_REF], see [START_REF] Candès | Unbiased risk estimates for singular value thresholding and spectral estimators[END_REF] for an expression similar to ours. The generalization of this result to other matrix-valued spectral function has been studied in [START_REF]Risk estimation for matrix recovery with spectral regularization[END_REF].

Application to matrix completion We illustrate the nuclear norm regularization on a matrix completion problem encountered in recommendation systems such as the popular Netflix problem [START_REF] Bennett | The Netflix prize[END_REF]. We therefore consider y ∈ R P with the forward model Y = Φx 0 + W , W ∼ N (0, σ 2 Id P ), where x 0 is a dense but low-rank (or approximately so) matrix and Φ is a binary masking operator.

We have taken (n 1 , n 2 ) = (1000, 100) and P = 25000 observed entries (i.e., 25%). The underlying matrix x 0 = V x 0 diag Λ x 0 U * x 0 has been chosen with V x 0 and U x 0 two realizations of the uniform distribution of orthogonal matrices and Λ x 0 = (k -1 ) 1 k n such that x 0 is approximately low-rank with a rapidly decaying spectrum. The binary masking operator is such that for i = 1, . . . , P , (Φx

) i = x Σ(i) 1 ,Σ(i) 2 where Σ : [1, . . . , n 1 × n 2 ] → [1, . . . , n 1 ] × [1, . . . , n 2 ]
is the realization of a random permutation of the n 1 × n 2 entries of x. The standard deviation σ has been set such that the resulting minimum least-square estimate x LS (y) = Φ * y has a relative error x LS (y)x 0 F / x 0 F = 0.9.

Figure 4.(a) and (b) depict the risk, the prediction risk and the SURE = SURE A (with A = Id) estimates * * as a function of λ obtained from a single realization of y and δ. In (a), SURE MC {x}(y, λ, δ) has been evaluated for 12 values of λ chosen in a suitable tested range using the algorithm given in Figure 2. Figure (b), shows the benefit of computing SURE FDMC {x}(y, λ, δ, ε) and SUGAR FDMC {x}(y, λ, δ, ε), as described in Figure 3, to realize a quasi-Newton optimization. The sequence of iterates λ n is represented as well as the sequence of the slopes of SURE FDMC {x}(y, λ n , δ, ε) given by SUGAR FDMC {x}(y, λ n , δ, ε). The BFGS algorithm reaches to the optimal value in 5 iterations only. One can also notice that SURE FDMC {x}(y, λ, δ, ε) and SURE MC {x}(y, λ, δ) are both good, and visually equivalent, estimators of the prediction risk. At the optimum value λ ⋆ minimizing the SURE, the true risk is not too far from its minimum showing that, in this case, the prediction risk is indeed a good objective in order to minimize the risk. In Figure 4.(c) a close in on the solution x(y, λ ⋆ ) is compared to x 0 and x LS (y), and their respective spectrum in Figure 4.(d). The solution x(y, λ ⋆ ) has a rank of 57 with a relative error of 0.45 (i.e., a gain of about a factor 2 w.r.t. the least-square estimator).

Total-Variation Regularization

We consider the recovery of a piece-wise constant two dimensional image x 0 ∈ R n 1 ×n 2 from an observation y of Y = Φx 0 + W ∈ R P , W ∼ N (0, σ 2 Id P ), where we have identified the image space R n 1 ×n 2 to the vector space R N with N = n 1 n 2 . To this end, we suggest using (isotropic) total variation regularization of the form

x ⋆ (y, λ) ∈ Argmin x 1 2 Φx -y 2 + λ ∇x 1,2 , (26)
where λ > 0 and ∇ : R N → R N ×2 is the two-dimensional discrete gradient operator. The ℓ 1 -ℓ 2 norm of a vector field t = (t i ) N i=1 ∈ R N ×2 , with t i ∈ R 2 , is defined as t 1,2 = i t i . Total-variation promotes the sparsity of the gradient field which turns out to be a prior that enforces smoothing while preserving edges. The parameter λ controls the regularity of the image. A large value of λ results to an image with large homogeneous areas while a small value results to an image with several small disconnected regions. The optimal value of λ is image and degradation dependent revealing the importance of automatic selection procedures. Problem ( 26) is a special instance of (17

) using x = (f, u) ∈ X = R N × R N ×2 , the parameter λ = θ ∈ Θ = R + , Q = 2 simple functionals, and for x = (f, u) F (x, y, λ) = 1 2 Φf -y 2 , G 1 (x, y, λ) = λ u 1,2 and G 2 (x, y, λ) = ι C (x) where C = x = (f, u) \ u = ∇f .
Hence the GFB algorithm can be used to solve (26) using

∇ 1 F (x, y, λ) = (Φ * (Φf -y), 0), Prox τ G 1 (x, y, λ) = (f, ST 1,2 (u, τ λ)) and Prox τ G 2 (x, y, λ) = ((Id + ∆) -1 (f + div u), ∇(Id + ∆) -1 (f + div u))
where ∆ is the Laplacian operator and div is the discrete divergence operator such that div = -∇ * . The operator ST 1,2 is the component-wise ℓ 1 -ℓ 2 soft-thresholding defined, for any dimensions N and D, t ∈ R N ×D and ρ > 0, by

ST 1,2 (t, ρ) i = 0 if t i ρ t i -ρ t i / t i otherwise , for all 1 i N. (27)
For D = 1, the component-wise ℓ 1 -ℓ 2 soft-thresholding reduces to [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]. Corollary 1 and 2 can then be applied, for any δ x = (δ f , δ u ) ∈ X and δ y ∈ Y, using the relations

∂ 1 {∇ 1 F }(x, y, λ)[δ x ] = (Φ * Φδ f , 0), ∂ 2 {∇ 1 F }(x, y, λ)[δ y ] = (-Φ * δ y , 0), ∂ 3 {∇ 1 F }(x, y, λ) = (0, 0), ∂ 1 {Prox τ G 1 }(x, y, λ)[δ x ] = (δ f , ∂ 1 ST 1,2 (u, τ λ)[δ u ]), ∂ 2 {Prox τ G 1 }(x, y, λ)[δ y ] = (0, 0), ∂ 3 {Prox τ G 1 }(x, y, λ) = (0, ∂ 2 ST 1,2 (u, τ λ))
and

∂ 1 {Prox τ G 2 }(x, y, λ)[δ x ] = ((Id + ∆) -1 (δ + div δ u ), ∇(Id + ∆) -1 (δ f + div δ u )), ∂ 2 {Prox τ G 2 }(x, y, λ)[δ y ] = (0, 0), ∂ 3 {Prox τ G 2 }(x, y, λ) = (0, 0)
where the weak derivatives of the component-wise ℓ 1 -ℓ 2 soft-thresholding are defined, for any dimensions N and D, t ∈ R N ×D , ρ > 0 and δ t ∈ R N ×D , by

∂ 1 ST 1,2 (t, ρ)[δ t ] i = 0 if t i ρ δ t,i -ρ t i P t i (δ t,i ) otherwise (28)
and

∂ 2 ST 1,2 (t, ρ) i = 0 if t i ρ -t i / t i otherwise
where P α is the orthogonal projector on α ⊥ for α ∈ R 2 .

Application to image deblurring We illustrate the total-variation regularization on an image deblurring problem. We therefore consider the forward model Y = Φx 0 + W ∈ R P , W ∼ N (0, σ 2 Id P ), where x 0 is a piece-wise constant (or approximately so) image and Φ is a discrete convolution matrix.

We have taken a cartoon-like image of size (n 1 , n 2 ) = (512, 512) and P = 512 2 observations corresponding to noisy observations of a convolution product with a discrete Gaussian kernel of radius 2 pixels. To ensure numerical stability of the pseudo-inverse (typically for the leastsquare estimate and the computation of the projection risk and its estimate), the kernel has been truncated in the Fourier domain such that too small contributions have been set to 0. The consequence is that around 80% of (high) frequencies are masked. The standard deviation of the noise has been set to σ = 10 (for an image x 0 with a range [0, 255]) such that the resulting minimum least-square estimate x LS (y) = Φ + y has a peak signal-to-noise ratio (PSNR) equals to 10 log 10 (255 2 / x LS (y)x 0 2 F ) = 21.02 dB. Figure 5.(a) and (b) display the risk, the projection risk and the GSURE = SURE A (with A = Π) estimates as a function of λ obtained from a single realization of y and δ. In (a), GSURE MC {x}(y, λ, δ) has been evaluated for 12 values of λ chosen in a suitable tested range using the algorithm given in Figure 2. Figure (b), shows the benefit of computing GSURE FDMC {x}(y, λ, δ, ε) and GSUGAR FDMC {x}(y, λ, δ, ε), as described in Figure 3, to realize a quasi-Newton optimization. The sequence of iterates λ n is represented as well as the sequence of the slopes of GSURE FDMC {x}(y, λ n , δ, ε) given by GSUGAR FDMC {x}(y, λ n , δ, ε).

The BFGS algorithm reaches to the optimal value in 4 iterations. The deviation of GSURE FDMC {x}(y, λ, δ, ε) from the projection risk is of the same order as the deviation of GSURE MC {x}(y, λ, δ). At the optimum value λ ⋆ minimizing the GSURE, the true risk is not too far from its minimum showing that relatively to the range of variation of the risk, in this case, the projection risk is indeed a good objective in order to minimize the risk. In Figure 5.(c-e) the solution x(y, λ ⋆ ) is compared to x 0 and y. The solution x(y, λ ⋆ ) has a PSNR of 24.98 dB (i.e., a gain of about +3.94 dB). Remark that given such noise level and convolution operator, masking 80% of (high) frequencies, the solution selected by minimizing the GSURE criterion is still a bit blurred and exhibit a staircasing effect. Using a larger regularization parameter λ would result in a more "cartoon" result with less blur. This would however entail a larger bias, which corresponds to a loss of contrast inherent to the convexity of the TV prior. This larger bias subsequently degrades the mean square error, which explains why it is not selected by the GSURE criterion.

Table 1: Illustration of the minimization of SURE FDMC in multi-scale regularization obtained for the three images of Figure 6 with different numbers of scales from J = 1 to J = 3 using either one global parameter or one parameter per scales. For each case, the obtained optimal parameters λ ⋆ are given. The associated value of SURE and the PSNR are compared to neighbors of λ ⋆ located at 0.75λ ⋆ and 1.25λ

Weighted ℓ 1 -analysis Wavelet Regularization

We focus on the recovery of a piece-wise regular image x 0 ∈ R n 1 ×n 2 from an observation y of Y = Φx 0 + W ∈ R P , W ∼ N (0, σ 2 Id P ), using a J-scale undecimated wavelet analysis regularization of the form

x ⋆ (y, λ) ∈ Argmin x 1 2 Φx -y 2 + Ψx 1,λ where Ψ =        Ψ h 1 Ψ v 1 . . . Ψ h J Ψ v J        (29)
and λ ∈ R + J and Ψ ∈ R 2JN ×N is the analysis operator of a two-orientation wavelet transform, where, for all scales 1 j J, Ψ h j , Ψ v j are defined such that, for x ∈ R N , u h j = Ψ h j x and u v j = Ψ v j x are respectively the vectors of undecimated wavelet coefficients of x in the horizontal and vertical directions at the decomposition level j. The weighted ℓ

1 -norm . 1,λ is Ψx 1,λ = J j=1 λ j Ψ h j x 1 + Ψ v j x 1 .
Multi-scale wavelet analysis promotes piece-wise regular images by enforcing smoothness while preserving sharp discontinuities at different scales and orientations. Each parameter λ j controls the regularity at scale j. A large value of λ j tends to over-smooth structures at scale j, while a small value leads to under-smoothing. As noted in several papers, see e.g. [START_REF] Chaux | A nonlinear stein-based estimator for multichannel image denoising[END_REF][START_REF] Pesquet | A new wavelet estimator for image denoising[END_REF], the optimal values λ j are also image and degradation dependent revealing again the importance of automatic selection procedures.

Problem ( 29) is a special instance of [START_REF] Deledalle | Proximal splitting derivatives for risk estimation[END_REF] where the parameter λ = θ ∈ Θ = R + J , and

H(x, y, λ) = 1 2 Φx -y 2 , G(u, y, λ) = u 1,λ and K(x) = Ψx .
Hence the primal-dual CP splitting can be used to solve [START_REF] Fazel | Matrix Rank Minimization with Applications[END_REF] using

Prox ξH (x, y, λ) = x + ξΦ * y -ξΦ * (Id + ξΦΦ * ) -1 Φ(x + ξΦ * y), Prox τ G * (u, y, λ) = u -τ ST(u/τ, λ/τ ) and K * (u, λ) = J j=1 (Ψ h j * u h j + Ψ v j * u v j )
where ST denotes here the multi-scale extension of the soft-thresholding operator [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF], such that, for t ∈ R 2JN and ρ ∈ R J , we have

ST(t, ρ) o j = ST(t o j , ρ j )
for all scales 1 j J and orientations o = v, h. Corollary 1 and 2 can then be applied using

∂ 1 {Prox ξH }(x, y, λ)[δ x ] = δ x + ξΦ * (Id + ξΦΦ * ) -1 Φδ x , ∂ 2 {Prox ξH }(x, y, λ)[δ y ] = ξΦ * δ y -ξ 2 Φ * (Id + ξΦΦ * ) -1 ΦΦ * δ y , ∂ 3 {Prox ξH }(x, y, λ) = 0 and ∂ 1 {Prox τ G * }(u, y, λ)[δ u ] = δ u -∂ 1 ST(u/τ, λ/τ )[δ u ], ∂ 2 {Prox τ G * }(u, y, λ)[δ y ] = 0, ∂ 3 {Prox τ G * }(u, y, λ) = -∂ 2 ST(u/τ, λ/τ ).
where the derivatives of the multi-scale soft-thresholding are defined, for any t ∈ R 2JN , ρ ∈ R J and δ t ∈ R 2JN , by

∂ 1 ST(t, ρ)[δ t ] o j = ∂ 1 ST(t o j , ρ j )[δ t o j ] and ∂ 2 ST(t, ρ) o j = ∂ 2 ST(t o j , ρ j ) (30)
for all scales 1 j J and orientations o = v, h.

Application to compressed sensing We illustrate the multi-scale wavelet ℓ 1 -analysis regularization on a compressed sensing problem. We therefore consider the forward model Y = Φx 0 + W ∈ R P , W ∼ N (0, σ 2 Id P ), where x 0 is a piece-wise multi-scale regular (or approximately so) image and Φ is a random matrix. Here the multi-scale transform W is constructed from undecimated Daubechies 4 wavelets [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF].

We have taken a uniformly randomized sub-sampling of a uniform random convolution, where (P/N = 0.5). The standard deviation has been set to σ = 10 (for an image x 0 with a range [0, 255]) such that the resulting minimum least-square estimate x LS (y) = Φ + y has a PSNR given by 10 log 10 (255 2 / x LS (y)x 0 2 F ) ≈ 16 dB. Table 1 and Figure 6 illustrates the multi-scale regularization obtained by minimizing the SURE = SURE A (with A = Id) for three different images x 0 , known as Mandrill, House and Cameraman, and a single realization of y and δ. Three levels of decomposition from J = 1 to 3 are considered. We consider also to use either one global regularization parameter or one parameter per scales. Table 1 gives the selected optimal vector of parameters λ ⋆ for each level of decomposition and their associated performance in terms of SURE and PSNR. We first observe that compared to the global approach, optimizing one parameter per scale indeed adapts better to the regularity of the image. For instance, the image Mandrill contains fine scales with more energy than House, and then the obtained penalization of the first scale is smaller for Mandrill than for House. Visual inspection of these results on Figure 6 illustrates this automatic adaptation. In the same vein, with three levels of decomposition, the penalization is less severe for Mandrill than for House and Cameraman. We next observe that increasing the level of decomposition improves the PSNR when using one parameter per scale, while this is not the case when a global parameter is used. The gap is more important between J = 1 and J = 2. To assess the minimization of SURE FDMC , we have compared the SURE and the PSNR values at 0.75λ ⋆ and 1.25λ ⋆ . At the optimal λ ⋆ , the SURE is as expected minimal. Furthermore, at λ ⋆ , the PSNR is either maximal or not too far from its maximal value, showing that, in this case, the prediction risk is indeed a good objective in order to maximize the PSNR.

Conclusion

We have proposed a methodology for optimizing multiple continuous parameters of a weakly differentiable estimator that attempts solving a linear ill-posed inverse problem contaminated by additive white Gaussian noise. The proposed method selects the parameters minimizing an estimate of the risk and is driven by an estimate of its gradient. Classical unbiased estimators of the risk are generally non-continuous functions of the parameters, so that, their local variations cannot be used to estimate the gradient of the risk. These estimators require estimating the degree of freedom by evaluating the variations of the estimator with respect to the observations. We have shown that estimating the degree of freedom by finite differences leads to a weakly differentiable risk estimator. By carefully choosing the finite differences step and by computing explicitly the (weak) gradient of this estimate, an asymptotically unbiased estimator of the gradient of the risk is obtained. This estimator is numerically smooth enough to apply a quasi-Newton method. An explicit strategy to compute this (weak) gradient is given for a large class of (iterative) weakly differentiable algorithms. We exemplified our methodology on several popular proximal splitting methods. Numerical experiments have demonstrated the wide applicability and scope of the approach.

Our choice of the finite differences step size was essentially guided by a careful analysis of the soft-thresholding estimator. Choosing this step size with theoretical guarantees (such as consistency or optimality) in more general cases remains an open question. Beyond consistency and optimality, the question of quantifying the influence of the finite differences step on the smoothness of the risk gradient estimates and then on the performance of quasi-Newton methods is still open. To deal with parameter space of higher dimensions, other accumulation Jacobian strategies could be explored following [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]. Improvements could also be achieved on the settings of the quasi-Newton methods. In particular, a drawback of our approach is the sensitivity to local minima of the risk with respect to the collection of parameters. In some settings, more elaborated optimization strategies could be employed. Future work could also focus on the extensions to non-weakly differentiable estimators and/or inverse problems with non-Gaussian noises.

A Proofs of Section 3

Proposition 1. This is a consequence of the chain rule and linearity of the weak derivative. Indeed, df A FD {µ}(y, θ, ε) is just the sum of P weakly differentiable functions, and hence is weakly differentiable with the weak derivative with respect to θ as given. Moreover, A(µ(y, θ)y) 2 = P i=1 ((A(µ(y, θ)y)) i ) 2 . Each term i is the composition of a weakly differentiable function (A(µ(y, •)y)) i and (•) 2 , where the latter is obviously continuously differentiable with bounded derivative, and takes 0 at the origin. It then follows from the chain rule [28, Theorem 4(ii), Section 4.2.2] that (A(µ(y, •)y)) i is weakly differentiable, and the weak derivative of A(µ(y, •)y) 2 with respect to θ is indeed 2∂ 2 µ(y, θ) * A * A(µ(y, θ)y) .

Theorem 1. The proof strategy consists in commuting in an appropriate order the different signs (limit, integration and derivation) while checking that our assumptions provide sufficient conditions for this to hold.

Let V a compact subset of Θ, and choose ϕ ∈ C 1 c (Θ) with support in V . We have

Θ R A {µ}(µ 0 , θ) ∂ϕ(θ) ∂θ i dθ = V R A {µ}(µ 0 , θ) ∂ϕ(θ) ∂θ i dθ = V E W A(µ(Y, θ) -y) 2 ∂ϕ(θ) ∂θ i dθ [Stein Lemma] (S.1) = V E W SURE A {µ}(Y, θ) ∂ϕ(θ) ∂θ i dθ [(12)] (S.2) = V E W lim ε→0 SURE A FD {µ}(Y, θ, ε) ∂ϕ(θ) ∂θ i dθ [Dominated convergence] (S.3) = lim ε→0 V E W SURE A FD {µ}(Y, θ, ε) ∂ϕ(θ) ∂θ i dθ [Fubini] (S.4) = lim ε→0 E W V SURE A FD {µ}(Y, θ, ε) ∂ϕ(θ) ∂θ i dθ [Weak differentiability, Proposition 1] (S.5) = -lim ε→0 E W V ∂ ∂θ i SURE A FD {µ}(Y, θ, ε)ϕ(θ)dθ [Proposition 1] (S.6) = -lim ε→0 E W V SUGAR A FD {µ}(Y, θ, ε) i ϕ(θ)dθ [Fubini] (S.7) = -lim ε→0 V E W SUGAR A FD {µ}(Y, θ, ε) i ϕ(θ)dθ [Dominated convergence] (S.8) = - V lim ε→0 E W SUGAR A FD {µ}(Y, θ, ε) i ϕ(θ)dθ = - Θ lim ε→0 E W SUGAR A FD {µ}(Y, θ, ε) i ϕ(θ)dθ .
From the definition of weak derivative, we get the claimed result on the asymptotic unbiasedness of SUGAR A FD . The asymptotic unbiasedness of the gradient of the finite difference DOF naturally follows with the same proof strategy by ignoring the two first terms in the decomposition SURE A FD {µ}(µ 0 , θ, ε) = A(µ(y, θ)y) 2σ 2 tr(A * A) + 2σ 2 df A FD {µ}(µ 0 , θ, ε). We now justify each of the steps (S.1)-(S.8). We denote g 1,σ the Gaussian probability density function of zero-mean and variance σ 2 , and g σ its P -dimensional version, i.e. g σ = (g 1,σ ) P .

(S.1) This is Stein lemma which applies owing to Assumption (A.1). Indeed, µ(•, θ) is Lipschitz, hence weakly differentiable and its derivative equals its weak derivative Lebesgue a.e. [28, Theorem 1-2, Section 6.2]. Moreover, we have for any θ

µ(y, θ) -µ(y ′ , θ) L 1 y -y ′ ⇒ µ i (y, θ) -µ i (y ′ , θ) L 1 y -y ′ , (31)
and thus, whenever the derivatives of µ i (•, θ) exist, they are bounded by L 1 . Consequently,

E W ∂µ i (Y ) ∂y i L 1 ,
i.e. the weak partial derivatives are essentially bounded.

(S.2) This is just [START_REF]Fixed-Point Algorithms for Inverse Problems in Science and Engineering[END_REF] and the arguments justifying hold owing to Assumption (A.1).

(S.3) Let f ε (y, θ) = SURE A FD {µ}(y, θ, ε). From ( 12), lim ε→0 f ε (y, θ) = SURE A {µ}(y, θ) exists Lebesgue a.e. Assumptions (A.1)-(A.2) give µ(y, θ)y y + µ(y, θ)µ(0, θ) (1 + L 1 ) y .

Combining this with (31) leads to

|f ε (y, θ)| = A(µ(y, θ) -y) 2 -σ 2 tr(A * A) + 2σ 2 1 ε P i=1 (A * A(µ(y + εe i , θ) -µ(y, θ))) i A 2 P σ 2 (1 + L 1 ) 2 y 2 P σ 2 + 1 + 2L 1 .
Note that the bound is independent of θ. Thus

E W A 2 P σ 2 (1 + L 1 ) 2 y 2 P σ 2 + 1 + 2L 1 = A 2 P σ 2 (1 + L 1 ) 2 µ 0 2 P σ 2 + 1 + 1 + 2L 1 < ∞ .
This bound together with the fact that ϕ is continuously differentiable with compact support in V means that f ε ∂ϕ ∂θ i is dominated by an integrable function on Y × V . The dominated convergence then applies which yields the claim.

(S.4) Fubini theorem surely applies in view of the integrability just shown at the end of (S.3). (S.7) Let f ε (y, θ) = (SUGAR A FD {µ}(y, θ, ε)) i and h(y, θ) = (2∂ 2 µ(y, θ) * A * A(µ(y, θ)y)) i . By the translation invariance of the convolution product, we have

E W [f ε (Y, θ)] = 2(g σ * h(•, θ))(µ 0 ) + 2σ 2 P j=1 g σ (• + εe i ) -g σ ε * (∂ 2 µ(•, θ) * A * Ae j ) i (µ 0 ) . Thus |(g σ * h(•, θ))(µ 0 )| 2 Y g σ (y -µ 0 ) |(∂ 2 µ(y, θ) * A * A(µ(y, θ) -y)) i | dy [Assumption (A.3)] 2L 2 A 2 Y g σ (y -µ 0 ) µ(y, θ) -y dy [Assumptions (A.1)-(A.2)] 2(1 + L 1 )L 2 A 2 Y g σ (y -µ 0 ) y dy 2(1 + L 1 )L 2 A 2 E W [ y ]dy [Jensen inequality] 2(1 + L 1 )L 2 A 2 E W [ y 2 ] 1/2 dy 2(1 + L 1 )L 2 A 2 µ 0 2 + P σ 2 1/2 < ∞ .
For the second term, we have

g σ (• + εe i ) -g σ ε * (∂ 2 µ(•, θ) * A * Ae j ) i (µ 0 ) Y g σ (y -µ 0 + εe i ) -g σ (y -µ 0 ) ε (∂ 2 µ(y, θ) * A * Ae j ) i dy [Assumption (A.3)] L 2 A 2 Y g σ (y -µ 0 + εe i ) -g σ (y -µ 0 ) ε dy L 2 A 2 R g 1,σ (t -(µ 0 ) i + ε) -g 1,σ (t -(µ 0 ) i ) ε dt [Taylor] L 2 A 2 R 1 0 g ′ 1,σ (t -(µ 0 ) i + τ ) dtdτ [Fubini] L 2 A 2 R g ′ 1,σ (t) dt < ∞ .
In view of these bounds, and since ϕ is compactly supported in V , integrability of f ε ϕ on Y × V is ensured, whence the claimed result follows.

(S.8) Let f ε defined as in (S.7).

We have just shown that the integrand in θ, i.e. E W [f ε (Y, •)) i ]ϕ, is dominated by a function that is integrable on V . It remains to check that its limit exists Lebesgue a.e. But this is yet again an application of the dominated convergence theorem to the sequence f ε as an integrand with respect to the Gaussian measure g σ (y)dy, which allows to deduce that lim

ε→0 E W [f ε (Y, θ)ϕ(θ)] = E W [lim ε→0 f ε (Y, θ)ϕ(θ)].
This completes the proof.

Proposition 1. For a fixed λ, it can be shown similarly to [START_REF] Evans | Measure theory and fine properties of functions[END_REF]Theorem 4(iii), Section 4.2.2], that ST(•, λ) is weakly differentiable and that its weak Jacobian h(y) = ∂ 2 ST(y, λ) is diagonal, with diagonal elements, for 1 i P ,

h(y) i =    +1 if y i -λ 0 if -λ < y i < λ -1 otherwise . ( 32
)
We next define, for a fixed λ, the quantity h ′ (y, ε) = ∇ 2 { df FD {ST}}(y, λ, ε). Using Proposition 1 and the fact that ε < 2λ give

h ′ (y, ε) = P i=1 h(y + εe i ) i -h(y) i ε = P i=1            0 if y i < -λ -ε -1/ε if -λ -ε < y i < -λ 0 if -λ < y i < λ -ε -1/ε if λ -ε < y i < λ 0 if λ < y i .
Computing the expectation and the variance of h ′ (Y, ε) in closed-form with truncated Gaussian statistics, and using the fact that h is separable in its arguments, give the proposed formula.

Theorem 2. For P big enough, ε(P ) < 2λ since lim P →∞ ε(P ) = 0. Using the notations in the proof of Lemma 1 leads to SUGAR FD {ST}(y, λ, ε) = 2h(y) * (ST(y, λ)y) + 2σ 2 h ′ (y, ε(P )) .

The Cauchy-Schwartz inequality implies that By assumption, we have lim P →∞ P -1 ε(P ) -1 = 0 and then the variance of 1 P h ′ (Y, ε(P )) vanishes to zero. Next, remark that h(y) * (ST(y, λ)y) = λ#{|y| > λ} where #{|y| > λ} denotes the number of entries of |y| greater than λ. We have #{|Y i | > λ} ∼ iid Bernoulli(p i ) whose variance is p i (1-p i ), where p i = 1 2 erf (µ 0 ) i +λ

V

B Regularity of the proximal operator of a gauge

We first provide a glimpse of gauges. We now summarize some key properties that will be needed in the main proof. We are now equiped to prove our regularity result.

where the first equality follows from positive homogeneity (Lemma 2(i)) and Definition 2, and the second equality is a consequence of Lemma 2(iv) and polarity.

Applying Moreau identity, we get that y -Prox θG (y) = y -Prox σ θC • (y) = Proj θC • (y) .

By virtue of Lemma 2(iii), there exists a constant L 2 > 0, independent of y, such that

where the last inequality follows from Lemma 2(ii) since obviously Proj θC • (y) ∈ θC • . Combining [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] and [START_REF] Attouch | A continuous dynamical Newton-like approach to solving monotone inclusions[END_REF], we get the desired result. † † The constant L2 can be given explicitly by bounding from below the support function of the inscribed ellipsoid of maximal volume, the so-called John ellipsoid. For symmetric convex bodies, L2 can be made tightest possible. For simplicity, we avoid delving into these technicalities here. where the required Lipschitz continuity assumptions w.r.t. θ hold true.

Corollary 1. We first notice that D

z k = ∂ 2 z (ℓ)
k (y, θ). Arguing as in the proof of Corollary 1, using now that J

  Bias-variance trade-off w.r.t. ε Bias-variance w.r.t. P (using ε(P )) MSE w.r.t. (ε, P ) Choice of σ, µ0 and λ

Figure 1 :

 1 Figure 1: Bias-variance trade-off of the gradient estimator of the DOF of soft-thresholding, (a) with respect to the step ε and (b) with respect to the dimension P when using a power decay function P → ε(P ). (c) Its mean squared error as a function of P and ε (in logarithmic scales). The solid line represents the pairs (ε ⋆ (P ), P ) where, for a fixed dimension P , ε⋆ (P ) minimizes the mean squared error. The function ε⋆ (P ) looks like a power function of the form of Cσ/P α with C > 0 and 0 < α < 1. The dashed lines represent respectively the power functions εinf (P ) = Cσ and εsup (P ) = Cσ/P outside which the mean squared error diverges when P increases. (d) Description of the settings of the experiments, i.e., the choice of σ, µ 0 and λ.

Figure 2 :

 2 Figure 2: Pseudo-algorithm for risk estimation of an iterative scheme. The symbols * indicate the lines corresponding to the computation of x. The others are dedicated to the computation of the estimated risk R A using Monte Carlo simulation. Even if computing the risk requires more operations, the global complexity of the algorithm is unchanged.

y

  (δ) are computed by evaluating iteratively

y

  (δ) are computed by evaluating iteratively

Figure 4 :

 4 Figure 4: (a-b) Risk, prediction risk and its SURE estimates * * as a function of the regularization parameter λ. (a) The 12 points where SURE MC {x}(y, λ, δ) has been evaluated by exhaustive search. (b) The 5 evaluation points of SURE FDMC {x}(y, λ, δ, ε) and SUGAR FDMC {x}(y, λ, δ, ε) required by BFGS to reach the optimal one. (c-d) Respectively, a close in and the spectrum of the underlying matrix x 0 , the least-square estimate x LS (y) and the solution x(y, λ) at the optimal λ.

Figure 5 :

 5 Figure 5: (a-b) Risk, projection risk and its GSURE estimates * * as a function of the regularization parameter λ. (a) The 12 points where GSURE MC {x}(y, λ, δ) has been evaluated by exhaustive search. (b) The 4 evaluation points of GSURE FDMC {x}(y, λ, δ, ε) and GSUGAR FDMC {x}(y, λ, δ, ε) required by BFGS to reach the optimal one. (c-d) Respectively, a close in of the underlying image x 0 , the observation y and the solution x(y, λ) at the optimal λ.

Figure 6 :

 6 Figure 6: From top to bottom, a close up on Mandrill, House and Cameraman. (a) Underlying image x 0 . (b) Least-square estimate x LS (y). (c) Result with λ ⋆ for one level of decomposition J = 1, (d) for three levels of decomposition J = 3 using one global parameter, and (e) for three levels of decomposition J = 3 using one parameter per scales.

(S. 5 )

 5 This is a consequence of Proposition 1 and definition of weak differentiability since µ(y, •) is Lipschitz continuous independently of y. (S.6) By definition of SUGAR A FD {µ} in Proposition 1.

√ 2σ -Proposition 2 .

 2σ2 erf (µ 0 ) i -λ √ 2σ . It follows that V W [#{|Y | > λ}] = P i=1 p i (1p i ) P and hence lim ) * (ST(Y, λ) -Y ) = lim P →∞ V W 1 P λ#{|Y | > λ} = 0 .Consistency (i.e. convergence in probability) follows from traditional arguments by invoking Chebyshev inequality and using asymptotic unbiasedness (Theorem 1) and vanishing variance. Developing the mean squared error in terms of bias and variance give the first part of the proposition. Lemma 1 and the fact that lim ε→0 E W ∇ 2 { df {ST}}(Y, λ, ε) = ∇ 2 {df {ST}}(µ 0 , λ) concludes the second part.

Definition 2 (Definition 3 (

 23 Gauge). Let C be a non-empty closed convex set containing the origin. The gauge of C is the function γ C defined byγ C (y) = inf {ω > 0 \ y ∈ ωC} .As usual, γ C (y) = +∞ if the infimum is not attained. Polar set). Let C be a non-empty convex set. The set C • given byC • = z ∈ R N \ z, x 1 for all x ∈ Cis called the polar of C. C • is a non-empty closed convex set containing the origin, and if C is closed and contains the origin as well, C •• = C.

Lemma 2 .

 2 Let C be a non-empty closed convex set containing the origin. The following assertions hold.(i) γ C is a non-negative, closed, convex and positively homogenenous function.(ii) C is the unique closed convex set containing the origin such thatC = {y ∈ Y \ γ C (y) 1} .(iii) γ C is bounded and coercive if, and only if, C is compact and contains the origin as an interior point.(iv) The gauge of C and the support functionσ C • (y) = max z∈C • y, z coincide, i.e. γ C = σ C • .Proof. (i)-(ii) follow from [36, Theorem V.1.2.5]. (iii) is a consequence of [36, Theorem V.1.2.5(ii) and Corollary V.1.2.6]. (iv) [36, Proposition V.3.2.4].

Proposition 5 .

 5 Let C be a compact convex set containing the origin as an interior point, i.e. a convex body, and G = γ C is its gauge. For any θ > 0, θ ′ > 0 and any y ∈ Y, the following holdsProx θG (y) -Prox θ ′ G (y) L 2 |θθ ′ | for some constant L 2 > 0 independentof y, i.e. for any y, θ → Prox θG (y) is Lipschitz continuous on ]0, +∞[. Proof. From [35, Proposition 2.3(ii)], we have that for any y, the function θ → Prox θG (y) is such that Prox θG (y) -Prox θ ′ G (y) |θθ ′ | y -Prox θG (y) /θ .(33)Now, we have θG(y) = γ C/θ (y) = σ θC • (y) ,

= ∂ 2 ξ

 2 ξ = ∂ 1 ξ (ℓ) (y, θ)[δ] and D (ℓ) z k = ∂ 1 z (ℓ) k (y, θ)[δ].Hence, applying again the chain rule [28, Theorem 4 and Remark, Section 4.2.2] to the the sequence of iterates and using the fact that all involved mappings are Lipschitz, and D (ℓ) (y, θ) and J (ℓ)

Corollary 3 .

 3 As before, but now withD ξ = ∂ 1 ξ (ℓ) (y, θ)[δ], D (ℓ) x = ∂ 1 x(ℓ) (y, θ)[δ], D(ℓ)u = ∂ 1 u (ℓ) (y, θ)[δ], and D

Corollary 4 .= ∂ 1 ξ

 41 As before, but now with J (ℓ) (y, θ),J (ℓ) x = ∂ 1 x(ℓ) (y, θ), J (ℓ) u = ∂ 1 u (ℓ) (y, θ) and J (ℓ) x = Γ (ℓ) a (J (ℓ) a ) = J (ℓ)ξ . The chain rule completes the proof.

  † † y -Prox θG (y) = Proj θC • (y) L 2 γ C • Proj θC • (y) .Applying[START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] to γ C • , we gety -Prox θG (y) L 2 θγ θC • Proj θC • (y) L 2 θ ,

* (•) + stands for the Moore-Penrose pseudo-inverse.

† Another case which is trivial corresponds to G being the indicator function of a non-empty closed convex set, in which case L2 = 0.

* * Without impacting the optimal choice of λ, the curves have been rescaled for visualization purposes.

* This work has been supported by the European Research Council (ERC project SIGMA-Vision)