Charles-Alban Deledalle

Nicolas Papadakis

Joseph Salmon

Samuel Vaiter

Block based refitting in 12 sparse regularisation

In many linear regression problems, including ill-posed inverse problems in image restoration, the data exhibit some sparse structures that can be used to regularize the inversion. To this end, a classical path is to use 12 block based regularization. While efficient at retrieving the inherent sparsity patterns of the data -the support -the estimated solutions are known to suffer from a systematical bias. We propose a general framework for removing this artifact by refitting the solution towards the data while preserving key features of its structure such as the support. This is done through the use of refitting block penalties that only act on the support of the estimated solution. Based on an analysis of related works in the literature, we introduce a new penalty that is well suited for refitting purposes. We also present a new algorithm to obtain the refitted solution along with the original (biased) solution for any convex refitting block penalty. Experiments illustrate the good behavior of the proposed block penalty for refitting solutions of Total Variation and Total Generalized Variation models.

Introduction

We consider linear inverse problems of the form y = Φx+w, where y ∈ R p is an observed degraded image, x ∈ R n the unknown clean image, Φ : R n → R p a linear operator and w ∈ R p a noise component, typically a zero-mean white Gaussian random vector with standard deviation σ > 0.

To reduce the effect of noise and the potential ill-conditioning of Φ, we consider a regularized least squares problem with a sparse analysis regularization term based on an 12 block pernalty of the form x ∈ argmin

x∈R n 1 2 Φx -y 2 2 + λ Γx 1,2 , (1)
where λ > 0 is a regularization parameter, Γ : R n → R m×b is a linear analysis operator mapping an image over m blocks of size b, and for z ∈ R m×b

