Charles-Alban Deledalle
email: charles-alban.deledalle@math.u-bordeaux.fr

Nicolas Papadakis
email: nicolas.papadakis@math.u-bordeaux.fr

Joseph Salmon
email: joseph.salmon@telecom-paristech.fr

Samuel Vaiter
email: samuel.vaiter@u-bourgogne.fr

CLEAR: Covariant LEAst-Square Refitting with Applications to Image Restoration

Keywords: inverse problems, variational methods, re-fitting, twicing, boosting, debiasing AMS subject classifications. 49N45, 65K10, 68U10

HAL is

Introduction.

Restoring an image from its single noisy and incomplete observation necessarily requires enforcing regularity or a model prior on the targeted properties of the sought solution. Regularity properties such as sparsity or gradient sparsity of an image are difficult to enforce in general, and notably lead to combinatorial and non-convex problems. When one is willing to guarantee such kinds of features, convex relaxation is a popular path. This is typically done using the 1 norm instead of the 0 pseudo-norm, as for the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] or the total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Nevertheless, such relaxations are well known to create solutions with a larger bias.

Typically, for the Lasso, using the 1 convex relaxation of the 0 pseudo-norm leads large coefficients to be shrunk towards zero.

For the total variation, the same relaxation on the jump amplitudes induces a loss of contrast in the recovered signal; see Figure 1.(a) for an illustration in this case.

In the Lasso case, a well known re-fitting scheme consists in performing a posteriori a least-square re-estimation of the non-zero coefficients of the solution. This post re-fitting technique became popular under various names in the literature: Hybrid Lasso [START_REF] Efron | Least angle regression[END_REF], Lasso-Gauss [START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF], OLS post-Lasso [START_REF] Belloni | Least squares after model selection in high-dimensional sparse models[END_REF], Debiased Lasso (see [START_REF] Lederer | Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions[END_REF][START_REF] Belloni | Least squares after model selection in high-dimensional sparse models[END_REF] for extensive details on the subject). For the anisotropic total-variation (aniso-TV), the same post re-fitting approach can be performed to re-estimate the amplitudes of the jumps, provided their locations have been correctly identified.

In this paper, we introduce a generalization of this re-fitting technique that aims at re-enhancing the estimation towards the data without altering the desired properties imposed by the model prior. The underlying reason is that if the user choose an estimator for which theoretical guarantees have been proven (as for the Lasso and TV), re-fitting should preserves these guarantees. To that end, we introduce the Covariant LEAst-square Re-fitting (CLEAR). Though this method was originally elaborated with 1 analysis problems in mind, it has the ability to generalize to a wider family, while in simple cases such as the Lasso or the aniso-TV, it recovers the classical post re-fitting solution described earlier. For instance, our methodology successfully applies to the Tikhonov regularization [START_REF] Tikhonov | On the stability of inverse problems[END_REF], the isotropic total-variation (iso-TV), the non-local means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], the block matching 3D (BM3D) [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF] and the Dual Domain Image Denoising (DDID) [START_REF] Knaus | Dual-domain image denoising[END_REF]. In common variational contexts, e.g., 1 -2 analysis [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] (encompassing the Lasso, the group Lasso [START_REF] Lin | Component selection and smoothing in multivariate nonparametric regression[END_REF][START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Obozinski | High-dimensional support union recovery in multivariate regression[END_REF], the aniso-and iso-TV), we show that our re-fitting technique can be performed with a complexity overload of about twice that of the original algorithm. In other cases, we introduce a scheme requiring about three times the complexity of the original algorithm.

While our covariant re-fitting technique recovers the classical post re-fitting solution in specific cases, the proposed algorithm offers more stable solutions. Unlike the Lasso post re-fitting technique, ours does not require identifying a posteriori the support of the solution, i.e., the set of non-zero coefficients. In the same vein, it does not require identifying the jump locations of the aniso-TV solution. Since the Lasso or the aniso-TV solutions are usually obtained through iterative algorithms, stopped at a prescribed convergence accuracy, the support or jump numerical identification might be imprecise (all the more for ill-posed problems). Such erroneous support identifications lead to results that strongly deviate from the sought re-fitting. Our covariant re-fitting jointly estimates the re-enhanced solution during the iterations of the original algorithm and, as a by product, produces more robust solutions.

This work follows a preliminary study [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF] that attempted to suppress the bias emerging from the choice of the method (e.g., 1 relaxation), while leaving unchanged the bias due to the choice of the model (e.g., sparsity). While the approach from [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF] hereafter referred to as invariant re-fitting -provides interesting results, it is however limited to a class of restoration algorithms that satisfy restrictive local properties.

In particular, the invariant re-fitting cannot handle iso-TV. In this case, the invariant re-fitting is unsatisfactory as it removes some desired aspects enforced by the prior, such as smoothness, and suffers from a significant increase of variance. A simple illustration of this phenomenon for iso-TV is provided in Figure 2.(d) where artificial oscillations are wrongly amplified near the boundaries.

While the covariant and the invariant re-fitting both correspond to the leastsquare post re-fitting step in the case of aniso-TV, the two techniques do not match for iso-TV. Indeed, CLEAR outputs a more relevant solution than the one from the invariant re-fitting. Figure 2.(e) shows the benefit of our proposed solution w.r.t. the (naive) invariant re-fitting displayed in Figure 2.

(d).

It is worth mentioning that the covariant re-fitting is also strongly related to boosting methods re-injecting useful information remaining in the residual (i.e., the map of the point-wise difference between the original signal and its prediction). Such approaches can be traced back to twicing [START_REF] Tukey | Exploratory data analysis[END_REF] and have recently been thoroughly investigated: boosting [START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF], Bregman iterations and nonlinear inverse scale spaces [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF][START_REF] Burger | Nonlinear inverse scale space methods[END_REF][START_REF] Xu | Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising[END_REF][START_REF] Osher | Sparse recovery via differential inclusions[END_REF], unsharp residual iteration [START_REF] Charest | On iterative regularization and its application[END_REF], SAIF-boosting [START_REF] Milanfar | A tour of modern image filtering: New insights and methods, both practical and theoretical[END_REF][START_REF] Talebi | How to saif-ly boost denoising performance[END_REF], ideal spectral filtering in the analysis sense [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF] and SOS-boosting [START_REF] Romano | Boosting of image denoising algorithms[END_REF] being some of the most popular ones. Most of these methods can be performed iteratively, leading to a difficult choice for the number of steps to consider in practice. Our method has the noticeable advantage that it is by construction a two-step one. Iterating more would not be beneficial (see, Subsection 4.2). Unlike re-fitting, these later approaches aim at improving the overall image quality by authorizing the re-enhanced result to deviate strongly from the original biased solution. In particular, they do not recover the aforementioned post re-fitting technique in the Lasso case. Our objective is not to guarantee the image quality to be improved but rather to generalize the re-fitting approach with the ultimate goal of reducing the bias while preserving the structure and the regularity of the original biased solution.

Interestingly, we have also realized that our scheme presents some similarities with the classical shrinking estimators introduced in [START_REF] Stein | Inadmissibility of the usual estimator for the mean of a multivariate normal distribution[END_REF], especially as presented in [START_REF] George | Minimax multiple shrinkage estimation[END_REF]. Indeed the step performed by CLEAR, is similar to a shrinkage step with a data-driven residual correction weight (later referred to as ρ in our approach, see Definition 14) when performing shrinkage as in [START_REF] George | Minimax multiple shrinkage estimation[END_REF]Section 3.1].

Last but not least, it is well known that bias reduction is not always favorable in terms of mean square error (MSE) because of the so-called bias-variance trade-off. It is important to highlight that a re-fitting procedure is expected to re-inject part of the variance, therefore it could lead to an increase of residual noise. Hence, the MSE is not expected to be improved by re-fitting techniques (unlike the aforementioned boostinglike methods that attempt to improve the MSE). We will show in our numerical experiments, that re-fitting is in practice beneficial when the signal of interest fits well the model imposed by the prior. In other scenarios, when the model mismatches the sought signal, the original biased estimator remains favorable in terms of MSE. Re-fitting is nevertheless essential in the latter case for applications where the image intensities have a physical sense and critical decisions are taken from their values.

Background models and notation.

We tackle the problem of estimating an unknown vector x 0 ∈ R p from noisy and incomplete observations (1)

y = Φx 0 + w ,
where Φ is a linear operator from R p to R n and w ∈ R n is the realization of a noisy random vector. This linear model is widely used in statistics and in imagery (e.g., for encoding degradations such as entry-wise masking, convolution, etc.). Typically, the inverse problem associated to ( 1) is ill-posed, and one should add additional information to recover at least an approximation of x 0 from y. We denote by x : y → x(y) the procedure that provides such an estimation of x 0 . In this paper, we consider a popular class of estimators relying on a variational formulation that involves a data fidelity term F (x, y) and a regularizing term G(x)1 :

(2)

x(y) ∈ argmin

x∈R p F (x, y) + G(x) .
Another kind of estimator can be defined as the output of an iterative algorithm (k, y) → x k , e.g., a numerical solver. For a chosen iteration k, we define the final estimator x(y) = x k . Such a framework includes for instance proximal splitting methods, as well as discretization of partial differential equations, though we do not investigate this latter road in details.

Notation. For a matrix

M , M + is its Moore-Penrose pseudo-inverse. For a (closed convex) set C, Π C is the Euclidean projection over C and ι C is its indicator function defined by ι C (u) = 0 if u ∈ C, +∞ otherwise.
For any integer d ∈ N * , we write [d] for the set {1, . . . , d}. For any subset I of [d], its complement in [d] is written I c . For any vector v, v I ∈ R |I| is the sub-vector whose elements are indexed by I ⊂ N and |I| is its cardinality. For any matrix M , M I is the sub-matrix whose columns are indexed by I. We denote respectively by Im[A] and Ker[A] the image space and the kernel space of an operator A. For any vector x ∈ R d , and q ∈ [1, +∞] we denote by ||x|| q its standard q norm, and by ||x|| 0 , the number of non-zero entries of x, i.e., ||x|| 0 = | {i ∈ [d] : The gray surface is a piece-wise affine mapping that models the evolution of x in an extended neighborhood of y. The light red affine plane is the model subspace, i.e., the set of images sharing the same jumps as those of the solution x(y). The red triangle is the restriction of the model subspace to images that can be produced by TV. Finally, the pink dot represents the re-fitting R inv x (y) as the orthogonal projection of y on M x(y).

x i = 0} |.

Main investigated estimators.

Here, we provide several canonical examples of estimators of the form (2) that help us illustrating our methodology.

The affine constrained least-squares, used when x 0 belongs to the affine subspace

C = b + Im[A] with b ∈ R p and A ∈ R p×n , are a particular case of (2) where F (x, y) = 1 2 ||y -Φx|| 2 2 and G(x) = ι C (x)
. The solution of minimum Euclidean norm is unique and given by: x(y) = b + A(ΦA) + (y -Φb) (see Appendix A.1). [START_REF] Tikhonov | On the stability of inverse problems[END_REF] (or Ridge regression [START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF]), used when ||Γx 0 || 2 is small with Γ ∈ R m×p , is an instance of (2) where F (x, y) = 1 2 ||y -Φx|| 2 2 and G(x) = λ 2 ||Γx|| 2 2 , for some parameters λ > 0. Provided Ker Φ ∩ Ker Γ = {0}, x(y) exists and is uniquely defined as: x(y) = (Φ Φ + λΓ Γ) -1 Φ y (see Appendix A.2).

The Tikhonov regularization

The hard-thresholding [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], used when Φ = Id and x 0 is sparse, is a solution of (2) where F (x, y) = 1 2 ||y -x|| 2 2 and G(x) = λ 2 2 ||x|| 0 , for some parameter λ > 0. The hard-thresholding operation writes: x(y) I = y I and x(y) I c = 0, where

I = {i ∈ [p] : |y i | > λ} (see Appendix A.3).
The soft-thresholding [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], used when Φ = Id and x 0 is sparse, is a solution of (2) where F (x, y) = 1 2 ||y -x|| 2 2 and G(x) = λ||x|| 1 for some parameter λ > 0. The soft-thresholding operation writes: x(y) I = y I -λ sign(y I ) and x(y) I c = 0, where I is defined as for the hard-thresholding (see Appendix A.4).

The 1 synthesis (or Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]), used when x 0 is sparse, is a solution of (2) where F (x, y) = 1 2 ||Φy -x|| 2 2 and G(x) = λ||x|| 1 for some parameter λ > 0. Provided the solution is unique (see for instance [START_REF] Tibshirani | The lasso problem and uniqueness[END_REF]), it reads [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] x(y) I = (Φ I ) + y I -λ((Φ I ) Φ I ) -1 s I and x(y)

I c = 0 ,
for almost all y and where I = supp(x(y)) = {i ∈ [p] : x(y) i = 0} is the support of x(y), and s I = sign(x(y) I ).

The 1 analysis, used when Γx 0 is sparse with Γ ∈ R m×p , is a solution of ( 2) where ( 4)

F (x, y) = 1 2 ||Φx -y|| 2 2 and G(x) = λ||Γx|| 1 , for some λ > 0 .
Provided Ker Φ ∩ Ker Γ = {0}, there exists a solution given implicitly, see [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF], as