z 1,2 = m i=1 z i 2 = m i=1 b j=1 z 2 i,j 1/2 , (2)
with z i = (z i,j ) b j=1 ∈ R b . Note that the terminology sparse analysis is used here by opposition to sparse synthesis models as discussed in the seminal work of [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. The first term in [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF] is a data fidelity term enforcing x to be close to y through Φ, while the second term enforces the so-called group sparsity on x (sometimes refered to as joint sparsity, block sparsity or structured sparsity) capturing the organization of the data as encoded by Γ, see for instance [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF].

Related examples

A typical example is the Lasso (Least absolute shrinkage and selection operator) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. The Lasso is a statistical procedure used for variable selection and relying on regularized linear least square regression as expressed in eq. [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF] in which Γ = Id, m = n and b = 1. The Lasso is known to promote sparse solutions, i.e., such that xk = 0 for most indices 1 ≤ k ≤ n. Since blocks are of size b = 1, the regularization term boils down to the classical 1 sparsity term that is unstructured as no interactions between the elements of x are considered. In this paper, we will focus instead on cases of block penalties where b > 1. The group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Lin | Component selection and smoothing in multivariate nonparametric regression[END_REF] is one of them, for which Γ is designed to reorganize the elements xk into groups ẑi = (Γx) i supposedly meaningful according to some prior knowledge on the data. The group Lasso is known to promote block sparse solutions, i.e., such that ẑi = 0 b for most of the groups 1 ≤ i ≤ m. Note that elements within a non-zero group ẑi = 0 b are not required to be sparse.

Regarding image restoration applications, the authors of [START_REF] Peyré | Group sparsity with overlapping partition functions[END_REF] use an 12 regularization term where Γ extracts m = n overlapping blocks of wavelet coefficients of size b or where Γ use a dyadic decomposition of the wavelet coefficients into blocks of variable size but non-overlapping [START_REF] Peyré | Adaptive structured block sparsity via dyadic partitioning[END_REF]. Such strategy was also used in audio processing [START_REF] Yu | Audio denoising by time-frequency block thresholding[END_REF] for denoising. Another example that we will investigate here is the one of the total-variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. We can distinguish two different forms of TV models. Anisotropic total-variation (TVaniso) [START_REF] Esedoḡlu | Decomposition of images by the anisotropic rudin-osher-fatemi model[END_REF], considers Γ the operator which concatenates the vertical and horizontal components of the discrete gradients into a vector of size m = 2n, hence b = 1. Isotropic total-variation (TViso) considers instead Γ = ∇ being the operator which extracts m = n discrete image gradient vectors of size b = 2. Unlike TVaniso, TViso jointly enforces vertical and horizontal components of the gradient to be simultaneously zero. Since TVaniso does not take into account interactions between both directions, it over favors vertical and horizontal structures while TViso behaves similarly in all directions, hence their name [START_REF] Esedoḡlu | Decomposition of images by the anisotropic rudin-osher-fatemi model[END_REF]. Both models promote sparsity of the discrete gradient field of the image, and, as a result, their solutions are piece-wise constant. A major difference is that TVaniso favors constant regions that are rectangular-like shaped and separated by sharp edges, while TViso favors constant regions that are rounded-like shaped and separated by fast but gradual transitions [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF].

As TV promotes piece-wise constant solutions, it is known to produce staircasing artifacts that are all the more harmful as the images contain shaded objects [START_REF] Dobson | Recovery of blocky images from noisy and blurred data[END_REF][START_REF] Chan | High-order total variation-based image restoration[END_REF]. To reduce this effect, the authors of [START_REF] Liu | Image restoration using total variation with overlapping group sparsity[END_REF] suggested using an 12 block sparsity term not only by grouping vertical and horizontal components of the gradient, but by grouping neighboring gradients in overlapping patches. An alternative to reduce staircasing, that we will also investigate here, is the second order Total Generalized Variation (TGV) model [START_REF] Bredies | Total generalized variation[END_REF] that promotes piece-wise affine solutions. As we will see, TGV is another example of models that falls into this type of least squares problems regularized with a sparse analysis term based on 12 block penalties. Sparsity in that case encodes that sought images are composed of few shaded regions with few variations of slopes and separated by edges.

Support of the solution

Solutions of the sparse analysis regularization model in (1) are known to be sparse [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Nam | The cosparse analysis model and algorithms[END_REF], i.e., such that (Γx) i = 0 b for most blocks 1 ≤ i ≤ m. It results that a key notion, central to all of these estimators, is the one of support, i.e., the set of non-zero blocks in Γx, defined as

Î = supp(Γx) = {1 ≤ i ≤ m : (Γx) i = 0 b } . (3)
For the group Lasso, the support is typically used to identify groups of covariates (columns of Φ) being explanatory variables for the dependent variable y (i.e., significantly correlated with y).

For TViso, the support is the set of pixel indices where transitions occur in the restored image.

In general, the support plays an important role as it captures the intrinsic structural information underlying the data. While being biased, in practice, the estimate x obtained by sparse analysis regularization (1) recover quite correctly the support supp(Γx) of the underlying sparse signal x. Under some additional assumptions, support recovery is even proven to be exact as proved in [START_REF] Vaiter | Robust sparse analysis regularization[END_REF] for b = 1 (anisotropic case) and [START_REF] Vaiter | Model consistency of partly smooth regularizers[END_REF] for b ≥ 1.

Bias of the solution

Though the support Î of x can be very close to the one of the sought image x, the estimated amplitudes xi suffers from a systematical bias. When Φ = Id, the Lasso corresponds to the soft-thresholding (ST) operator

xi = ST(y i , λ) = max(y i -λ sign y i , 0) (4)
for which all non-zero elements of the solutions are shrinked towards 0 by a shift ±λ resulting to under-and over-estimated values. With TViso, this bias is reflected by a loss of contrast in the image since the amplitudes of some regions are regressed towards the mean of the image [START_REF] Strong | Edge-preserving and scale-dependent properties of total variation regularization[END_REF][START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF][START_REF] Vaiter | The degrees of freedom of partly smooth regularizers[END_REF]. In TGV, not only a loss of contrast results from this bias, but one can observe that the slopes in areas of transitions are often over estimated.

Boosting approaches

Given the artifacts induced by the 12 sparse regularization, many approaches have been developed to re-enhance the quality of the solutions, e.g., to reduce the loss of contrast and staircasing of TViso. We refer to these approaches as boosting. Most of them consist in solving (1) iteratively based on the residue Φx -y, or a related quantity, obtained during the previous iterations. Among them, the well-known Bregman iterations [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] is often considered to recover part of the loss of contrast for TViso. Other related procedures are twicing [START_REF] Tukey | Exploratory data analysis[END_REF], boosting with the 2 loss [START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[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]. While these approaches reduce the bias in the estimated amplitudes, the support Î of the original solution is not guaranteed to be preserved in the boosted solution, even though this one may correspond to the support of the sought image x.

Projection on the support

Given the key role of the support of solutions of (1), we believe that it is of main importance that a re-enhanced solution x preserves it, i.e., such that supp(Γx) ⊆ Î. For this reason, we focus on re-fitting strategies that, unlike boosting, reduce the bias while preserving the support of the original solution. In the Lasso (Γ = Id, m = n and b = 1) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], 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 statistical 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). Such approaches consists in approximating y through Φ by an image sharing the same support as x:

xsupp ∈ argmin x; supp(Γx)⊆ Î 1 2 Φx -y 2 2 , (5)
where Î = supp(Γx). While this strategy works well for blocks of size b = 1, e.g., for the Lasso or TVaniso, it suffers from an excessive increase of variance whenever b ≥ 2 (see [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF] for illustrations on TViso). This is due to the fact that solutions do not only present sharp edges, but may involve gradual transitions. To cope with this issue, additional features of x than its support must be also preserved by a refitting procedure.

Advanced refitting strategies

For the Lasso, it has been observed that a pointwise preservation of the sign of xi onto the support improves the numerical performances of the refitting [START_REF] Chzhen | On lasso refitting strategies[END_REF]. For b = 2 and TViso like models, the joint projection on the support with conservation of the direction (or orientation) of (Γx) i has been proposed in [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF]. Extension to second order regularization such as TGV [START_REF] Bredies | Total generalized variation[END_REF] are investigated in [START_REF] Burger | Convergence rates and structure of solutions of inverse problems with imperfect forward models[END_REF] in the context of partially order spaces and approximate operators Φ. In a parallel line of research, it has been proposed in [START_REF] Weiss | Contrast invariant SNR and isotonic regressions[END_REF] to respect the inclusion of the level lines of x in the refitting by solving an isotonic regression problem. All these models are constrained to respect exactly the orientation (Γx) i of the biased solution on elements of the support, i.e., when i ∈ Î.

In [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF][START_REF] Pierre | Luminance-guided chrominance denoising with debiased coupled total variation[END_REF], an alternative approach, based on the preservation of covariant information between x and y, aims only at preserving the orientation (Γx) i to some extent. While also respecting the support of x, this gives more flexibility for the refitted solution to correct x and adapt to the data content y. This model is nevertheless insensitive to the direction and it involves a quadratic penalty that tends to promote over smoothed refittings.

Outline and contributions

In Section 2, we present a general framework for refitting solutions promoted by 12 sparse regularization (1) that extends a preliminary version of this work [START_REF] Deledalle | Refitting solutions promoted by 12 sparse analysis regularizations with block penalties[END_REF]. Our variational refitting method relies on the use of block penalty functions that act on the support of the biased solution x. We introduce the Soft-penalized Direction model (SD), while discussing suitable properties a refitting block penalty should satisfy.

In Section 3, we propose stable algorithms to compute our refitting strategy for any convex refitting block penalty.

We show in Section 4 how our model relates and inherits the advantages of other methods such as Bregman iterations [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] or de-biasing approaches [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF][START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF].

Experiments in Section 5 exhibit the practical benefits for the SD refitting for imaging problems involving TViso, a variant of TViso for color images and TGV based regularization.

Refitting with block penalties

The refitting procedure of a biased solution x of (1) is expressed in the following general framework

xφ ∈ argmin x; supp(Γx)⊆ Î 1 2 Φx -y 2 2 + i∈ Î φ((Γx) i , (Γx) i ) , (6)
where φ : R b × R b → R is a block penalty (b ≥ 1 is the size of the blocks) promoting Γx to share information with Γx in some sense to be specified. To compute global optimum of the refitting model ( 6), we only consider in this paper refitting block penalties such that z → φ(z, ẑ) is convex.

To refer to some features of the vector Γx, let us first define properly the notions of relative orientation, direction and projection between two vectors. Definition 1. Let z and ẑ be vectors in R b , we define

cos(z, ẑ) = z z 2 , ẑ ẑ 2 = 1 z 2 ẑ 2 b j=1 z j ẑj , (7)
and

P ẑ (z) = z, ẑ ẑ 2 ẑ ẑ 2 = z 2 ẑ 2 cos(z, ẑ)ẑ , (8)
where P ẑ (z) is the orthogonal projection of z onto Span(ẑ) (i.e., the orientation axis of ẑ). We say that z and ẑ share the same orientation (resp. direction), if |cos(z, ẑ)| = 1 (resp. cos(z, ẑ) = 1).

We also consider that cos(z, ẑ) = 1 in case of null vectors z = 0 b and/or ẑ = 0 b .

Thanks to Definition 1, we can now introduce our refitting block penalty designed to preserve the desired features of ẑ = Γx in a simple way. We call our block penalty the Soft-penalized Direction (SD) penalty which reads as

φ SD (z, ẑ) = λ z 2 (1 -cos(z, ẑ)) . (9)
We also introduce five other alternatives, the Hard constrained Orientation (HO) penalty

φ HO (z, ẑ) = ι {z∈R b : |cos(z,ẑ)|=1} (z) , (10)
where ι C is the 0/+∞ indicator function of a set C, the Hard-constrained Direction (HD) penalty

φ HD (z, ẑ) = ι {z∈R b : cos(z,ẑ)=1} (z) , (11)
the Quadratic penalized Orientation (QO) penalty

φ QO (z, ẑ) = λ z 2 2 2 ẑ 2 (1 -cos 2 (z, ẑ)) , (12)
the Quadratic penalized Direction (QD) penalty

φ QD (z, ẑ) = λ 2 z 2 2 ẑ 2 (1 -cos 2 (z, ẑ)) if cos(z, ẑ) ≥ 0 λ 2 z 2 2 ẑ 2 otherwise ( 13
)
and the Soft-constrained Orientation (SO) penalty

φ SO (z, ẑ) = λ z 2 1 -cos 2 (z, ẑ) . ( 14
)
We will see in Section 4 that HD, HO and QO lead us to retrieve existing refitting models known respectively in the litterature as ICB (Infimal Convolution betweem Bregman distances) debiasing [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF], Bregman debiaising [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF], and CLEAR (Covariant LEAst square Refitting) [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF].

Desired properties of refitting block penalties

We now introduce properties a block penalty φ should satisfy for refitting purposes, for any ẑ:

(P1) φ is convex, non negative and φ(z, ẑ) = 0, if cos(z, ẑ) = 1 or z 2 = 0, (P2) φ(z , ẑ) ≥ φ(z , ẑ) if z 2 = z 2 and cos(z , ẑ) ≥ cos(z , ẑ), (P3) z → φ(z, ẑ) is continuous, Table 1: Properties satisfied by the considered block penalties. e √ √ √ √ (P4) φ(z, ẑ) ≤ C z 2 , for C > 0.
Property (P1) stipulates that no configuration can be more favorable than z and ẑ having the same direction. Hence, the direction of the refitted solution should be encouraged to follow the one of the biased solution. Property (P2) imposes that for a fixed amplitude, the penalty should be increasing w.r.t. the angle formed with ẑ. Property (P3) enforces refitting that can continuously adapt to the data and be robust to small perturbations. Property (P4) claims that a refitting block penalty should not penalize more some configurations than the original penalty . 1,2 , at least up to some multiplicative factor C > 0.

The next proposition, whose proof is given in Appendix A, provides additional properties.

Proposition 1. Properties of block penalties lead to the following implications.

(a) (P1) ⇒ Φx φ -y 2 2 ≤ Φx -y 2 2 . (b) (P1) ⇒ φ is non decreasing with respect to z 2 for a fixed angle (z, ẑ), (c) (P2) ⇒ φ is symmetric with respect to the orientation axis induced by ẑ,

(d) (P1) + (P2) ⇒ φ(z , ẑ) ≥ φ(z, ẑ), if z 2 ≥ z 2 and cos(z, ẑ) = cos(z , ẑ), (e) (P1) + (P3) ⇒ φ(z, ẑ) → 0 when cos(z, ẑ) → 1,

Properties of considered block penalties

The properties of the previously introduced refitting block penalties are synthesized in Table 1. The proposed SD model is the only one satisfying all the desired properties. As illustrated in Fig. 1, it is a continuous penalization that increases continuously with respect to the absolute angle between z and ẑ.

Other block penalties are insensitive to directions (HO, QO and SO), completely intolerant (HD and HO) or too tolerant (QD) to small changes of orientations, hence not satisfying. These drawbacks will be illustrated in our experiments conducted in Section 5.

When b = 1, the orientation-based penalties (QO, HO and SO) have absolutely no effect while the direction-based penalties HD and QD preserve the sign of (Γx) i . In this paper, when b ≥ 1, we argue that the direction of the block (Γx) i carries important information that is worth preserving when refitting, at least to some extent.

Refitting in practice

We now introduce a general algorithm aiming to jointly solve the original problem (1) and the refitting one [START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF] for any refitting block penalty φ. This framework has been extended from the stable projection onto the support developped in [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF] and later adapted to refitting with the Quadratic Orientation penalty in [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF]. Given x solution of (1), a posterior refitting can be obtained by solving [START_REF] Bühlmann | Boosting with the L2 loss: regression and classification[END_REF] for any refitting block penalty φ. To that end, we write the characteristic function of support preservation as i∈ Îc ι {0} (ξ i ), where ι {0} (ξ) = 0 if z = 0 and +∞ otherwise. By introducing the convex function

ω φ (ξ, Γx, Î) = i∈ Îc ι {0} (ξ i ) + i∈ Î φ(ξ i , (Γx) i ) , (15)
the general refitting problem ( 6) can be expressed as

xφ ∈ argmin x∈R n 1 2 Φx -y 2 2 + ω φ (Γx, Γx, Î) . ( 16
)
We now describe two iterative algorithms that can be used for the joint computation of x and x.

Primal-dual formulation

We first consider the primal dual formulation of the problems ( 1) and ( 16) that reads

min x∈R n max ξ∈R b 1 2 Φx -y 2 2 + Γx, ξ -ι B λ 2 (ξ) , (17)
min x∈R n max ξ∈R b 1 2 Φx -y 2 2 + Γx, ξ -ω * φ (ξ, Γx, Î) , (18)
where ι B λ 2 is the indicator function of the 2 ball of radius λ (that is 0 if z i 2 ≤ λ for all i ∈ [m] and +∞ otherwise) and

ω * φ (ξ, Γx, Î) = sup ζ∈R b ξ, ζ -ω φ (ζ, Γx, Î) (19)
is the convex conjugate, with respect to the first argument of ω φ (•, Γx, Î). Two iterative primal-dual algorithms are used to solve these problems. They involve the biased variables (ẑ k , xk ) and the refitted ones (z k , xk ). Let us now present the whole algorithm, defined for parameters κ > 0, τ > 0 and θ ∈ [0, 1] as: [START_REF] Efron | Least angle regression[END_REF] with the operator Φ - τ = (Id + τ Φ t Φ) -1 and Π, prox κω * φ and Ψ are functions that will be detailed bellow. As we will see in Proposition 2, the quantity β > 0 is used to guarantee that we estimate the support of x correctly. In practice, we choose β as the smallest available non-zero floating number. The process for the biased variable involves the orthogonal projection on the 2 ball of radius

                                 νk+1 = ẑk + κΓv k νk+1 = zk + κΓṽ k ξk+1 = Π(ν k+1 , λ) Îk+1 = i ∈ [m] : νk+1 i 2 > λ + β ξk+1 = prox κω * φ (ν k+1 , Ψ(ν k+1 ), Îk+1 ) xk+1 = Φ - τ xk + τ (Φ t y -Γ t ξk+1 ) xk+1 = Φ - τ xk + τ (Φ t y -Γ t ξk+1 ) vk+1 = xk+1 + θ(x k+1 -xk ) ṽk+1 = xk+1 + θ(x k+1 -xk ),
λ in R b Π( ξ, λ) i = λ ξi max(λ, ξi 2 ) .
wheras the refitting algorithm relies on the proximal operator of ω * φ . We recall that the proximal operator of a convex function ψ at point ξ 0 reads

prox κψ (ξ) = argmin ξ 1 2κ ξ -ξ 0 2 2 + ψ(ξ) . (21)
From the block structure of the function ω φ defined in [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF], the computation of its proximal operator may be realized pointwise. Since ι {0} (ξ) * = 0, we have

prox κω * φ (ξ 0 , Γx, Î) i = prox κφ * (ξ 0 i , (Γx) i ), if i ∈ Î , ξ 0 i , otherwise . (22)
Table 2 gives the expressions of the dual functions φ * with respect to their first variable and their related proximal operators prox κφ * for the refitting block penalties considered in this paper. All details are given in the Appendix B. Following [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], for any positive scalars τ and κ satisfying τ κ Γ t Γ 2 < 1 and θ ∈ [0, 1], the estimates ( ξk , xk , vk ) of the biased solution converge to ( ξ, x, x), where ( ξ, x) is a saddle point of [START_REF] Deledalle | Refitting solutions promoted by 12 sparse analysis regularizations with block penalties[END_REF]. When the last arguments of the function ω * φ are the converged Γx and its support Î, the refitted variables converge to a saddle point of [START_REF] Dobson | Recovery of blocky images from noisy and blurred data[END_REF]. However, we did not succeed to show convergence for the refitting process, since the quantities Γx are Î are only estimated from the biased variables at the current iteration, as explained in the next paragraphs.

Online support identification. Estimating supp(Γx) from an estimation xk is not stable numerically: the support supp(Γx k ) can be far from supp(Γx) even though xk is arbitrarily close to x. As in [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF], we rather consider the dual variable ξk to estimate the support. From relations in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], we indeed expect at convergence ξk to saturate on the support of Γx and to satisfy the optimality condition ξk i = λ (Γx)i (Γx)i 2 . In practice, the norm of the dual variable ξk i saturates to λ relatively fast onto Î. As a consequence, it is far more stable to detect the support of Γx with the dual variable ξk than with the vector Γx k itself. The next propostion, adapted from [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF] and whose sketch of proof is given in Appendix A, shows that this approach can indeed converge towards the support Î.

Proposition 2. Let α > 0 be the minimum non zero value of (Γx) i 2 , and choose β such that ακ > β > 0. Then, denoting νk+1 = ξk + κΓv k , [START_REF] Efron | Least angle regression[END_REF] converges in finite time to the true support Î = supp(Γx) of the biased solution 1 .

Îk+1 = i ∈ [m] : νk+1 i 2 > λ + β in Algorithm
Online direction and norm identification. In algorithm [START_REF] Efron | Least angle regression[END_REF], the function Ψ(ν k+1 ) aims at approximating Γx . As for the estimation of the support, instead of directly considering Γx k , we rather rely on νk+1 = ξk + κΓv k to obtain a stable estimation ẑ of the vector Γx.

Table 2 shows that for all the considered block penalties, the computation of the proximal operator of φ * involves the normalized vector ẑi / ẑi 2 on the support Î. This direction is approximated at each iteration by normalizing νk+1 i . The amplitude ẑi 2 is just required for the models QO and QD.

In Proposition 7, we will show that the model of [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF] is related to the QO block penalty. This method considers the algorithmic differentiation of the biased process to obtain a refitting algorithm. Given its good numerical performance, we leverage this approach to define the function Ψ as shown in the next proposition whose proof can be foun in Appendix A.

Proposition 3. The approach of [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF] corresponds to refit with the QO block penalty and the following convergent estimation Ψ(ν k+1 i ) of (Γx) i .