(5) x(y) = U (ΦU ) + y -λU (U Φ ΦU ) -1 U (Γ ) I s I ,
for almost all y and where I = supp(Γx(y)) = {i ∈ [m] : (Γx(y)) i = 0} is called the Γ-support of the solution, s I = sign((Γx(y)) I ), U is a matrix whose columns form a basis of Ker[Γ I c ] and ΦU has full column rank. The anisotropic Total-Variation (aniso-TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is an instance of (4) where x 0 ∈ R p is identified to a b-dimensional signal, for which Γ = ∇ : R p → R p×b is the discrete gradient operator and ||∇x|| 1 = p i=1 ||(∇x) i || 1 . Aniso-TV promotes piecewise constant solutions with large constant regions and few sharp transitions. Here, I is the set of indexes where the solution has discontinuities. The 1 norm of the gradient field induces an anisotropic effect by favoring the jumps to be aligned with the canonical directions (so it favors squared like structures rather than rounded ones).

The 1 -2 analysis [START_REF] Lin | Component selection and smoothing in multivariate nonparametric regression[END_REF][START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Obozinski | High-dimensional support union recovery in multivariate regression[END_REF], used when Γx 0 is block sparse with Γ ∈ R p → R m×b , is a solution of ( 2) with ( 6)

F (x, y) = 1 2 ||Φx -y|| 2 2 and G(x) = λ||Γx|| 1,2 ,
for some parameter λ > 0 and where

||Γx|| 1,2 = m i=1 ||(Γx) i || 2 .
The isotropic Total-Variation (iso-TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is a particular instance of [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF] where

x 0 ∈ R p can be identified to a b-dimensional signal, for which Γ = ∇ : R p → R p×b and ||∇x|| 1,2 = p i=1 ||(∇x) i || 2 .
Like aniso-TV, it promotes solution with large constant regions, but some transition regions can be smooth (see, e.g., [START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF]), typically those with high curvature in the input image y, see Figure 2.(a)-(c). A major difference is that the 1 -2 norm induces an isotropic effect by favoring rounded like structures rather than squared ones. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], used when the image x 0 ∈ R p (with p = p 1 × p 2 ) is composed of many redundant patterns, is the solution of (2) for F (x, y) = 1 2 i,j w i,j ||P i x -P j y|| 2 2 and G(x) = 0, where P i is the linear operator extracting a patch (i.e., a small window) at pixel i of size (2b + 1) 2 . The index i ∈ [p 1 ] × [p 2 ] spans the image domain and j -i ∈ [-s, s] 2 spans a limited search window. The weights w i,j are defined as w i,j = ϕ ||P i y -P j y|| 2 2 where the kernel ϕ : R + → [0, 1] is a decreasing function. Assuming periodical boundary conditons, its solution is given as [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] x(y) i = j wi,j y j j wi,j with wi,j =

The (block-wise) non-local means

k w i-k,j-k ,
where k ∈ [-d, d] 2 spans the patch domain (see Appendix A.5).

Limitations.

It is important to note that some of the previously introduced estimators are known to suffer from systematic drawbacks.

For instance, the Lasso contracts large coefficients towards 0 by the shift λ((Φ I ) Φ I ) -1 s I whose expression simplifies to λ sign(y I ) for the soft-thresholding. In the same vein, the 1 analysis contracts the signal towards Ker[Γ] by the shift λU (U Φ ΦU ) -1 U (Γ ) I s I . In the particular case of aniso-TV, this last quantity is the well known systematic loss of contrast of TV: a shift of intensity on each piece depending on its surrounding and the ratio between its perimeter and its area (see [START_REF] Strong | Edge-preserving and scale-dependent properties of total variation regularization[END_REF] for a thorough study for the 1D case). A simple illustration is provided in Figure 1 where TV is used for denoising a 1D piece-wise constant signal in [0,192] and damaged by additive white Gaussian noise (AWGN) with a standard deviation σ = 10. Even though TV has perfectly retrieved the support of ∇x 0 with one more extra jump, the intensities of some regions are systematically under-or overestimated. Iso-TV also suffers from a systematic loss of contrast, as illustrated in Figure 2 (c).

In the following we investigate possible solutions to reduce such artifacts.

3. Invariant least-square re-fitting. As mentioned earlier, practitioners have realized that a systemic contraction affect estimators like the Lasso and the aniso-TV. In the Lasso case, a simple remedy (presented in the introduction) is to perform a posteriori a least-square re-fitting step of the non-zero coefficients, i.e., constrained to the support I of the Lasso solution x(y), given by argmin

x ; supp(x) ⊆ I 1 2 ||Φx -y|| 2 2 . ( 8
)
In this section, we present a re-fitting procedure, discussed in [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF], that generalizes this approach to a broad family of estimators.

3.1.

Re-fitting through model subspace least-squares. We investigate a re-fitting procedure well suited for estimators almost everywhere differentiable (a.e. differentiable). It relies on the notion of model subspace, which requires Jacobian matrix computations. From now on, we consider only estimators y → x(y) from R n to R p and a.e. differentiable.

Many estimation procedures rely on a structural prior of the data. Such structures include smoothness, sparsity, auto-similarity or the fact that the signal is piece-wise constant. Such priors can be captured by the following notion of model subspace.

Definition 1. The model subspace associated to an a.e. differentiable estimator x is defined at almost all points y ∈ R n by the affine subspace of R p [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] M x(y) = x(y) + Im [J x(y)] ,

where J x(y) ∈ R p×n is the Jacobian matrix of x taken at y.

The model subspace captures what is linearly invariant through x w.r.t. small perturbations of y, typically, for the Lasso, it will encode the set of signals sharing the same support. In order to generalize the re-fitting step, it is thus natural to cast it as a constrained optimization procedure preserving the model subspace.

Definition 2. The invariant re-fitting associated to an a.e. differentiable estimator y → x(y) is given for almost all y ∈ R n by [START_REF] Cornillon | Recursive bias estimation for multivariate regression smoothers[END_REF] R inv x (y) = x(y) + J(ΦJ) + (y -Φx(y)) ∈ argmin

x∈M x(y) 1 2 ||Φx -y|| 2 2 ,
where J = J x(y) is the Jacobian matrix of x at the point y. In the following, we use the notation J when no ambiguity is possible.

Remark 3. We only consider F (x, y) = The affine constrained least-squares have everywhere the same Jacobian matrix J = A(ΦA) + and its affine model subspace is M x(y) = b+Im[A(ΦA) ] (as Im[M + ] = Im[M ]). Taking C = R p , n = p, A = Id and b = 0, leads to x(y) = Φ + y whose model subspace is M x(y) = Im[Φ ], reducing to R p when Φ has full column rank. In this case, the invariant re-fitting is R inv x (y) = x(y). The Tikhonov regularization has everywhere the same Jacobian matrix J = (Φ Φ + λΓ Γ) -1 Φ and everywhere the same affine model subspace M x(y) = Im[J]. If follows that R inv x (y) = J(ΦJ) + y. In particular, when Φ has full column rank, M x(y) = R p and R inv x (y) = Φ + y. The soft-and hard-thresholding share a.e. the same Jacobian matrix given by J = Id I ∈ R p×|I| . Their model subspace reads as M x(y) = Im[Id I ] = Im[J] and the invariant re-fitting is for both the hard-thresholding.

The 1 synthesis (or Lasso) has for Jacobian matrix J = Id I (Φ I ) + a.e. and when Φ I has full column rank, it shares the same model subspace M x(y) = Im[Id I ] = Im[J] as the hard-and soft-thresholding. Its invariant re-fitting is in this case R inv x (y) = Id I (Φ I ) + y. While the Lasso systematically underestimates the amplitude of the signal by a shift λId I ((Φ I ) t Φ I ) -1 s I , the re-fitting R inv x (y) is free of such a contraction. The 1 analysis has for Jacobian matrix J = U (ΦU ) + a.e. (since the Γ-support, I and the sign s I are a.e. stable w.r.t. small perturbations [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF]). It results that the model subspace reads as M x(y) = Ker[Γ I c ] = Im [J]. The generalization of the Lasso re-fitting leads to the least-square estimator constrained on the Γ-support of x(y): [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF] R inv x (y) = U (ΦU ) + y . For aniso-TV denoising (i.e., Φ = Id and Γ = ∇), the model subspace is the space of images whose jumps are included in those of the solution x(y). The re-fitting procedure R inv x is the projector Π Im[J] = U U + that performs a piece-wise average of its input on each plateau of the solution (see Figure 1.(a) for an illustration in 1D).

The 1 -2 analysis admits a.e. the following Jacobian operator [START_REF] Vaiter | The degrees of freedom of partly smooth regularizers[END_REF]

J = U (U Φ ΦU + λU Γ ΩΓU ) -1 U Φ y (12)
where Ω :

z ∈ R m×b → 1 ||(Γx(y))i||2 z i -z i , (Γx(y))i ||(Γx(y))i||2 (Γx(y))i ||(Γx(y))i||2 if i ∈ I, 0 otherwise,
with U and I defined exactly as for the 1 analysis case. Note that Eq. ( 12) is well founded as soon as Ker[Γ I c ] ∩ Ker Φ = {0}. For weaker conditions, see [50, example 26]. In the particular case where ΦU has full column rank, the model subspace matches with the one of the 1 analysis, i.e., M x(y) = Ker[Γ I c ] = Im[J], and the re-fitting is also given by Eq. ( 11), hence R inv x (y) = U (ΦU ) + y. As a consequence, for the iso-TV denoising, i.e., when Φ = Id and Γ = ∇, the invariant re-fitting step consists again in performing a piece-wise average on each plateau of the solution.

The non-local means has a Jacobian matrix with a complex structure [START_REF] Van De Ville | Nonlocal means with dimensionality reduction and surebased parameter selection[END_REF][START_REF] Duval | A bias-variance approach for the nonlocal means[END_REF]. In particular, computing the projection on the model subspace is challenging in this case and so is the computation of the invariant re-fitting. Note that a greedy procedure was proposed in [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF] to compute the invariant re-fitting. of all plateaus towards the data. While the invariant re-fitting acts properly for the 1 analysis estimator, it is however less appealing in other scenarios. Figure 4,(a)-(c) gives another illustration of the invariant re-fitting of a 2D Tikhonov masking example (with Φ a diagonal matrix with 0 or 1 elements on the diagonal and Γ = ∇). While the dynamic of the Tikhonov solution x(y) has been strongly reduced, the re-fitting R inv x (y) re-fits clearly the solution towards the original intensities. However, such a re-fitting is not satisfactory as it does not preserve the smoothness of the solution x(y).

Results and limitations.

In fact, the model subspace captures only what is linearly invariant through x w.r.t. small perturbations of y. This includes the support of the solution for the iso-TV, and the absence of variations inside Im[J] ⊥ for the Tikhonov regularization. In particular, it fails at capturing some of the desirable relationships between the entries of y and the entries of x(y), what we call the covariants. These relationships typically encode some of the local smoothness and non-local interactions between the entries of the solution x(y). Such crucial information is not encoded in the linear model subspace, but interestingly the Jacobian matrix captures by definition how much the entries of x linearly varies w.r.t. all the entries of y. This is at the heart of the covariant re-fitting defined in the next section and, for comparison, it produces the results given in Figure 2.

Covariant LEast-square Re-fitting (CLEAR).

The objective of this section is to present our main contribution, the introduction of the covariant re-fitting procedure (CLEAR). We particularly aim at solving the issues raised in Subsection 3.3. Toward this goal, we put a stronger emphasis on the first-order behavior of the original estimator by imposing the conservation of its Jacobian, at least locally. To define the re-fitting procedure we first need to introduce a procedure which, loosely speaking, takes as input the original estimator and a guess of Φx 0 and outputs a new estimator with some targeted properties. We next define our Covariant LEast-square re-fitting (CLEAR) by choosing a specific guess for Φx 0 , see Definition 14.

4.1.

Local approach and desired properties for a suitable re-fitting. In this subsection, our objective is to define, from the original estimator x and a guess z ∈ R n of Φx 0 , a new estimator D x,z : R n → R p that satisfies several desirable properties and shares with x some first-order properties. After-wise, we will consider the choice z = y, and the resulting estimator is going to be our covariant re-fitting R x. We are now equipped to introduce such a guess based re-fitting.

Definition 5. Let x : R n → R p be differentiable at z ∈ R n . An estimator D x,z : R n → R p is a guess based covariant least-square re-fitting of x for z, if (13) D x,z ∈ argmin h∈H ||Φh(z) -z|| 2 2 ,
where H is the set of maps h : R n → R p satisfying, for all y ∈ R n , 1. Affine map:

h(y) = Ay + b for some A ∈ R p×n , b ∈ R p , 2. Covariant preserving: J h (z) = ρJ x(z) for some ρ ∈ R, 3. Coherent map: h(Φx(z)) = x(z).
Definition 5 is natural as it states that a guess based re-fitting of x for z should be, in prediction, as close as possible to z. Of course, it should satisfy some extra conditions. First, the estimator should be easy to compute, and so we choose a first order approximation, leading to a locally affine estimator. Second, as highlighted in Subsection 3.3, the relative variation of the original estimator w.r.t. the input should be preserved to capture, not only the invariant features of the estimator but also its first-order behavior, capturing both its singularities and smoothness. Third, applying a re-fitting step to the prediction obtained by the original estimator at z should not modify it. The purpose of re-fitting is to be close to y, while also preserving the structure of x(z). Hence, if y = Φx(z), the result should be unaltered.

The next theorem provides a unique closed form expression for D x,z (y).

Theorem 6. Let x be an estimator from R n to R p differentiable at z ∈ R p . Then, for δ = z -Φx(z), the guess based covariant least-square re-fitting, defined in Definition 5, exists, is unique if ΦJδ = 0, and is given by

(14) D x,z (y) = x(z) + ρJ(y -Φx(z)) where ρ =    ΦJδ, δ ||ΦJδ|| 2 2 if ΦJδ = 0 , 1 otherwise ,
where J = J x(z) is the Jacobian matrix of x at the point z.

Proof. Let h be a mapping satisfying properties 1., 2. and 3. in the previous definition. Observe that properties 1. and 2. of the set H together ensures that the estimator is of the form h(y) = ρJy +b for some ρ ∈ R and b ∈ R p . Plugging condition 3. gives that b = (Id -ρJΦ)x(z), hence h(y) = x(z) + ρJ(y -Φx(z)). Reciprocally, it is easy to see that any estimator of the form h(y) = x(z) + ρJ(y -Φx(z)) satisfies properties 1., 2. and 3. It thus remains to find ρ.

Remark that Problem (13) can be recast as a one-dimensional problem

min ρ∈R ||Φ(x(z) + ρJ(z -Φx(z))) -z|| 2 2 = ||(Id -ρΦJ)(Φx(z) -z)|| 2 2 (15)
whose unique solution, if ΦJ(Φx(z) -z) = 0 is given by ( 14) and ρ = 1 otherwise.

The case where ΦJ is an orthogonal projector leads to interesting properties for instance when x is associated to the constrained least-squares, the Lasso or aniso-TV.

Proposition 7. If ΦJ is an orthogonal projector, then ρ = 1.
Proof. If ΦJ is an orthogonal projector, it follows that ||ΦJδ|| 2 2 = ΦJδ, ΦJδ = (ΦJ) ΦJδ, δ = ΦJδ, δ . Injecting this in [START_REF] Deledalle | Contrast re-enhancement of total-variation regularization jointly with the Douglas-Rachford iterations[END_REF] gives ρ = 1.

Statistical interpretation. For a random vector with expectation Φx 0 and finite second order moment, bias and covariance of D x,z are given in closed form.

Proposition 8. Let Y be a random vector in R n such that E[Y ] = Φx 0 , and Cov[Y ] = Σ ∈ R n×n . Then y → D x,z (y) satisfies E[D x,z (Y )] -x 0 = (Id -ρJΦ)(x(z) -x 0 ) , (16) Cov[D x,z (Y ), Y ] = ρJΣ , (17) Cov[D x,z (Y )] = ρ 2 JΣJ , ( 18
)
where the cross covariance is Cov[X, Y ] = E[XY ] -E[X]E[Y ] ,

for any random column vectors X and Y (not necessarily of the same size), and Cov

[Y ] = Cov[Y, Y ].
Proof. The first equality is a direct consequence of the linearity of the expectation operator. The second equality arises from the following

E[(x(z) + ρJ(Y -Φx(z)))Y ]-E[(x(z) + ρJ(Y -Φx(z)))]E[Y ] (19) = ρJ E[Y Y ] -E[Y ]E[Y ]
since J and ρ are constant w.r.t. y as they depend only on the guess z. The third equation follows the same sketch by expanding the expression of Cov[D x,z (Y )].

Proposition 8 provides a closed form expression for the bias, the cross-covariance and the covariance of D x,z . These quantities are much more intricate to derive for a non-linear estimator x. Nevertheless, the next corollary shows how these quantities relate to those of the first order Taylor expansion of the original estimator x.

Corollary 9. Let T x,z (y) be the tangent estimator of x at z ∈ R n defined as

(20) T x,z (y) = x(z) + J(y -z) .
Let Y be a random vector in

R n such that E[Y ] = Φx 0 and Cov[Y ] = Σ. Then y → T x,z (y) and y → D x,z (y) satisfy E[T x,z (Y )] -x 0 = (x(z) -x 0 ) + J(Φx 0 -z) , (21) Cov[D x,z (Y ), Y ] = ρCov[T x,z (Y ), Y ] , (22) Cov[D x,z (Y )] = ρ 2 Cov[T x,z (Y )] , (23) Corr[D x,z (Y ), Y ] = Corr[T x,z (Y ), Y ] , (24) Corr[D x,z (Y )] = Corr[T x,z (Y )] , ( 25
)
where the cross correlation matrix is defined as Corollary 9 is essential in this work as it states that, by preserving the Jacobian structure, D x,z (Y ) cannot depart from the tangent estimator of x at z in terms of (cross-)correlations. As a consequence, one can expect that they only differ in terms of expectation, i.e., in terms of bias. The next propositions state that when ΦJ is a projector, the bias in prediction is guaranteed to be reduced by our re-fitting.

Corr[X, Y ] i,j = Cov[X, Y ] i,j / Cov[X] i,i Cov[Y ] j,j ,
Proposition 10. Let Y be a random vector of R n such that E[Y ] = Φx 0 . Assume ΦJ is an orthogonal projector, then y → D x,z (y) satisfies ||Φ(E[D x,z (Y )] -x 0 )|| 2 ||Φ(E[T x,z (Y )] -x 0 )|| 2 . ( 26
)
Proof. As ΦJ is an orthogonal projector, by virtue of Proposition 7, ρ = 1, then

||Φ(E[T x,z (Y )] -x 0 )|| 2 2 = ||Φ(x(z) + J(Φx 0 -z) -x 0 )|| 2 2 (27) = ||Φ(x(z) + JΦ(x 0 -x(z)) -x 0 )|| 2 2 + ||ΦJ(Φx(z) -z)|| 2 2 = ||Φ(E[D x,z (Y )] -x 0 )|| 2 2 + ||ΦJδ|| 2 2
which concludes the proof.

Proposition 10 is a bit restrictive as it requires ΦJ to be a projector. Nevertheless, this assumption can be relaxed when z satisfies a more technical assumption as shown in the next proposition.

Proposition 11. Let Y be a random vector of R n such that E[Y ] = Φx 0 . Let ρ 0 = δ0, ΦJδ0 ||ΦJδ0|| 2 2 and δ 0 = Φ(x 0 -x(z)). Assume there exists α ∈ [0, 1] such that ρ -ρ 0 ρ 0 √ 1 -α , ( 28
)
and ||ΦJ(δ -δ 0 )|| 2 2 + 2 δ 0 , ΦJ(δ -δ 0 ) -α δ 0 , ΦJδ 0 2 ||ΦJδ 0 || 2 2 . ( 29
)
Then, y → D x,z (y) satisfies ||Φ(E[D x,z (Y )] -x 0 )|| 2 ||Φ(E[T x,z (Y )] -x 0 )|| 2 . ( 30
)
Proof. It follows from Proposition 8 and Corollary 9 that

||Φ(E[D x,z (Y )]-x 0 )|| 2 = ||(Id -ρΦJ)δ 0 || 2 and ||Φ(E[T x,z (Y )] -x 0 )|| 2 = ||δ 0 + ΦJ(δ -δ 0 )|| 2 . Subsequently, we get that Equation (30) holds true if ||(Id -ρΦJ)δ 0 || 2 2 ||δ 0 + ΦJ(δ -δ 0 )|| 2 2 , (31) i.e., ρ 2 ||ΦJδ 0 || 2 2 -2ρ δ 0 , ΦJδ 0 ||ΦJ(δ -δ 0 )|| 2 2 + 2 δ 0 , ΦJ(δ -δ 0 ) . ( 32
)
Using Assumption (29), a sufficient condition for Equation [START_REF] Lin | Component selection and smoothing in multivariate nonparametric regression[END_REF] to hold is

ρ 2 ||ΦJδ 0 || 2 2 -2ρ δ 0 , ΦJδ 0 2 + α δ 0 , ΦJδ 0 2 ||ΦJδ 0 || 2 2 0 . ( 33
)
The roots of this polynomial are given by (1 ± √ 1 -α)ρ 0 , which concludes the proof.

Remark 12. Remark that requiring [START_REF] Lebrun | An analysis and implementation of the bm3d image denoising method[END_REF] is quite natural as it states that ρ should be close enough to the optimal ρ 0 minimizing the discrepancy with regards to Φx 0 (i.e., minimizing ||Φh(z) -Φx 0 || 2 2 for h ∈ H defined as in Definition 5). While the condition (29) sounds more technical, it however holds true in several interesting cases. For instance, when z = Φx 0 , Assumption (29) holds true as it would read 0 -α δ 0 , ΦJδ 0 2 /||ΦJδ 0 || 2 2 (since δ = δ 0 and ρ = ρ 0 ). Another case of interest is when ΦJ is an orthogonal projector for which [START_REF] Lederer | Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions[END_REF] holds true as it would read

||ΦJδ|| 2 2 -||ΦJδ 0 || 2 2 -||ΦJδ 0 || 2 2 (using that ρ = ρ 0 = 1, •, ΦJ• = ||ΦJ • || 2 2
, and choosing α = 1). Hence, Proposition 11 recovers Proposition 10.

Remark 13. Using the same sketch of proof as Proposition 11, the condition

|(ρ -ρ 0 )/ρ 0 | 1 is sufficient to get ||Φ(E[D x,z (Y )] -x 0 )|| 2 ||Φ(x(z) -x 0 )|| 2 .
In other words, even though ρ has a relative error of 100% w.r.t. ρ 0 , the estimator y → D x,z (y) still reduces the bias of the constant estimator y → x(z). This result remains valid when comparing y → D x,z (y) to the pseudo-oracle estimator y → x(z) + J(y -Φx 0 ), with the notable difference that they moreover share the same correlation structure.

While it is difficult to state a general result, we can reasonably claim from Proposition 10, Proposition 11 and Remark 12 that the bias tends to be reduced by our re-fitting providing ΦJ is almost a projector (i.e., has eigenvalues concentrated around 0 and 1) and/or z is not too far from Φx 0 . In such cases, the estimator D x,z can be considered as a debiasing procedure of x, in the sense that it reduces the bias of T x,z while preserving its correlation structure (according to Corollary 9).

Definitions and properties. Using

) 34
where δ = y -Φx(y) and J = J x(y) is the Jacobian matrix of x at the point y.

Figure 5 gives a geometrical interpretation of CLEAR for a denoising task in dimension p = 3. One can observe that if Y varies isotropically, so will its projection on the model subspace. Contrarily, the tangent estimator at a guess z can present an anisotropic behavior along the model subspace, and the guess based re-fitting, which is closer to z, will respect this anisotropy in order to capture the local regularity of x. Finally, the covariant re-fitting is obtained from the guess based re-fitting at z = y. For clarity, we assumed that M x(z) = M x(y) on this illustration.

Remark 15. The covariant re-fitting performs an additive correction of x(y) with a fraction of the directional derivative Jδ in the direction of the residual δ.

Remark 16. Observe that in Definition 14, the value of ρ varies when y varies, contrary to the map y → D x,z (y) for which ρ is constant. Note that the mapping y → D x,z (y) is affine, but not the map y → R x(y). Note that, as a consequence, the statistical interpretations given in the previous section do not hold for R x(y) even though they shed some light on its behavior.

Why not iterating the procedure? The re-fitting procedure computes x2 = x1 + ρJ(y -Φx 1 ), with x1 = x(y). One may wonder if it is beneficial to iterate the process as xk+1 = xk + ρJ(y -Φx k ) (in the same vein as [START_REF] Tukey | Exploratory data analysis[END_REF][START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF][START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF][START_REF] Talebi | How to saif-ly boost denoising performance[END_REF][START_REF] Romano | Boosting of image denoising algorithms[END_REF]). Consider a denoising problem Φ = Id with Tikhonov or iso-TV, for which J is symmetrical and x(y) ∈ Im x (y) (by virtue of Remark 4), i.e., the invariant re-fitting. Reminding the artifacts illustrated in Figure 3.(e), this is not satisfying.

An interesting property of R x is the fact that it belongs to the model subspace of x as stated in the following proposition.

Proposition 17. Let y → x(y) be an a.e. differentiable estimator. Then for almost all y ∈ R n , one has R x(y) ∈ M x(y).

Proof. As M x(y) = x(y)+Im J and ρJ(y-Φx(y)) ∈ Im J, the proposition holds.

The case where ΦJ is an orthogonal projector, leads to interesting properties that will be of interest regarding some estimators considered in Subsection 2.2.

Proposition 18. Suppose that ΦJ is an orthogonal projector. Then, R x(y) = x(y) + J(y -Φx(y)), and, ΦR x(y) = ΦR inv x (y). Proof. By virtue of Proposition 7, ρ = 1 and then R x(y) = x(y) + J(y -Φx(y)). The fact that ΦR x(y) = ΦR inv x (y) comes from the fact that ΦJ(ΦJ) + = ΦJ. The next proposition provides, when JΦ satisfies a fixed point formulation, an expression of R x(y) that will be useful for efficient computations of the re-fitting as discussed in Section 5, a notable example being the iso-TV regularization.

Proposition 19. Assume that JΦx(y) = x(y). Then, the covariant re-fitting reads R x(y) = (1 -ρ)x(y) + ρJy.

Proof. We have R x(y) = x(y) + ρJ(y -Φx(y)) = x(y) + ρJy -ρJΦx(y), and since JΦx(y) = x(y) by assumption, this concludes the proof.

Interestingly, the next theorem shows that the condition JΦx(y) = x(y) is met provided x(y) is solution of a variational problem with a 1-homogeneous regularizer.

Theorem 20. Let x(y) be the unique a.e. differentiable solution of [START_REF] Osher | Sparse recovery via differential inclusions[END_REF] x(y) = argmin The affine constrained least-squares, the 1 synthesis, the 1 -2 analysis, aniso-TV and iso-TV, are solutions of a variational problem with F being differentiable and G being 1-homogeneous. As a consequence, Theorem 20 shows that the aforementioned methods satisfy JΦx(y) = x(y), and hence R x(y) = (1 -ρ)x(y) + ρJy.

Examples of re-fitting procedures.

We now exemplify the previous definitions for the wide class of variational estimators introduced in Subsection 2.2.

The affine constrained least-squares have for Jacobian matrix J = A(ΦA) + . In this case, ΦJ = ΦA(ΦA) + is an orthogonal projector, ρ = 1 and the covariant re-fitting coincides with the invariant one and reads R x(y) = R inv x (y) = x(y). The Tikhonov regularization has for Jacobian J = (Φ Φ + λΓ Γ) -1 Φ and in this case ρ depends on the residual δ and the re-fitting reads as the weighted sum R x(y) = (1 + ρ)x(y) -ρJΦx(y).

The soft-and hard-thresholding, used with Φ = Id, have the Jacobian matrix J = Id I . As a consequence ΦJ = Id I is an orthogonal projection and the covariant re-fitting coincides with the invariant one, namely the hard-thresholding itself.

The 1 synthesis has for Jacobian matrix J = Id I (Φ I ) + where Φ I has full column rank. As for the thresholding, ΦJ = ΦU (ΦU ) + is an orthogonal projection and the covariant re-fitting reads R x(y) = R inv x (y). The 1 analysis has for Jacobian matrix J = U (ΦU ) + . Again, ΦJ = ΦU (ΦU ) + is an orthogonal projection and the covariant re-fitting reads R x(y) = R inv x (y). The 1 -2 analysis has the Jacobian operator given in Eq. [START_REF] Darbon | Fast nonlocal filtering applied to electron cryomicroscopy[END_REF], which applied to a vector d ∈ R n is a solution of the following problem

Jd ∈ argmin x ; supp(Γx) ⊆ I 1 2 ||Φx -d|| 2 2 + λ 2 ω(Γx) , ( 36
)
where

ω : z ∈ R m×b → i∈I 1 ||(Γx(y)) i || 2 ||z i || 2 2 -z i , (Γx(y)) i ||(Γx(y)) i || 2 2 .
Remark that ω(ΓJd) = 0 only if (ΓJd) i is co-linear to (Γx(y)) i , for all i ∈ I. For iso-TV, it means that the level lines of Jd must be included in the ones of x(y). Moreover, by virtue of Theorem 20, one has JΦx(y) = x(y) and hence R x(y) = (1-ρ)x(y)+ρJy. As a consequence, unlike the invariant re-fitting of x(y), the covariant re-fitting is constrained to be faithful to the regularity of x(y), since it enforces the discontinuities of Jd to be co-linear to (Γx(y)) I . This is especially important where the iso-TV solution presents transitions with high curvature. Such an appealing behavior of the covariant re-fitting explains the results observed in Figure 2.(e) and Figure 3.(f).

The non-local means has an intricate Jacobian matrix. Nevertheless, its directional derivative has a simpler expression given, for any direction d ∈ R n , by [START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF] Jd = j w i,j y j + wi,j d j -x(y) i j w i,j j wi,j with w i,j = k w i-k,j-k , where x(y) is defined in Eq. ( 7) and [START_REF] Romano | Boosting of image denoising algorithms[END_REF] w i,j = 2 P i y -P j y, P i d -P j d ϕ ||P i y -P j y|| 2 , with ϕ the a.e. derivative of the kernel function ϕ. Subsequently, the covariant refitting is obtained from its general form with two steps, by computing first x(y), and applying next the Jacobian to the direction d = y -x(y).

Covariant re-fitting in practice.

This section details the computation of CLEAR for standard algorithms. We first recall some properties of two different differentiation techniques that allow computing some image of J jointly with x(y).

Algorithmic differentiation.

Following [START_REF] Deledalle | Stein unbiased gradient estimator of the risk (SUGAR) for multiple parameter selection[END_REF], we consider restoration algorithms whose solutions x(y) = x k are obtained via an iterative scheme of the form [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] x k = γ(a k ) , a k+1 = ψ(a k , y) .

Here, a k ∈ A is a sequence of auxiliary variables, ψ : A × R n → A is a fixed point operator in the sense that a k converges to a , and γ : A → R p is non-expansive (i.e., ||γ(a 1 ) -γ(a 2 )|| ||a 1 -a 2 ||, ∀a 1 , a 2 ∈ A) entailing x k converges to x = γ(a ).

As a result, for almost all y and for any direction d ∈ R n , the directional derivatives D k x = J xk (y)d and D k a = J a k (y)d can be jointly obtained with x k and a k as (40)

         x k = γ(a k ) , a k+1 = ψ(a k , y) , D k x = Γ a D k a , D k+1 a = Ψ a D k a + Ψ y d ,
where Γ a = ∂γ(a) ∂a a k , Ψ a = ∂ψ(a,y) ∂a a k and Ψ y = ∂ψ(a k ,y) ∂y y

. Interestingly, in all considered cases, the cost of evaluating Γ a , Ψ a and Ψ y is about the same as the one of evaluating γ and ψ. As a result, the complexity of ( 40) is of about twice the complexity of [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. In practice, Γ a , Ψ a and Ψ y can be implemented either thanks to their closed form expression or in a black box manner using automatic differentiation. The later has been well studied and we 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 study.

Finite difference based differentiation.

Another strategy is to approximate directional derivatives by finite differences, for any d ∈ R n and ε > 0, as

J x(y)d ≈ x(y + εd) -x(y) ε . ( 41
)
As a result, the complexity of evaluating ( 41) is also twice the complexity of (39) since x must be evaluated at both y and y + εd. The main advantage of this method is that x can be used as a black box, i.e., without any knowledge on the underlying algorithm that provides x(y). For ε small enough, it performs as well as the approach described in (40) (with x(y) = x k ) that requires the knowledge of the derivatives. Indeed, if y → x(y) is Lipschitz-continuous, then (41) converges to [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF] when ε → 0 (by virtue of Rademacher's theorem and [22, Theorem 1-2, Section 6.2]). This implies that the value ε can be chosen as small as possible (up to machine precision) yielding an accurate approximation of J x(y)d. This finite difference approach has been used in many fields, and notably for risk estimation, see e.g., [START_REF] Ye | On measuring and correcting the effects of data mining and model selection[END_REF][START_REF] Shen | Adaptive model selection[END_REF][START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF]. We will apply this black box strategy on state-of-the-art algorithms in Subsection 6.5.

Two-step computation for the general case.

In the most general case, the computation of the covariant re-fitting, given by

(42) R x(y) = x(y) + ρJ(y -Φx(y)) with ρ = ΦJδ, δ ||ΦJδ|| 2 2
and δ = y -Φx(y) , requires to evaluate sequentially x(y) and J(y -Φx(y)). With finite difference differentiation, two steps are required. First x(y) must be computed with the original algorithm and next J(y -Φx(y)) is obtained by finite difference (41) on the direction of the residual d = y -Φx(y). Once J(y -Φx(y)) is computed, ρ can be evaluated and subsequently [START_REF] Stein | Inadmissibility of the usual estimator for the mean of a multivariate normal distribution[END_REF]. The overall complexity is about twice that of the original algorithm producing x(y).

With algorithmic differentiation, as J(y -Φx(y)) depends on x(y), the original iterative scheme (39) must be run first. In the second step, J(y -Φx(y)) is obtained with the differentiated version (40) on the direction of the residual d = y -Φx(y). As a result, x(y) is computed twice, first by [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], next by [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF]. It leads to an overall complexity about three times the one of the original algorithm. Nevertheless, in several cases, one can avoid the first step by running (40) only once. The computations of x(y) and Jy are then sufficient to compute the re-fitting R x(y).

One-step computation

As a result, in the case of algorithmic differentiation, (40) can be run once to get R x(y) since using d = y provides directly x(y), Jy and subsequently ρ. Compared to the two step approach, the complexity of the re-fitting reduces to about twice the one of the original step from [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Recall, that the condition JΦx(y) = x(y) is met for the Lasso, the Generalize Lasso, aniso-TV and iso-TV, hence they can be re-enhanced with a complexity being twice the one of their original algorithm.

Example on a primal-dual solver for 1 analysis.

In this section we instantiate Algorithm [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF] to the case of the primal-dual sequence of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. By dualizing the 1 analysis norm x → λ||Γx|| 1 , the primal problem (4) can be reformulated, with x = x(y), as the following saddle-point problem [START_REF] Talebi | How to saif-ly boost denoising performance[END_REF] (z , x ) ∈ arg max z∈R m min

x∈R p 1 2 Φx -y 2 2 + Γx, z -ι B λ (z) ,
where z ∈ R m is the dual variable, and

B λ = {z ∈ R m : ||z|| ∞ λ} is the ∞ ball.
Problem [START_REF] Talebi | How to saif-ly boost denoising performance[END_REF] can be efficiently solved using the primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. By taking στ < 1/ Γ 2 2 , θ ∈ [0, 1] and initializing (for instance,)

x 0 = v 0 = 0 ∈ R p , z 0 = 0 ∈ R m , the algorithm reads (45)    z k+1 = Π B λ (z k + σΓv k ) , x k+1 = (Id + τ Φ Φ) -1 x k + τ (Φ y -Γ z k+1 ) , v k+1 = x k+1 + θ(x k+1 -x k ) ,
where the projection of z over B λ is done component-wise as

(46) Π B λ (z) i = z i if |z i | λ , λ sign(z i ) otherwise .
The sequence x k converges to a solution x of the 1 analysis problem [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

It is easy to check that the primal-dual sequence defined in ( 45) can be written in the general form considered in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], see for instance [START_REF] Deledalle | Stein unbiased gradient estimator of the risk (SUGAR) for multiple parameter selection[END_REF]. As a result, we can use the algorithmic differentiation based strategy described by [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF] as follows: for the initialization x0 = ṽ0 = 0 ∈ R p , z0 = 0 ∈ R m , and for β = 0, as

               z k+1 = Π B λ (z k + σΓv k ) , x k+1 = (Id + τ Φ Φ) -1 x k + τ (Φ y -Γ z k+1 ) , v k+1 = x k+1 + θ(x k+1 -x k ) , zk+1 = Π z k +σΓv k (z k + σΓṽ k ) , xk+1 = (Id + τ Φ Φ) -1 xk + τ (Φ y -Γ zk+1 ) , ṽk+1 = xk+1 + θ(x k+1 -xk ) , ( 47
)
where Π z (z

) i = zi if |z i | λ + β , 0 otherwise .
Recall that the re-fitting is R x k (y) = xk , since JΦ is an orthogonal projector.

Remark that the algorithmic differentiation of ( 45) is exactly [START_REF] Tikhonov | On the stability of inverse problems[END_REF] for β = 0, hence, xk = R x k (y). However, if one wants to guarantee the convergence of the sequence xk towards R x(y), one needs a small β > 0 as shown in the next theorem. In practice, β can be chosen as the smallest available positive floating number.

Theorem 21. Assume that x satisfies (5) with ΦU full-column rank 2 . Let α > 0 be the minimum non zero value 3

of |Γx | i for all i ∈ [m]. Choose β such that ασ > β > 0. Then, the sequence xk = R x k (y) defined in (47) converges to the re-fitting R x(y) of x(y) = x .
The proof of this theorem is postponed to Appendix C.

A similar result was obtained in [START_REF] Deledalle | Contrast re-enhancement of total-variation regularization jointly with the Douglas-Rachford iterations[END_REF] when solving the 1 analysis problem (4) with the Douglas-Rachford splitting algorithm described in [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF].

Example for a primal-dual solver for

1 -2 analysis. The algorithm for the 1 -2 analysis regularization can be derived with the exact same considerations as for the 1 analysis case. The only difference in the application of the primal dual algorithm comes from the non linear operation (46) that now reads, for z ∈ R m×b , as

(48) Π iso B λ (z) i = z i if ||z i || 2 λ , λ zi ||zi||2 otherwise .
Algorithm Non-local means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] and its directional derivative

Inputs:

noisy image y, direction d, noise standard-deviation σ Parameters:

half search window width s, half patch width b, kernel function ϕ Outputs:

x = x(y) and x = J x(y)d

Initialize W ← ϕ(2σ 2 (2b + 1) 2 ) 1p 1 ×p 2 (add weights for the central pixels [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF]) Initialize Wy ← ϕ(2σ 2 (2b + 1) 2 ) y (accumulators for the weighted sum)

Initialize W ← 0p 1 ×p 2 Initialize W y ← ϕ(2σ 2 (2b + 1) 2 ) d for k ∈ [-s, s] 2 \ {0, 0} do Compute e ← [ (y -S k (y))
Compute e ← [ 2(y -S k (y))(d -S k (d)) ] κ Compute w ← [ e ϕ (e) ] κ Update W ← W + w Update W y ← W y + w S k (f ) + wS k (d) end for Compute x ← Wy/W (weighted mean) Compute x ← (W y -W x )/W
Fig. 6. Pseudo-algorithm for the computation of the non-local means and its Jacobian in a direction d. All arithmetic operations are element wise, S k is the operator that shift all pixels in the direction k, is the discrete convolution operator, and

κ ∈ R p 1 ×p 2 is such that κ i,j = 1 if (i, j) ∈ [-b, b] 2 , 0 otherwise.

It follows that the algorithmic differentiation strategy reads as

               z k+1 = Π iso B λ (z k + σΓv k ) , x k+1 = (Id + τ Φ Φ) -1 x k + τ (Φ y -Γ (z k+1 )) , v k+1 = x k+1 + θ(x k+1 -x k ) , zk+1 = Π iso z k +σΓv k (z k + σΓṽ k ) , xk+1 = (Id + τ Φ Φ) -1 xk + τ (Φ y -Γ zk+1 ) , ṽk+1 = xk+1 + θ(x k+1 -xk ) , ( 49
)
where Π iso z (z

) i = zi if ||z i || 2 λ + β , λ ||zi||2 zi -zi , zi ||zi||2 zi ||zi||2 otherwise .
Unlike the 1 case, the re-fitted solution is not xk itself, but following Subsection 5.4, it can be obtained at iteration k as

R x k (y) = (1 -ρ)x k + ρx k with ρ = Φ(x k -x k ), y -Φx k ||Φ(x k -x k )|| 2 2 . ( 50
)

Example for the non-local means.

In this section, we specify the update rule [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF] to an acceleration of the non-local means inspired from [START_REF] Darbon | Fast nonlocal filtering applied to electron cryomicroscopy[END_REF][START_REF] Deledalle | Non-local methods with shape-adaptive patches (NLM-SAP)[END_REF]. We use the procedure of [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF] to correctly handle the central pixel. Again, one can check that this implementation can be written in the general form considered in the update rule [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], where the fixed point solution is obtained directly at the first iteration. The pseudo-code obtained with the algorithmic differentiation scheme, as described in [START_REF] Salmon | On two parameters for denoising with non-local means[END_REF] is given in Figure 6. All variables with suffix correspond to the directional derivative obtained by using the chain rule on the original variables. This fast computation relies on that all convolutions can be computed with integral tables, leading to a global complexity in O(s 2 n), for both the computation of the estimator x(y) and its directional derivative Jd. Recall that the covariant re-fitting is obtained from its general form with two steps, by computing first x(y), and applying next the proposed pseudo-code in the direction d = y -x(y).

Numerical experiments and comparisons with related approaches.

Here, we first give illustrations of our CLEAR method on toy image restoration problems. Then, through quantitative results in terms of PSNR 4 and SSIM [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], we numerically evaluate the re-fitting, discuss its benefit and limitations in several scenarios and compare this method with popular approaches from the literature.

Denoising with isotropic total-variation (iso-TV).

Figure 7 gives an illustration of our covariant re-fitting of the 2D iso-TV, where λ has been chosen large enough in order to highlight the behavior of the re-fitting. We apply it for the denoising (i.e., Φ = Id) of an 8bits piece-wise constant image damaged by AWGN with standard deviation σ = 20, known as the Shepp-Logan phantom. As discussed earlier, iso-TV introduces a significant loss of contrast [START_REF] Strong | Edge-preserving and scale-dependent properties of total variation regularization[END_REF], typically for thin detailed structures, which are re-enhanced in our result. The residuals x(y) -x 0 and R x(y) -x 0 highlight that our re-fitting technique re-enhances efficiently the attenuated structure while leaving the lost structures unchanged. Nevertheless, after re-fitting, some small residuals around the edges appear. In fact, in the vicinity of edges, iso-TV finds (barely visible) discontinuities that are not in accordance with the underlying image. This creates an overload of small constant regions. When re-fitting is performed, such regions are re-fitted to the noisy data, and they become barely visible artifacts. In other words, the re-fitting has re-enforced the presence of a modeling problem, resulting to an increase of residual variance, that iso-TV had originally compensated by attenuating the amplitudes.

De-masking with isotropic total-variation (iso-TV).

Figure 8 gives another illustration of our covariant re-fitting of the 2D iso-TV used for the restoration of an approximate 8bits piece-wise constant image damaged by AWGN with σ = 20, known as Boat. The operator Φ is a random mask removing 40% of pixels. Again, iso-TV introduces a significant loss of contrast, typically for thin details such as the contours of the objects, which are re-enhanced in our re-fitting result.

Inspecting the map of residuals in Figure 8.(e)-(f) illustrates that our re-fitting technique eliminates most of the bias to the price of a small variance increase. This is clearer when looking at the mast and the ropes. While the mast was preserved by iso-TV and re-enhanced by our re-fitting, the ropes remain lost for both methods. 9 gives an illustration of our re-fitting procedure for the non-local means algorithm used in a denoising problem of the 8bits image Pirate, enjoying many repetitive patterns and damaged by AWGN with σ = 20. We choose a regularizing kernel ϕ(•) = exp(•/h), h > 0 that leads to strong smoothing in order to highlight the behavior of the re-fitting. Our refitting technique provides favorable results: many details are enhanced compared to the standard method. This reveals that the non-local means is actually able to well capture the repetitions of many patterns but this information is not used properly to create a satisfying result. The re-fitting produces a sharper result, by enforcing the correct use of all the structures identified by patch comparisons. The maps of residuals Figure 9.(e)-(f) highlight that our re-fitting technique suppresses efficiently this dull effect while it preserves the model originally captured by patch redundancy. Again the suppression of this phenomenon is counter balanced by an increase of the residual variance, prominent where the local patch redundancy assumption is violated.

Denoising with non-local means. Figure

In these examples, the overall residual norm is clearly reduced by the re-fitting because the amount of reduced bias surpasses the increase of residual variance. This favorable behavior depends on the internal parameters of the original estimator acting on the bias-variance trade-off, as we investigate next.

6.4.

A bias-variance analysis for the covariant re-fitting. Previous experiments have revealed that while CLEAR tends to reduce the bias, it increases (as expected) the residual variance. It is therefore important to understand under which conditions the bias-variance trade-off is in favor of our re-fitting technique.

Figure 10 illustrates the evolution of performance, measured in terms of mean squared error (MSE), of both aniso-TV and its re-fitting version as a function of the regularization parameter λ. Two images are considered: Cameraman, an approximate piece-wise constant image (top), and a truly piece-wise constant image (bottom).

This experiment highlights that optimal results for both approaches are not reached at the same λ value. Visual inspection of the optima shows that due to the bias, the optimal solution of aniso-TV is reached for a λ value promoting a model subspace that is not in accordance with the underlying signal: typically the presence of an overload of (barely visible) transitions in homogeneous areas. These transitions become clear when looking at the re-fitted version where each small region is re-fitted on the noisy data, revealing an excessive residual variance. Conversely, the optimal λ value for the re-fitting seems to retrieve the correct model, i.e., with transitions that are closely in accordance with the underlying signal. Comparing their relative performance, when both are used at their own optimal λ, reveals that our re-fitting brings a significant improvement if the underlying image is in fact piece-wise constant.

Figure 11 provides a similar illustration of the evolution of performance for the non-local means and its re-fitted version as a function of the smoothing parameter h of the kernel function ϕ(•) = exp(•/h). Two images are considered: Lady, a crude approximation of an image with redundant patterns (top), and Fingerprint, an image with many redundant patterns (bottom). Similar conclusions can be made from this

.t. h. Two values of h are selected corresponding to (c) re-fitting for a sub-optimal

h, (d) original for a sub-optimal h, (f) original for the optimal h, (g) re-fitting for the optimal h. experiment. In particular, comparing their relative performance, when both are used at their own optimal h value, seems to demonstrate that the re-fitting brings an improvement when most patches of the underlying image are redundant.

While it is difficult to make a general statement, we can reasonably claim from these experiments that re-fitting is all the more relevant in terms of MSE than the underlying image x 0 is in agreement with the retrieved subspace model M x(y). In other words, re-fitting is safe when the original restoration technique was chosen appropriately w.r.t. the underlying image of interest. Beyond MSE performance, the Two values of γ are selected corresponding to (c) re-fitting for a sub-optimal γ, (d) original for a sub-optimal γ, (f) original for its optimal γ, (g) re-fitting for its optimal γ. re-fitted results at their optimal parameters choice might nevertheless be preferable (as assessed by the SSIM values): intensities and contrasts being recovered better.

6.5. Behavior on more sophisticated filters. We focus here on two filters BM3D [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF] and DDID [START_REF] Knaus | Dual-domain image denoising[END_REF] reaching state-of-the-art results in As mentioned in the previous paragraph, re-fitting is all the more relevant that the first estimate is obtained with high smoothing strength. It is thus important to compare the original filter with its re-fitted version for varying smoothing parameters. For simplicity, we have considered only one smoothing parameter γ > 0 for these two algorithms.

For DDID, we use the authors implementation, and choose to multiply by γ the two inner parameters γ f and γ r [START_REF] Knaus | Dual-domain image denoising[END_REF]. For BM3D, we use the implementation of [START_REF] Lebrun | An analysis and implementation of the bm3d image denoising method[END_REF], and choose to multiply by γ the inner parameters λ hard 3D (in the hard thresholding step) and σ (in the Wiener filtering step). Unlike previous experiments, the re-fitted results are obtained here by finite difference as discussed in Subsection 5.2.

Figure 12 illustrates the evolution of performance, measured in terms of MSE, of BM3D, DDID and their re-fitted versions as a function of the smoothing parameter γ. BM3D is studied on the Monarch image and DDID on the Lena image. Similar conclusions can be made from this experiment, re-fitting reaches its optimal performance for a larger smoothing parameter. Because the original estimator is nearly unbiased, re-fitting becomes challenging only when the original solution is over-smoothed, otherwise the gain in terms of bias is too small to compensate the loss in terms of variance. Given this high smoothing strength, some tiny structures have been lost, and thus the optimal re-fitting does not reach as good performance as the optimal original solution.

However, even tough the MSE is not necessarily improved, our re-fitting solutions present less artifacts (as confirmed by the very small gap in terms of SSIM values), see for instance the stripes of the Monarch or the left cheek of Lena. In fact, in order to recover details with low signal to noise ratio, the optimal original estimators authorize the apparition of low contrasted oscillating features. Nevertheless, a few of these oscillations tends to amplify some noise structures, hence, explaining these artifacts. In contrast, the optimal re-fitting prefers loosing such details rather than taking the risk of creating arbitrary structures and is thus more reliable.

Regarding our previous discussion, we believe that re-fitting would be nevertheless beneficial in terms of MSE for the class of images that are promoted by such restoration techniques. Characterizing this class of solutions for these methods is a very challenging topic out of the scope of this paper.

6.6.

Comparisons with other techniques devoted to the 1 case. We detail hereafter two different alternative strategies devoted to re-enhance the solution of the 1 analysis regularization.

Iterative hard-thresholding. As shown earlier, the hard-thresholding is the refitted version of the soft-thresholding. Given an iterative solver (k, y) → x k composed of linear and soft-thresholding (such as primal-dual algorithms), one could consider replacing all soft-thresholding by hard-thresholding while keeping linearities unchanged: a technique often referred to as "iterative hard-thresholding" [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF]. Unfortunately, such techniques only provide convergence to a local minimum of the 0 -regularized problem, and they do not converge to the sought re-fitting R x(y).

Co-support identification based post re-fitting. Another solution referred to as post re-fitting, and studied in, e.g., [START_REF] Efron | Least angle regression[END_REF][START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF][START_REF] Belloni | Least squares after model selection in high-dimensional sparse models[END_REF][START_REF] Lederer | Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions[END_REF][START_REF] Belloni | Least squares after model selection in high-dimensional sparse models[END_REF], consists in identifying the (co-)support I = {i : (Γx(y)) i = 0} and solving a least-square problem constrained to {x : (Γx) I c = 0}, typically with conjugate gradient descent. However, x(y) is usually obtained thanks to a converging sequence x k , and unfortunately, supp(Γx k ) can be far from supp(Γx(y)) even though x k can be made arbitrarily close to x(y). Such erroneous support identifications lead to results that strongly deviates from R x(y).

Figure 13 provides a comparison of our re-fitting method with the two other approaches mentioned earlier for aniso-TV used on an 8bits image (peppers), damaged by AWGN with σ = 20. The iterative hard-thresholding approach does not preserve the model space: transitions are not localized at the same positions as in the original solution and suspicious oscillations are created. The post re-fitting and our approach have both improved x(y) by enhancing each piece and preserving the location of transitions. Our method is nevertheless more stable than support identification which produces many errors due to wrong co-support identification.

Figure 14 provides another illustration highlighting the problem of support identification in an ill-posed problem. It consists of an 8bits image (Cameraman) damaged by a Gaussian blur of 2px and AWGN with σ = 20. Again, while aniso-TV reduces the contrast, the re-fitting recovers the original amplitudes and leaves discontinuities unchanged. Post re-fitting offers comparable results to ours except for suspicious oscillations due to wrong co-support identification.

In contrast with the support identification, CLEAR does not require the identification of the co-support, nor the identification of the model subspace M x(y). This is appealing since the co-support of x(y) is difficult to identify, particularly in the analy-sis context. Being computed along the iterations of the original algorithm, i.e., jointly with x(y), our re-fitting strategy also provides more stable solutions.

Comparisons with boosting strategies.

We detail hereafter other popular alternatives designed to re-enhance results of an arbitrary estimator.

Twicing and boosting. The boosting iterations introduced in [START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF] is a simple approach that consists in re-injecting to the current solution xk a filtered version of its residual y -Φx k . The idea is that if parts of the signal where lost at iteration k, they might be retrieved in the residual. Given x0 = 0, the iterations reads

xk+1 = xk + x(y -Φx k ) . ( 51
)
The first iterate is x1 = x(y), and x2 is known as the twicing estimate [START_REF] Tukey | Exploratory data analysis[END_REF]. Such approaches are popular in non-parametric statistics, e.g., for kernel smoothing [START_REF] Cornillon | Recursive bias estimation for multivariate regression smoothers[END_REF].

When k increases, its bias tends to decreases while its variance increases, see [START_REF] Talebi | How to saif-ly boost denoising performance[END_REF]. In denoising (i.e., when Φ = Id) with a linear estimator x(y) = W y (e.g., the Tikhonov regularization), we get xk = (Id -(Id -W ) k )y, for k > 0. In particular, the twicing reads as x2 = (2W -W 2 )y. In this linear case, the covariant re-fitting reads as R x(y) = (W + ρW -ρW 2 )y and coincides with the twicing when ρ = 1.

The Unsharp Residual Iteration (URI) is a related approach studied in [START_REF] Charest | On iterative regularization and its application[END_REF], that reads xk+1 = xk + x(y) -x(Φx k ) and which is also equivalent in the linear case.

Iterative Bregman refinement. In [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF], the authors proposed an iterative procedure, originally designed to improve iso-TV results, given by

xk+1 = x y + k i=1 (y -Φx k ) . ( 52
)
Unlike boosting that iteratively filters the residual, the idea is to filter a modified version of the input y amplified by adding the sum of the residuals. When Φ = Id and x(y) = W y, the iterative Bregman refinement reads as xk = (Id -(Id -W ) k )y and coincides with boosting (we refer to [START_REF] Burger | Nonlinear inverse scale space methods[END_REF][START_REF] Xu | Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising[END_REF][START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF][START_REF] Osher | Sparse recovery via differential inclusions[END_REF] for related approaches).

SOS-boosting.

In [START_REF] Romano | Boosting of image denoising algorithms[END_REF], the authors follow a similar idea by iteratively filtering a strengthened version of the input y. Their method, named Strengthen Operate Substract boosting (SOS-boosting), originally proposed for Φ = Id, performs iteratively the following update

xk+1 = τ x(y + αΦx k ) -(τ α + τ -1)x k . ( 53
)
where α and τ are two real parameters. The first one controls the emphasis of the solution (and the convergence of the sequence), while the second one controls the rate of convergence. When Φ = Id and x(y) = W y, the SOS refinement with τ = 1 reads as xk = W y+α(W -Id)x k-1 , and in particular, for k = 2, we get x2 = (W -αW +αW 2 )y which coincides with our covariant re-fitting for the choice α = -ρ. For all considered estimators, we have always observed ρ > 0, contrarily to [START_REF] Romano | Boosting of image denoising algorithms[END_REF], where α > 0 is implicitly assumed. Hence, we cannot concludes that the two models match in a specific linear setting. Another difference is that while we provide an automatic way to compute ρ (see Equation ( 43)), the α parameter of the SOS-boosting must be tuned by the practitioner, a possibly cumbersome task, e.g., when using cross-validation on a fix dataset of images and/or for varying noise levels. (e,f,g,h) From top to bottom, boosting [START_REF] Tukey | Exploratory data analysis[END_REF][START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF], Bregman iterations [START_REF] Burger | Nonlinear inverse scale space methods[END_REF] and SOS-boosting [START_REF] Romano | Boosting of image denoising algorithms[END_REF] at respectively 1, 2, 5 and 20 iterations. SAIF-boosting. As described in [START_REF] Milanfar | A tour of modern image filtering: New insights and methods, both practical and theoretical[END_REF], the diffusion of a filter consists in iteratively re-applying the filter to the current estimate xk+1 = x xk . The authors of [START_REF] Talebi | How to saif-ly boost denoising performance[END_REF] noticed that, unlike the boosting method of [START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF], the bias of this estimator increases and its variance decreases with k. As a consequence, the authors suggest mixing the two approaches by deciding at each iteration between performing a boosting or a diffusion step. To that end, they proposed a plug-in risk estimator that crudely estimates the MSE from a pre-filtered version of y. This approach is in fact applied locally on image patches, and is referred to as the Spatially Adapted Iterative Filter (SAIF)-boosting. Unlike the other techniques, the SAIF-boosting cannot be used as a black-box. Indeed, it requires to perform an eigen decomposition of x locally for each patch of y. This can be efficiently done for some kernel-based averaging filters, but can be very challenging for arbitrary estimators, such as for instance iso-TV.

Though we have compared the expressions of boosting approaches with our refitting in the case of linear estimators, it is worth mentioning that boosting approaches are scarcely used in this case. Boosting appears more relevant in the non-linear case, since the successive re-application of a non-linear estimator x allows to recover parts of the signal that were lost at former iterations. Nevertheless, to boost the solution, the internal parameters of x may need to be re-adapted at each iteration, leading to cumbersome tuning of parameters in practice. Unlike boosting methods, a re-fitting approach should not modify the regularity and the structure of the first estimate.

where R a diagonal matrix with diagonal elements R ii = λ/|y i | for i ∈ I and 0 otherwise. One can check that if 1 n ∈ Ker[Γ] (which holds true for aniso-TV), then W is row stochastic. However, as seen in Equation ( 5), the quantity U (U Φ ΦU ) -1 U (Γ ) I R is the term for the systematic contraction of the 1 analysis regularization (this simplifies to λ/|y i | for the soft-thresholding). As a consequence, the bias cannot be corrected by a single application of W . The Jacobian J = U (ΦU ) + , which is again row stochastic, is free of this contraction term. Therefore our covariant re-fitting gets rid of this bias term after one single application of the Jacobian J.

Figure 15 and Figure 16 provide a comparison of our covariant re-fitting with the three first aforementioned boosting approaches, on two 8bits images (Flinstones and Barbara) damaged by AWGN with σ = 20, respectively for iso-TV and for non-local means. Note that the inner parameters of iso-TV and non-local means are fixed for all approach. The α and τ parameters of the SOS-boosting approach have been tuned to offer the most satisfying results, even though, we did not observe a significant impact in the iso-TV case. As expected, our covariant re-fitting provides results re-enhanced towards the amplitudes of the noisy inputs. In contrast, boosting approaches do not systematically re-fit towards the original amplitudes. While CLEAR preserves the structural content and smoothness of the original solution, the boosting approaches re-inject structural contents that were not originally preserved, and can present a large amount of residual noise. In these experiments, the original estimators were significantly biased, but for smaller parameter choices, re-injecting structural contents would have improved more the PSNR than re-fitting the amplitudes only. In fact, boosting techniques are more relevant to improve the quality of near unbiased estimates, while re-fitting techniques corrects only for the in-accuracy of biased but precise estimates (i.e., with low estimation variance).

Figure 17 illustrates this point on a 8bits image (House), damaged by AWGN with σ = 20 and σ = 60 for iso-TV. Unlike the previous comparison, the approaches are compared for different values of the regularization parameter λ. We have chosen to perform one single iteration for all boosting methods such that they share the same computation time with our re-fitting. Studying the MSE as a function of λ, reveals that SOS-boosting and Twicing indeed improve on the MSE for lower values of the regularization parameter (i.e., , for near unbiased estimators). The Bregman iteration and our re-fitting improve on the biased estimate for higher values of the regularization parameter λ. The optimal re-fitting has higher MSE than the biased estimate and all other boosting approaches (tuned at their own optimal λ value). This was expected as Fig. 10, 11 and 12 have already shown that the MSE of the optimal original estimator is better than the one of the optimal CLEAR. As a consequence any "boosting" approach that barely modifies the original estimator will necessarily present lower MSE than CLEAR, in particular for Twicing and SOS-boosting in this context. However the SSIM values are really close for all compared methods.

Quantifying the quality of a debiasing is a difficult task. The main issue is that the noisy image itself has a null bias, and any method that over-fits to the noise tends to have a small bias too. Then debiasing cannot be only defined in terms of bias reduction, otherwise x(y) = Φ + y would be an optimal rule. Debiasing should be defined as a technique minimizing the bias while preserving properties of the original solution, e.g., the model subspace (invariant refitting), or the cross-correlations (covariant refitting). Because of the lack of such a quantitative criterion, we have displayed the results obtained by each method at their optimal λ value for visual comparisons. At a relatively low noise level σ = 20, all approaches are on par. For σ = 60, fundamental differences are revealed. SOS-boosting and Twicing appear close to the original biased estimate, while the Bregman iteration and our re-fitting lead to a significant contrast re-enhancement with larger homogeneous regions.

While in practice the results obtained by the Bregman iteration may look preferable, our re-fitting has a methodological benefit. If the user chooses a biased estimator for which theoretical guarantees on the model space have been proven (e.g., Lasso, TV, etc.), our CLEAR solution enjoys the same properties, whereas they can be lost with boosting. For instance, for TV, one can characterize the preserved/lost objects with respect to their perimeter, area and the chosen λ. Because the model subspace is preserved, the same objects can be preserved/lost after applying CLEAR. This cannot be guaranteed for boosting approaches even though they can offer a lower approximation error. This behavior in terms of MSE is well known for the Lasso (for instance, soft-thresholding is known to achieve lower MSE than hard-thresholding) but even though, refitting is commonly performed because it preserves the guarantees proved for the original estimator while reaching lower bias. Given the large efforts made by the community to develop theoretical results for solutions of variational problems, we believe it is important that our scheme inherits such guarantees.

Conclusion.

We have introduced a generalization of the popular least-square re-fitting, originally aimed at reducing the systematic contraction of the Lasso.

Together with this generalization, a generic implementation has been given for a wide class of ill-posed linear problem solvers. This implementation requires a computational overload of at most a factor three compared to the original solver; a factor that can even be reduced to two for most popular estimators in image processing.

While the classical re-fitting is inspired by the standard Lasso debiasing step (i.e., least-square re-fitting on the estimated support), our generalization leverages the Jacobian of the estimator and does not rely on the notion of support, and its unstable identification. In particular, the proposed implementation only requires chain rules and differentiating the considered solver, and in practice it has also the benefit of increased stability compared to naive implementations.

For estimators such as Tikhonov regularization, total-variation non-local means, BM3D or DDID, numerical experiments have demonstrated the efficiency of the CLEAR technique in retrieving correct intensities while respecting the structure of the original biased estimator. Moreover, it has been shown in practice that re-fitting is beneficial when the underlying signal structure is well captured by the original estimator. Otherwise, re-fitting leads to too simplistic approximations, typically reflecting an inaccurate prior model. In other words, if the considered estimator is adequate with respect to the application context, then re-fitting is recommended.

We have highlighted the importance in distinguishing boosting approaches from the re-fitting one. In particular, re-fitting should be preferred in applications where the content of the original solution must be preserved. While boosting approaches are mostly used to enhance near unbiased estimators (typically coming from combinatorial or non-convex problems), the re-fitting is all the more relevant for estimators that present biases. For instance, re-fitting is essential for estimators solution of a convex problem that require a large bias correction to accurately retrieve the content of the signal of interest, a canonical example being the isotropic total-variation. Nevertheless, we believe that the notion of Jacobian based re-fitting could be of interest for boosting applications and we leave this to future work.

Acknowledgments.

The authors would like to acknowledge the financial support of the GDR 720 ISIS (Information, Signal, Image et viSion). This study has also been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux -CPU (ANR-10-IDEX-03-02).

Fig. 1 .

 1 Fig. 1. (a) Solutions of 1D-TV and our re-fitting on a noisy signal. (b) Illustration of the invariant re-fitting in a denoising problem of dimension p = 3.The gray surface is a piece-wise affine mapping that models the evolution of x in an extended neighborhood of y. The light red affine plane is the model subspace, i.e., the set of images sharing the same jumps as those of the solution x(y). The red triangle is the restriction of the model subspace to images that can be produced by TV. Finally, the pink dot represents the re-fitting R inv x (y) as the orthogonal projection of y on M x(y).

Fig. 2 .

 2 Fig. 2. (a) A piece-wise constant signal. (b) Its noisy version. (c) Solution of iso-TV with λ = 10 on the noisy signal. (d) Solution of the invariant re-fitting of iso-TV. (e) Solution of the covariant re-fitting of iso-TV. Red points indicate locations where the discrete gradient is non-zero.

( a )Fig. 3 .

 a3 Fig. 3. (a) A noisy image y = x 0 + w where x 0 contains six overlapping squares. (b) The solution x(y) of aniso-TV with λ = 10, and(c) its invariant re-fitting R inv x (y) (which coincides in this case with the covariant one R x(y)). (d) The solution x(y) of iso-TV with λ = 10, (e) its invariant re-fitting R inv x (y) and (f) its covariant one R x(y).

Figure 1 .

 1 (a) illustrates the invariant re-fitting in the case of a 1D total-variation denoising example ( 1 analysis estimator). It recovers the jumps of the underlying signal (adding an extra one), but systematically under-estimates their amplitudes. As expected, re-fitting re-enhances the amplitudes

Fig. 4 .

 4 Fig. 4. (a) A noisy and incomplete image y = Φx 0 + w where Φ is a masking operator encoding missing pixel values (in red) and x 0 is a smooth signal. (b) Tikhonov regularization x(y) with Γ = ∇ and λ = 20, (c) its invariant re-fitting R inv x (y) and (d) its covariant one R x(y).

Figure 1 .

 1 (b) gives a geometrical interpretation in dimension p = 3 of the model subspace and the invariant re-fitting. The model subspace is represented as the tangent plane of x at y and its re-fitting is the projection of y on this plane. All elements of this plane share the same jumps with the solution x(y). Figure3.(a)-(c) gives a similar illustration for a 2D aniso-TV denoising example.

Figure 2 .

 2 (c),(d) and Figure3.(d),(e) give two illustrations of the invariant re-fitting of a 2D iso-TV denoising example. As for aniso-TV, the invariant re-fitting is the projection of y on the space of signals whose jumps are located at the same position as those of x(y). But unlike the anisotropic case, x(y) is not piece-wise constant. Instead of being composed of distinct flat regions, it reveals smoothed transitions with dense supports (referred to as extended supports in[START_REF] Chambolle | Geometric properties of solutions to the total variation denoising problem[END_REF]), see Figure2.(c). As a consequence, the invariant re-fitting re-introduces a large amount of the original noisy signal in these smooth but non-constant areas, creating the artifacts observed on Figure2.(d) and Figure3.(e).

  (e), Figure 3.(f) and Figure 4.(d).

  for any random column vectors X and Y (not necessarily of the same size), and Corr[Y ] = Corr[Y, Y ]. Proof. The first relation holds from the expression of T x,z and that J does not depend on y. It follows that Cov[T x,z (Y ), Y ] = JΣ and Cov[T x,z (Y )] = JΣJ . These, jointly with Proposition 8, conclude the proof.

  D x,z defined in Theorem 6, we can now give an explicit definition of CLEAR as R x(y) = D x,y (y). Definition 14 (CLEAR). The Covariant LEast-square Re-fitting associated to an a.e. differentiable estimator y → x(y) is, for almost all y ∈ R n , given by R x(y) = x(y) + ρJ(y -Φx(y)) with ρ =

Fig. 5 .

 5 Fig. 5. Geometrical illustration of the covariant re-fitting in a denoising problem of dimension p = 3. We assume that M x(z) = M x(y) for the sake of clarity. The gray surface is the manifold modeling the evolution of x in an extended neighborhood of y. The light red affine plane is the model subspace tangent at z. The ellipses represent the positive-definite symmetric covariance matrices of some random vectors, as defined in Proposition 8 and Corollary 9.

  [J] (see examples of Subsection 3.2). The sequence will converge if and only if J(y -xk ) vanishes, i.e., xk must converge to J + Jy + ζ with ζ ∈ Ker[J]. By construction, xk ∈ Im[J] = Ker[J ] ⊥ = Ker[J] ⊥ , hence ζ = 0. Moreover, as J is symmetrical and x(y) ∈ Im[J], the quantity J + Jy coincides with JJ + y = R inv

xF

  (y -Φx) + G(x) , with F , G being convex and G being 1-homogeneous. Then, JΦx(y) = x(y) a.e. . The proof of Theorem 20 is postponed to Appendix B.

  for specific cases. When x(y) fulfills the assumption JΦx(y) = x(y) of Proposition 19, the covariant re-fitting reads as (43) R x(y) = (1 -ρ)x(y) + ρJy with ρ = Φ(Jyx(y)), y -Φx(y) ||Φ(Jy -x(y))|| 2 2 .

PSNR: 23 Fig. 7 .

 237 Fig. 7. (a) Noisy y = x 0 + w and (d) noise-free x 0 . (b) Iso-TV x(y) with λ = 3 and (e) the residual x(y) -x 0 . (c) Our re-fitting R x(y) and (f) the residual R x(y) -x 0 .

PSNR: 9 . 51 ,Fig. 8 .

 9518 Fig. 8. (a) Partial and noisy y = Φx 0 + w (red indicates missing pixels) and (d) noise-free x 0 . (b) Iso-TV x(y) with λ = 3 and (e) the residual x(y) -x 0 . (c) Our re-fitting R x(y) and (f) the residual R x(y) -x 0 .

PSNR: 22 Fig. 9 .

 229 Fig. 9. (a) Noisy y = x 0 + w and (d) noise-free x 0 . (b) Non-local means x(y) with h = 3 and (e) the residual x(y) -x 0 . (c) Our re-fitting R x(y) and (f) the residual R x(y) -x 0 .

Fig. 10 .

 10 Fig. 10. Experiment with aniso-TV: (top) poorly piece-wise constant case. (bottom) pure piecewise constant case. (a) Noise-free x 0 . (b) Noisy y = x 0 + w. (e) MSE of x(y) and its re-fitting R x(y) w.r.t. λ. Two values of λ are selected corresponding to (c) re-fitting for a sub-optimal λ, (d) original for a sub-optimal λ, (f) original for the optimal λ, (g) re-fitting for the optimal λ.

Fig. 11 .

 11 Fig. 11. Experiment with non-local means: (top) moderate patch redundancy case. (bottom) high patch redundancy case. (a) Noise-free x 0 . (b) Noisy y = x 0 + w. (e) MSE of x(y) and refitting R x(y) w.r.t. h. Two values of h are selected corresponding to (c) re-fitting for a sub-optimalh, (d) original for a sub-optimal h, (f) original for the optimal h, (g) re-fitting for the optimal h.

Fig. 12 .

 12 Fig. 12. (top) Experiments with BM3D [11, 28] and with DDID [27] (bottom). (a) Noise-free x 0 . (b) Noisy y = x 0 + w. (e) MSE of x(y) and its re-fitting R x(y) w.r.t. γ.Two values of γ are selected corresponding to (c) re-fitting for a sub-optimal γ, (d) original for a sub-optimal γ, (f) original for its optimal γ, (g) re-fitting for its optimal γ.

Fig. 13 .

 13 Fig. 13. (a) Noise-free x 0 . (b) Noisy y = x 0 + w. (c) Original aniso-TV x(y) with λ = 1.2. Enhanced results by (d) iterative hard-thresholding, (e) post re-fitting with support identification, and (f) our proposed covariant re-fitting.

20 ,

 20 SSIM: 0.785 PSNR: 25.44, SSIM: 0.794 PSNR: 25.99, SSIM: 0.809 PSNR: 26.71, SSIM: 0.793 Bregman PSNR: 27.01, SSIM: 0.828 PSNR: 24.68, SSIM: 0.718 PSNR: 22.70, SSIM: 0.638 PSNR: 22.42, SSIM: 0.627 SOS-boosting PSNR: 24.85, SSIM:

Fig. 15 .

 15 Fig. 15. (a) Cartoon image x 0 . (b) Noisy version y = x 0 +w. (c) Iso-TV x(y) with λ = 1.5. (d) Our covariant re-fitting R x(y).(e,f,g,h) From top to bottom, boosting[START_REF] Tukey | Exploratory data analysis[END_REF][START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF], Bregman iterations[START_REF] Burger | Nonlinear inverse scale space methods[END_REF] and SOS-boosting[START_REF] Romano | Boosting of image denoising algorithms[END_REF] at respectively 1, 2, 5 and 20 iterations.

CFig. 17 .

 17 Fig. 17. Comparison of re-fitting and boosting with iso-TV for varying parameter value λ in AWGN contexts with (top) standard deviation σ = 20 and (bottom) σ = 60. (a) Noisy image y = x 0 + w, (b) Original estimate x(y), (c) Our re-fitting x(y), (d) Twicing, (f) One Bregman iteration, (g) SOS-boosting with one iteration. (b-d,f-g) All of this results are obtained at their own respective optimal λ value. (e) MSE of x(y), re-fitting x(y) and the three boosting approaches.

.2. Re-fitting examples. We

  now exemplify the previous definitions for the various variational estimators introduced in Subsection 2.2.

	1 2 ||Φx -y|| 2 2 , but extensions to a general F (e.g., logistic regression) is straightforward and reads R inv x∈M x (y) x (y) ∈ argmin F (x, y).
	Remark 4. When x(y) ∈ Im[J], then M x(y) = Im[J] and R inv x (y) = J(ΦJ) + (y).
	3

Often, the solution of (2) is non unique, but for simplicity we only consider a selection of such solutions, and we assume that the selected path x : y → x(y) is differentiable almost everywhere.

This could be enforced as shown in[START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF].

If |Γx | i = 0 for all i ∈ [m], the result remains true for any α > 0.

PSNR = 10 log 10 255 2 / 1 p ||x(y) -x 0 || 2 2 for an image ranging on [0, 255].

analysis described in Equation (5), a candidate for W is (56)W = U (ΦU ) + -U (U Φ ΦU ) -1 U (Γ ) I R ,

 (e,f,g,h) From top to bottom, boosting [START_REF] Tukey | Exploratory data analysis[END_REF][START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF], Bregman iterations [START_REF] Burger | Nonlinear inverse scale space methods[END_REF] and SOS-boosting [START_REF] Romano | Boosting of image denoising algorithms[END_REF] at respectively 1, 2, 5 and 20 iterations. This is why CLEAR only considers the linearization of x at y through the Jacobian.

It is important to have in mind that most of the theoretical results regarding boosting methods are often well grounded in the case where, even though x is nonlinear, it acts locally as an averaging filter. In other words, locally, there exists a row stochastic linear operator W , i.e., W 1 n = 1 n , (or even bi-stochastic, symmetric or independent of y) such that x(y) = W y. Such theoretical results could not be applied to the soft-thresholding (nor to more advanced methods we have considered). Indeed, for the soft-thresholding, a candidate for W is the diagonal matrix defined as

which is not row stochastic. In this context, the matrix W is not the Jacobian, which is given in this case by the diagonal matrix

which, unlike W , is row stochastic and locally a projector. A second limitation is that even though W is row stochastic, it might still encode a bias part. Typically, for the Appendix A. Sketch of proofs. This section details how to retrieve closedform expressions of some of the estimators studied in the paper.

A.1. Retrieving the least-square solution. We aim at retrieving here a minimizer of ||Φx -y|| 2 2 + ι C (x) where

The initial problem can be recast as (57) argmin

The right hand side problem being differentiable and convex, its first order optimality conditions give A Φ ΦAt = A Φ (y -Φb). In particular t = (ΦA) + (y -Φb) is a solution, and hence

A.2. Retrieving the Tikhonov solution. We consider the minimization problem, defined for Γ ∈ R m×p and λ > 0, of

The objective being differentiable and convex, its first order optimality conditions give

Provided Ker Φ ∩ Ker Γ = {0} and λ > 0, the quantity Φ Φ + λΓ Γ is invertible and x = (Φ Φ + λΓ Γ) -1 Φ y is the unique solution.

A.3. Retrieving the hard-thresholding solution. We consider the minimization problem, defined for λ > 0, of

The problem is separable meaning that argmin

otherwise , which is reached by setting x i = 0 when |y i | λ and x i = y i otherwise.

A.4. Retrieving the soft-thresholding solution. We consider the minimization problem, defined for λ > 0, of E(x, y) = 1 2 ||x -y|| 2 2 + λ||x|| 1 , which is, as for the hard-thresholding, separable with

By convexity, a minimum is reached when zero belongs to its sub-differential. Hence,

which holds true by setting x i = 0 when |y i | λ and x i = y i -λ sign(y i ) otherwise.

A.5. Retrieving the (block-wise) non-local means.

We assume periodical boundary conditions such that all quantities q indexed by i = (i 1 , i 2 ) ∈ Z 2 satisfies q i = q (i1,i2) = q (i1+k1p1,i2+k2p2) for all (k 1 , k 2 ) ∈ Z 2 . This leads, for l ∈ [-b, b] 2 , to the following relationship

For all i ∈ [p 1 ] × [p 2 ], studying the first optimality conditions gives (63) ∂F (x, y) ∂x i = 0 ⇔ j wi,j (x i -y j ) = 0 ⇔ x i = j wi,j y j j wi,j .

Appendix B. Proof of Theorem 20.

Before turning to the proof of this theorem, let us introduce a first lemma.

Lemma 22. For all y, let x(y) be a solution of

with F , G two convex functions and G being 1-homogeneous. Then for all ε ∈ [0, 1], the following holds

Proof. Note that if G is a convex and 1-homogeneous function, then G is subadditive, i.e., G(a) + G(b) G(a + b). Next, assume that, for some ε ∈ [0, 1], Eq. (65) does not hold, so that there exists v such that

It follows that

F (y -εΦx(y) -Φv) + G(v) + εG(x(y)) < F (yx(y)) + G(x(y)) . (67)

We also have

), since G is 1homogeneous and sub-additive. Hence, for w = v + εx(y), we get

which contradicts x(y) ∈ argmin F (y -Φx) + G(x), and then concludes the proof.

Proof of Theorem 20. By virtue of Lemma 22 and definition of x(y), we have (1 -ε)x(y) = x(y -εΦx(y)) since x(y) is supposed to be the unique solution for all y. Now, recall that the linear Jacobian operator applied to Φx(y) is the directional derivative of x(y) in the direction Φx(y), then, for almost all y, we get

which concludes the proof.

Choose ε > 0 sufficiently small such that σα -ε(1 + σ) > β. From the convergence of the primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], the sequence (z k , x k , v k ) converges to (z , x , x ). Therefore, for k large enough, |z -