Ψ(ν k+1 i ) = νk+1 i 2-λ κ νk+1 i 2 ν k+1 i . ( 23
)
1 As in [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF], the extended support ξi 2 = λ can be tackled by testing νk+1 i ≥ λ.

Discussion. This joint-estimation considers at every iteration k different refitting functions ω * φ (., Ψ(ν k+1 ), Îk+1 ) in [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF]. Then, unless b = 1 (see [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF]), we cannot show the convergence of the refitting scheme. As in [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF], we nevertheless observe convergence and a stable behavior for this algorithm.

In addition to its better numerical stability, the running time of joint-refitting is more interesting than the posterior approach. In Algorithm [START_REF] Efron | Least angle regression[END_REF], the refitted variables at iteration k require the biased variables at the same iteration and the whole process can be realized in parallel without significantly affecting the running time of the original biased process. On the other hand, posterior refitting is necessarily sequential and the running time is doubled in general.

Douglas-Rachford formulation

An alternative to obtain solutions x of (1) and x of ( 16) is to consider the splitting TViso reformulation, as proposed in [START_REF] Combettes | A douglas-rachford splitting approach to nonsmooth convex variational signal recovery[END_REF], and given by min

x∈R n ξ∈R b 1 2 Φx -y 2 2 + ι {x,ξ;Γx=ξ} (x, ξ) + λ ξ 1,2 , (24)
min

x∈R n ξ∈R b 1 2 Φx -y 2 2 + ι {x,ξ;Γx=ξ} (x, ξ) + ω φ (ξ, Γx, Î) . ( 25
)
This problem can be solved with the Douglas-Rachford algorithm [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF]. Introducing the parameters α ∈ (0, 2) and τ > 0, the iterates read

                                       υk+1 = Γ -(2x k -μk + Γ t (2 ξk -ζk )), υk+1 = Γ -(2x k -μk + Γ t (2 ξk -ζk )), μk+1 = μk + α(υ k+1 -xk ), μk+1 = μk + α(υ k+1 -xk ), ζk+1 = ζk + α(Γυ k+1 -ξk ), ζk+1 = ζk + α(Γυ k+1 -ξk ), xk+1 = Φ - τ (μ k+1 + τ Φ t y), xk+1 = Φ - τ (μ k+1 + τ Φ t y), ξk+1 = ST( ζk+1 , τ λ), Îk+1 = i ∈ [m] : ζk+1 i 2 > τ λ + β , ξk+1 = prox τ ω φ (•,Υ( ζk+1 ), Îk+1 ) ( ζk+1 ) (26) with Γ -= (Id + Γ t Γ) -1 , the block Soft Thresholding (ST) operator ST( ζ, λ) i = 0 if ζi 2 ≤ λ, ζi -λ ζi ζi 2 otherwise .
and the proximal operator of ω φ that can be written pointwise as

prox τ ω φ (ζ 0 , Γx, Î) i = ζ 0 -τ prox φ * /τ (ζ 0 i /τ, (Γx) i ), if i ∈ Î , 0, otherwise . (27)
The estimates (x k , ξk , ζk ) of the biased solution converge to (x, ξ, ζ) = (x, Γx, (1 + τ λ)(Γx)), from the optimality conditions of problem [START_REF] Lederer | Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions[END_REF]. On the other hand, there is again no convergence guarantee for the refitted variables since the last two arguments of the function ω φ are potentially modified at each iteration.

z 0 -P ẑ (z 0 ), if z 0 , ẑ ≥ 0 z 0 , otherwise QO ẑ 2 2λ z 2 2 , if cos(z, ẑ) = 0 +∞, otherwise λ λ+κ ẑ 2 (z 0 -P ẑ (z 0 )) QD ẑ 2 2λ z 2 2 if cos θ(z, ẑ) ≤ 0 +∞ otherwise λ λ+κ ẑ 2 z 0 -P ẑ (z 0 ) if z 0 , ẑ ≥ 0 z 0 otherwise SO 0, if cos(z, ẑ) = 0 and z 2 ≤ λ +∞, otherwise λ z0-P ẑ (z0) max(λ, z0-P ẑ (z0) 2) SD 0, if z + λ ẑ ẑ 2 ≤ λ † +∞, otherwise λ z0+λ ẑ ẑ 2 max(λ, z0+λ ẑ ẑ 2 2)
ẑ ẑ 2 † : note that the condition implies that cos(z, ẑ) ≤ 0.

The support of Γx is estimated in line from auxiliary variables as Îk+1 = i ∈ [m] : ζk+1 i 2 > τ λ + β , where, as shown in [START_REF] Deledalle | Contrast re-enhancement of total-variation regularization jointly with the douglas-rachford iterations[END_REF], β has to be taken such that 0 < β < λ min i∈ Î (Γx) i 2 to have convergence in finite time of Îk+1 to the true support Î = supp(Γx). Again the quantity β > 0 is chosen in practice as being the smallest available non-zero floating number.

Considering the stable algorithmic differentiation strategy of [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF] suggested in the previous subsection, the vector Γx is approximated on the support Îk+1 at each iteration with the function

Υ( ζ) = ζ 2-λτ ζ 2 ζ.

Related re-fitting works

We now review some other related refitting methods. First we will discuss of refitting methods based on Bregman divergences, next the refitting approch developped in [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF], and then we will see how these techniques are related to the block penalites introduced in Section 2.

Bregman-based Refitting

Bregman divergence of 12 structured regularizers

In the literature [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF][START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF], Bregman divergences have proven to be well suited to measure the discrepancy between the biased solution x and its refitting x. We recall that for a convex, proper and lower semicontinuous function ψ, the associated (generalized) Bregman divergence between x and x is, for any subgradient p ∈ ∂ψ(x):

D p ψ (x, x) = ψ(x) -ψ(x) -p, x -x ≥ 0 . ( 28
)
If ψ is an absolutely 1-homogeneous function, i.e., ψ(αx) = |α|ψ(x), ∀α ∈ R, then

p ∈ ∂ψ(x) ⇒ ψ(x) = p, x , (29)
and the Bregman divergence simplifies into

D p ψ (x, x) = ψ(x) -p, x . (30)
As an example, let us consider ψ(x) = x 2 . Since ψ is 1-homogeneous, it follows that

D p • 2 (x, x) = x 2 -p, x (31)
where

p ∈ ∂ • 2 (x) ⇔ p = x x 2 if x = 0 , p 2 ≤ 1 otherwise . ( 32
)
For regularizers of the form ψ(x) = Γx 1,2 , we introduce the following notations

Ω(x) = ∂ Γ • 1,2 (x) , (33)
∆(x) = η : Γ t η ∈ ∂ Γ • 1,2 (x) . (34)
We have, see for instance [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF], for all p ∈ R n :

p ∈ Ω(x) ⇔ ∃η ∈ ∆(x) s.t. p = Γ t η ( 35
)
or in short Ω(x) = Γ t ∆(x). Interestingly ∆(x) enjoys a separability property in terms of all subgradients associated to the 2 norms of the blocks

∆(x) = {η ∈ R m : ηi ∈ ∂ • 2 ((Γx) i ), ∀i ∈ [m]} . (36)
Additionally, the next proposition shows that there is a similar separability property for the Bregman divergence (refer to Appendix A for the proof).

Proposition 4. Let η ∈ ∆(x) and p = Γ t η. Then

D p Γ• 1,2 (x, x) = m i=1 D ηi • 2 ((Γx) i , (Γx) i ) . (37)
In the following we consider x and its support Î to be fixed, and we denote by D η i (Γx) the following

D η i (Γx) = D ηi • 2 ((Γx) i , (Γx) i ) . (38)
The next proposition shows that such a divergence measures the fit of directions between (Γx) i and (Γx) i , but also partially captures the support Î (refer to Appendix A for the proof).

Proposition 5. Let η ∈ ∆(x). We have for i ∈ Î

D η i (Γx) = 0 ⇔ ∃α i ≥ 0 s.t. (Γx) i = α i (Γx) i . ( 39
)
and we have for i ∈ Îc

D η i (Γx) = 0 ⇔ (Γx) i = 0 b or ηi = (Γx) i (Γx) i 2 . ( 40
)
As noticed in the image fusion model in [START_REF] Ballester | A variational model for p+ xs image fusion[END_REF], minimizing D η i (Γx) enforces the alignment of the direction (Γx) i with ηi and it is an efficient way to avoid contrast inversion.

In the following sections, unless stated otherwise, we will always consider η ∈ ∆(x) and p = Γ t η ∈ Ω(x).

Iterative Bregman regularization

The Bregman process [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] reduces the bias of solutions of (1) by successively solving problems of the form xl+1 ∈ argmin

x∈R n 1 2 Φx -y 2 2 + λD pl Γ• 1,2 (x, xl ) . ( 41
)
with pl ∈ Ω(x l ). We consider a fixed λ, but different strategies can be considered with decreasing parameters λ l as in [START_REF] Scherzer | Inverse scale space theory for inverse problems[END_REF][START_REF] Tadmor | A multiscale image representation using hierarchical (BV,L2) decompositions[END_REF]. For l = 0, setting x0 = 0 n , and taking p0 = 0 n ∈ Ω(x 0 ) so that

D p0 Γ• 1,2 (x, x0 ) = Γx 1,2
, the first step exactly gives the biased solution of (1) with x1 = x. We denote by xIB(p) = x2 the solution obtained after 2 steps of the Iterative Bregman (IB) procedure [START_REF] Tukey | Exploratory data analysis[END_REF]:

xIB(p) = x2 ∈ argmin x∈R n 1 2 Φx -y 2 2 + λD p Γ• 1,2 (x, x) . ( 42
)
As underlined in relation [START_REF] Talebi | How to saif-ly boost denoising performance[END_REF], by minimizing

D p Γ• 1,2 (x, x) = m i=1 D η i (Γx)
, one aims at preserving the direction of Γx on the support Î, without ensuring supp(Γx IB( p) ) ⊆ Î.

For the iterative framework, the support of the previous solution may indeed not be preserved ( (Γx l ) i 2 = 0 (Γx l+1 ) i 2 = 0) and can hence grow. The support of Γx 0 for x0 = 0 n is for instance totally empty whereas the one of x = x1 may not (and should not) be empty. For l → ∞, the process actually converges to some x such that Φx = y. Because the IB procedure does not preserve the support of the solution, it cannot be considered as a refitting procedure and is more related to boosting approaches as discussed in Section 1.

Bregman based refitting

In order to respect the support of the biased solution x and to keep track of the direction Γx during the refitting, the authors of [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF] proposed the following model:

xB(p) ∈ argmin x; p∈Ω(x) 1 2 Φx -y 2 2 . ( 43
)
This model enforces the Bregman divergence to be 0, since, from eq. ( 29), we have:

p ∈ Ω(x) ⇒ D p Γ• 1,2 (x, x) = 0 . ( 44
)
We see from [START_REF] Talebi | How to saif-ly boost denoising performance[END_REF] that for i ∈ Î, the direction of (Γx) i is preserved in the refitted solution. From [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], we also observe that the absence of support is also preserved for any i ∈ Îc where ηi 2 < 1. Note that extra elements in the support Γx B( p) may be added at coordinates i ∈ Îc where ηi 2 = 1.

Infimal Convolutions of Bregman (ICB) distances based refitting

To get rid of the direction dependency, the ICB (Infimal Convolutions of Bregman distances) model is also proposed in [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF]:

xICB(p) ∈ argmin x;± p∈Ω(x) 1 2 Φx -y 2 2 . ( 45
)
The orientation model may nevertheless involve contrast inversions between biased and refitted solutions. In practice, relaxations are used in [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF] by solving, for a large value γ > 0,

xB(p) γ ∈ argmin x∈R n 1 2 Φx -y 2 2 + γD p Γ• 1,2 (x, x) . ( 46
)
The main advantage of this refitting strategy is that no support identification is required since everything is implicitly encoded in the subgradient p. This makes the process stable even if the estimation of x is not highly accurate. The support of Γx is nevertheless only approximately preserved, since the constraint D p Γ• 1,2 (x, x) = 0 can never be ensured numerically with a finite value of γ.

The best of both Bregman worlds

From relations [START_REF] Strong | Edge-preserving and scale-dependent properties of total variation regularization[END_REF] and [START_REF] Tadmor | A multiscale image representation using hierarchical (BV,L2) decompositions[END_REF], the refitting models given in ( 42) and ( 43) can be reexpressed as a function of η xIB(η) ∈ argmin

x∈R n 1 2 Φx -y 2 2 + λ m i=1 D η i (Γx) , (47)
xB(η) ∈ argmin

x∈R n 1 2 Φx -y 2 2 s.t. D η i (Γx) = 0, ∀i ∈ [m] . ( 48
)
Alternatively we can defined a mixed model, that we coin Best of Both Bregman (BBB), as

xBBB(η) ∈ argmin x∈R n 1 2 Φx -y 2 2 + λ i∈ Î D η i (Γx) , (49)
s.t. D η i (Γx) = 0, ∀i ∈ Îc . (50)
With such reformulations, connections between refitting models ( 47), ( 48) and ( 49) can be clarified. The solution xIB [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] is too relaxed, as it only penalizes the directions (Γx) i using (Γx) i , without aiming at preserving the support of x. The solution xB [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF] is too constrained: the direction within the support is required to be preserved exactly. The proposed refitting xBBB lies in-between: it preserves the support, while authorizing some directional flexibility, as illustrated by the sharper square edges in Figure 1(g). An important difference with BBB is that we consider local inclusions of subgradients of the function λ • 1,2 at point (Γx) i instead of the global inclusion of subgradients of the function λ Γ • 1,2 at point x as in [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF] and [START_REF] Vaiter | Model consistency of partly smooth regularizers[END_REF]. Such a change of paradigm allows to adapt the refitting locally by preserving the support while including the flexibility of the original Bregman approach [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF].

Covariant LEAst Square Refitting (CLEAR)

We now describe an alternative way for performing variational refitting. When specialized to 1,2 sparse analysis regularization, CLEAR, a general refitting framework [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF], consists in computing

xCLEAR ∈ argmin x; supp(Γx)⊆ Î 1 2 Φx -y 2 2 + i∈ Î λ 2 (Γx)i 2 (Γx) i -P (Γx)i ((Γx) i ) 2 , (51)
where we recall that P (Γx)i (.) is the orthogonal projection onto Span(Γx) i . This model promotes refitted solutions preserving to some extent the orientation Γx of the biased solution. It also shrinks the amplitude of Γx all the more that the amplitude of Γx are small. This penalty does not promote any kind of direction preservation, and as for the ICB model, contrast inversions may be observed between biased and refitted solutions. The quadratic term also over-penalizes large changes of orientation.

Equivalence

In the next proposition we show that, for a given vector η, the concurrent Bregman based refitting techniques coincides with two of the refitting block penalties introduced in Section 2, namely HD, HO, and that SD is nothing else than the proposed Best of Both Bregman based penalty.

Proposition 6. Let η ∈ R m be defined as

η = arg min η∈∆(x) η 2 . (52)
Then, for any y, we have the following equalities

1. xB(η) = xHD , 2. xICB(η) = xHO , 3. xBBB(η) = xSD .
where the equalities have to be understood as an equality of the corresponding sets of minimizers.

The proof of Proposition 6 can be found in Appendix A. Additionaly, we can show that CLEAR corresponds to the QO refitting block penalty also introduced in Section 2.

Proposition 7. For any y, we have xCLEAR = xQO where the equality has to be understood as an equality of the sets of minimizers.

We invite the reader to refer to Appendix A for a proof of Proposition 7.

Experiments and Results

Toy experiments with TViso

We first consider TViso regularization of grayscale images on denoising problems of the form y = x + w where w is an additive white Gaussian noise with standard deviation σ. We start with a simple toy example of a 128 × 128 image2 composed of elementary geometric shapes, and we chose σ = 50. The corresponding image and its noisy version are given in Figs. 2.(a) and (b) respectively. To highlight the different behaviors between the six refitting block penalties, we chose to set the regularisation parameter of TViso to a large value λ = 750 leading to a strong bias in the solution. Fig. 2.(c) depicts this solution in which not only some structures are lost (the darkest disk), but a large loss of contrast can be observed (the white large square became gray). Unlike boosting, the purpose of refitting is not to recover the lost structures, but only to recover the correct amplitudes of the reconstructed objects without changing their geometrical aspect. In a second scenario, we consider the standard cameraman image of size 256 × 256, and we set σ = 20 and λ = 36. The corresponding images are given in Figs. 3.(a-c).

We applied the iterative primal-dual algorithm with our joint-refitting (Algorithm (20)) for 4, 000 iterations, with τ = 1/ ∇ 2 , κ = 1/ ∇ 2 and θ = 1, and where ∇ 2 = 2 √ 2. Results of refitting with HO, HD, QO, QD, SO and SD are given for the two scenarios on Figs. 2.(d-i) and Figs. 3.(d-i), respectively. As a quantitave measure of refitting between an estimate x and the underlying image x, we consider the Peak Signal to Noise Ratio (PSNR) defined in this case as PSNR = -10 log 10 1

n x -x 2 2 . (53)
The higher the PSNR, the better the refitting is supposed to be. In Fig. 2, we first observe that block penalties that are only based on orientation (HO, QO and SO) lead to refitted solution with contrast inversions compared to the original solution of TViso (see the banana dark shapes at the interface between the two squares). In that scenario, HD, QD and SD provides statisfying visual quality, the lattest achieving highest PSNR. In the context of Fig. 3, the image being more complex, the TViso solution presents more level lines than in the previous scenario. The block penalties HO and HD must preserve exactly these level lines. The refitted problem becomes too constrained and this lack of flexibility prevents HO and HD to deviate much from TViso. Their refitted solution remains significantly biased. Comparing now the quadratic penalties (QO and QD) against the soft ones (SO and SD) reveals that quadratic penalties does not allow to reenhance as much the contrast of some objects (see, e.g., the second tower on the right of the image). Again SD achieved the highest PSNR in that second scenario.

Experiments with color TViso

We now consider TViso regularization of degraded color images with three channels Red, Green, Blue (RGB). We defined blocks obtained by applying Γ

= (∇ 1 R , ∇ 2 R , ∇ 1 G , ∇ 2 G , ∇ 1 B , ∇ 2 B )
, where We first focused on a denoising problem y = x + w where x is a color image and w is an additive white Gaussian noise with standard deviation σ = 20. We next focused on a deblurring problem y = Φx + w where x is a color image, Φ is a convolution simulating a directional blur, and w is an additive white Gaussian noise with standard deviation σ = 2. In both cases, we chose λ = 4.3σ. We apply the iterative primal-dual algorithm with our joint-refitting (Algorithm (20)) for 1, 000 iterations, with θ = 1, τ = 1/ Γ 2 and κ = 1/ Γ 2 , where, as for the grayscale case, Γ 2 = 2 √ 2. Note that in order to use the primal-dual algorithm, we have to implement Γ t , defined for z = (z R , z G , z B ) ∈ R n×6 where z R , z G , z B are 2d vector fields, as given by

Γ t z = -div z R -div z G -div z B . (54)
Results are provided on Fig. 4 and5. Comparisons of refitting with our proposed SD block penalty, HO, HD, QO, QD and SO are provided. Using our proposed SD block penalty offers the best refitting performances in terms of both visual and quantitative measures. The loss of contrast of TViso is well-corrected, amplitudes are enhanced while smoothness and sharpness of TViso is preserved. Meanwhile, the approach does not create artifacts, invert contrasts, or reintroduce information that were not recovered by TViso.

Experiments with second order TGV

We finaly consider the second order Total Generalized Variation (TGV) regularization [START_REF] Bredies | Total generalized variation[END_REF] of grayscale images that can be expressed, for an image x ∈ R n and regularization parameters λ > 0 and ζ ≥ 0, as

x ∈ argmin x∈R n 1 2 Φx -y 2 2 + λ min z∈R n×2 ∇x -ζz 1,2 + E(z) 1,F (55)
where ∇ = (∇ 1 , ∇ 2 ) : R n → R n×2 , ∇ d denotes the forward discrete gradient in the direction d ∈ {1, 2}, and E : R n×2 → R n×2×2 is a symmetric tensor field operator defined for a 2d vector field z = (z 1 , z 2 ) as

E(z) = ∇1 z 1 1 2 ( ∇1 z 2 + ∇2 z 1 ) 1 2 ( ∇1 z 2 + ∇2 z 1 ) ∇2 z 2 (56)
where ∇d denotes the backward discrete gradient in the direction d ∈ {1, 2}, and • 1,F is the pointwise sum of the Frobenius norm of all matrices of a field. One can observe that for ζ = 0 this model is equivalent to TViso. Interestingly the solution of the second order TGV can be obtained by solving a regularized least square problem with an 12 sparse analysis term as

X ∈ argmin X∈R n×3 1 2 ΞX -y 2 2 + λ ΓX 1,2 (57)
where for an image x ∈ R n and a vector field z ∈ R n×2 , we consider X = x, z ∈ R n×3 and Ξ : (x, z) → Φx and we define Γ :

R n×3 → R 2n×3 as Γ(x, z) = ∇ 1 x -ζz 1 ∇ 2 x -ζz 2 0 ∇1 z 1 ∇2 z 2 1 √ 2 ( ∇1 z 2 + ∇2 z 1 ) . ( 58
)
After solving (57), the solution of (55) is obtained as x = Ξ X. We use Chambolle Pock algorithm for which we also need to implement the adjoint Γ t : R 2n×3 → R n×3 of Γ given for fields z ∈ R n×2 and e ∈ R n×3 by

Γ t z 1 , z 2 , • e 2 , e 2 , e 3 = -   ( ∇1 z 1 + ∇2 z 2 ) t (ζz 1 + ∇ 1 e 1 + 1 √ 2 ∇ 2 e 3 ) t (ζz 2 + ∇ 1 e 2 + 1 √ 2 ∇ 2 e 3 ) t   t . (59)
We focused on a denoising problem y = x + w where x is a simulated elevation profile, ranging in [0, 255], and w is an additive white Gaussian noise with standard deviation σ = 2. We chose λ = 15 and ζ = 0.45. We applied the iterative primal-dual algorithm with our joint-refitting Refitting comparisons between our proposed SD block penalty, and the five other alternatives are provided on Fig. 6. Using the SD model offers the best refitting performances in terms of both visual and quantitative measures. The bias of TGV is well-corrected, elevations (see the chimneys) and slopes (see the sidewalks) are enhanced while smoothness and sharpness of TGV is preserved. Meanwhile, the approach does not create artifacts, invert contrasts, or reintroduce information that was not recovered by TGV.

Conclusion

We have presented a block penalty formulation for the refitting of solutions obtained with 12 sparse regularization. Through this framework, desirable properties of refitted solutions can easily be promoted. Different block penalty functions have been proposed and their properties discussed. We have namely introduced the SD block penalty that interpolates between Bregman iterations [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] and direction preservation [START_REF] Brinkmann | Bias-reduction in variational regularization[END_REF].

Based on standard optimization schemes, we have also defined stable numerical strategies to jointly estimate a biased solution and its refitting. In order to take advantage of our efficient jointrefitting algorithm, we underline that it is important to consider simple block penalty functions, which proximal operator can be computed explicitly. This is the case for all the presented block penalties, more complex ones having been discarded from this paper for this reason.

Initially designed for TViso regularization, the approach has finally been extended to a generalized TGV model. Experiments show how the block penalties are able to preserve different structures of the biased solutions, while recovering the correct signal amplitude.

For the future, a challenging problem would be to define similar simple numerical schemes able to consider global properties on the refitting, such as the preservation of the tree of shapes [START_REF] Weiss | Contrast invariant SNR and isotonic regressions[END_REF] of the biased solution.

A Proof of Propositions Proposition 1. • (a) As xφ is solution of ( 6), the relation is obtained by observing that φ(ẑ, ẑ) = 0 with (P1).

• (b) Looking at a ray [0, z) for any vector z, this is a consequence of the convexity of φ and the fact that φ(z, ẑ) = 0 for z 2 = 0.

• (c) Since cos(z, ẑ) = cos(-z, ẑ), one can combine (z , z ) = (z, -z) and (z , z ) = (-z, z) in (P2) to obtain (c).

• (d) Direct consequence of points (b) and (c).

• (e) This is due to the continuity of the function φ(., ẑ) that is 0 on the ray [0, ẑ) from (P1). Proposition 2. We just give a sketch of the proof. More details can be found in [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF]. On the support, one has ( ξk

i , xk i , vk i ) → (λΓx k i / (Γx k ) i 2 , xk i , xk i ).
Then for k sufficiently large, ξk i + κ(Γv k ) i ≤ λ + β if and only if i ∈ I c . Proposition 3. The method in [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF] realizes the refitting through an algorithmic differentiation of the projection of the biased process. This leads to the Algorithm 20 except for the update of the dual variable ξi on the support i ∈ Îk+1 that reads in [START_REF] Deledalle | CLEAR: Covariant least-square refitting with applications to image restoration[END_REF]:

ξk+1 i = λ νk+1 i 2 νk+1 i -P νk+1 i (ν k+1 i ) , instead of ξk+1 i = prox κφ * QO (•,Ψ(ν k+1 i ), Îk+1 ) (ν k+1 i )
in Algorithm [START_REF] Efron | Least angle regression[END_REF]. Using the proximal operator given by the QO block penalty in Table 2 gives : νk+1

i 2 = λ + κ Ψ(ν k+1 i
) 2 which leads to the function [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF] for estimating ẑk i . We deduce that Ψ(ν k+1 i ) = Ψ( ξk +κΓv k ) → (Γx) i from the convergence of the biased variables ( ξk i , xk i , vk i ) → (λΓx k i / (Γx k ) i 2 , xk i , xk i ). Proposition 4. The proof is straigthforward and reads as follow

D p Γ• 1,2 (x, x) = Γx 1,2 -Γ t η, x (60)
where we used that if α 0 > 0, taking β = 0 and letting α → ∞ leads to φ * QD (α 0 , β 0 ) = +∞, otherwise, we used the optimality conditions on β and α giving us β = β0 ẑ 2 λ and α = α0 ẑ 2 λ .

• [SO] The convex conjugate φ * SO (α 0 , β 0 ) is sup

α,β α 0 α + β 0 , β -λ β 2 = ι {0}×B λ 2 (α 0 , β 0 ) , (75)

Figure 1 :

 1 Figure 1: Illustration of continuous block penalties QO and SD: (left) 2D level lines of φ for z = (z 1 , z 2 ) = A cos θ -sin θ sin θ cos θ ẑ, (middle) evolution regarding θ and (right) A.

Figure 2 :

 2 Figure 2: Comparison of alternative refitting approaches with the proposed SD model on a synthetic image.

Figure 3 :

 3 Figure 3: Comparison of alternative refitting approaches with the proposed SD model on the Cameraman image.

Figure 4 :

 4 Figure 4: (a) A color image. (b) A corrupted version by Gaussian noise with standard deviation σ = 20. (b) Solution of TViso. Debiased solution with (d) HO, (e) HD, (f) QO, (g) QD, (h) SO and (i) SD. The Peak Signal to Noise Ratio (PSNR) is indicated in brackets bellow each image.

Figure 5 :

 5 Figure 5: (a) A color image. (b) A corrupted version by a directional blur and Gaussian noise with standard deviation σ = 2. (c) Solution of TViso. Debiased solution with (d) HD, (e) QO and (f) SD. The PSNR is indicated in brackets bellow each image.

Figure 6 :

 6 Figure 6: (a) An elevation profile (range 0 to 255). (b) A corrupted version by Gaussian noise with standard deviation σ = 0.5. (c) Solution of TGV. Debiased solution with (d) HO, (e) HD, (f) QO, (g) QD, (h) SO and (i) SD. The PSNR is indicated in brackets bellow each profile.

= Γx 1 , 2 -( 2 ( 2 α 2

 12222 Γx) i 2 -ηi , (Γx) i D ηi • (Γx)i,(Γx)i).(62) also observe for instance that |cos(z, ẑ)| = 1 ⇔ β = 0 b-1 . All the block penalties φ(z, ẑ) can thus be expressed as φ(α, β). The convex conjugate readsφ * (α 0 , β 0 ) = sup α,β α 0 α + β 0 , β -φ(α, β) . (67) • [HO] It results that the convex conjugate φ * HO (α 0 , β 0 ) is sup α,β α 0 α + β 0 , β -ι {0 b-1 } (β) = ι {0} (α 0 ) .(68)• [HD] The convex conjugate φ * HD (α 0 , β 0 ) is sup α,β α 0 α + β 0 , β -ι R * + ×{0 b-1 } (α, β) = ι R-(α 0 ) .(69)•[QO] The convex conjugate φ * QO (α 0 , β 0 ) is sup α,β α 0 α + β 0 , β -λ condition on β gives us β = β0 ẑ 2 λ .• [QD] The convex conjugate φ * QD (α 0 , β 0 ) is sup α,β α 0 α + β 0 , β -

since the optimality condition on β gives us β 0 2 ≤ 2 + 2 + 2 +•λ 2 (

 22222 λ if β = 0, andβ 0 = λ β β 2 otherwise. • [SD] The convex conjugate φ * SD (α 0 , β 0 ) is sup α,β α 0 α + β 0 , βα 0 + λ/2) 2 ≤ λ/2 , +∞ otherwise ,(78)since ifβ 0 2 2 + (α 0 + λ/2) 2 > λ/2, then letting α → sign (α 0 + λ/2) × ∞ and β → sign β 0 × ∞ leads to φ * SD (α 0 , β 0 ) = +∞.B.2 Computing prox κφ *We here give the computation of the proximal operator prox κφ * (α 0 , β 0 ) of the different φ * that is given at point (α 0 , β 0 ) by argmin φ * (α, β).(79)• [HO] The proximal operator prox κφ * HO (α 0 , β 0 ) is given byargmin ι {0} (α) = (0, β 0 ).(80)• [HD] The proximal operator prox κφ * HD (α 0 , β 0 ) is given byargmin ι R-(α) = (min(0, α 0 ), β 0 ).(81)• [QO] The proximal operator prox κφ * QO (α 0 , β 0 ) is given by argmin condition gives λ(β -β 0 ) + κ ẑ 2 β = 0.• [QD] The proximal operator prox κφ * QD (α 0 , β 0 ) is given by argmin ẑ 2 (min(0, α 0 ), β 0 ) .(84)• [SO] The proximal operator prox κφ * QD (α 0 , β 0 ) is given by argmin [SD] The proximal operator prox κφ * QD (α 0 , β 0 ) is given byargmin α 0 + λ/2, β 0 ) max(λ/2, β 0 2 2 + (α 0 + λ/2) 2 ) -(λ/2, 0) ,(87)which just corresponds to the projection of the 2 ball of R b of radius λ/2 and center (-λ/2, 0).

Table 2 :

 2 Convex conjugates and proximal operators of the studied block penalties φ.

	φ		φ * (z, ẑ)	prox κφ * (z 0 , ẑ)
	HO	0, +∞,	if cos(z, ẑ) = 0 otherwise	z 0 -P ẑ (z 0 )
	HD	0, +∞,	if cos(z, ẑ) ≤ 0 otherwise

In this paper, we always consider images whose values are in the range [0, 255].

Acknowledgements This work was supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 777826.

 Proposition 5.• For i ∈ Î, i.e. (Γx) i 2 > 0, we get ηi = (Γx)i (Γx)i 2 from [START_REF] Scherzer | Inverse scale space theory for inverse problems[END_REF]. We conlude from the definition D η i (Γx) = (Γx) i 2 -ηi , (Γx) i = 0. • For i ∈ Îc and using [START_REF] Scherzer | Inverse scale space theory for inverse problems[END_REF] we distinguish two cases from [START_REF] Peyré | Adaptive structured block sparsity via dyadic partitioning[END_REF]. If ηi 2 < 1, then D η i (Γx) = 0 iff (Γx) i = 0 b . Otherwise one necessarily gets ηi = (Γx)i (Γx)i 2 .

Proposition 6. First notice that from [START_REF] Peyré | Adaptive structured block sparsity via dyadic partitioning[END_REF], η can be written explicitly as ηi =

• [1.] This is a direct consequence of relation [START_REF] Vaiter | Robust sparse analysis regularization[END_REF] and Proposition 5. Using (39), we have for

This is also valid for potential vanishing components (Γx) i with the considered convention cos(0 b , ẑ) = 0. We thus recover the penalty function HD in [START_REF] Charest | On iterative regularization and its application[END_REF]. Next, for all i ∈ Îc , we have ηi = 0 by assumption, hence according to eq. ( 40):

It follows that xB(η) = xHD . This case Îc is the same for all next points.

] With the ICB model, we have ±p ∈ Ω(x). For i ∈ Î, it gives with (39) that

corresponds to the penalty function HO in [START_REF] Chan | High-order total variation-based image restoration[END_REF].

• [3.] For all i ∈ Î, we have

that gives the penalty function SD in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. We then get that xBBB(η) = xSD .

Proposition 7. The case Îc is treated in Proposition 6. For i ∈ Î with the CLEAR model, we just have to observe that (Γx) i -

, which corresponds to the penalty function QO in [START_REF] Chzhen | On lasso refitting strategies[END_REF].

B Proximity operators of block penalties B.1 Convex conjugates φ *

We here compute the convex conjugate φ * of the different block penalties φ(z, ẑ) that only depends on z ∈ R b and where ẑ ∈ R b is a given fixed non null vector. We consider the following representation of the vectors z with respect to the ẑ axis: z = α ẑ ẑ 2 + β ẑ⊥ ẑ 2 . This expression is valid for the case b = 2. If b = 1 then z is only parameterized by α. When b > 2, ẑ⊥ must be understood as a subspace S of dimension b-1 and β as a vector of b-1 components corresponding to each dimension of S. With this change of variables, we have z 2 2 = α 2 + β 2 2 (since β is of dimension b -1), cos(z, ẑ) = α/ α 2 + β 2 2 , P ẑ (z) = αẑ/ ẑ 2 and z -P ẑ (z) = β ẑ⊥ / ẑ 2 . We