CLEAR: Covariant LEAst-square Re-fitting
with applications to image restoration

arXiv:1606.05158v2 [math.ST] 14 Sep 2016

Charles-Alban Deledalle∗ , Nicolas Papadakis∗ , Joseph Salmon† , AND Samuel Vaiter‡
Abstract. In this paper, we propose a new framework to remove parts of the systematic errors
affecting popular restoration algorithms, with a special focus for image processing tasks. Generalizing
ideas that emerged for `1 regularization, we develop an approach re-fitting the results of standard
methods towards the input data. Total variation regularizations and non-local means are special cases
of interest. We identify important covariant information that should be preserved by the re-fitting
method, and emphasize the importance of preserving the Jacobian (w.r.t. the observed signal) of the
original estimator. Then, we provide an approach that has a “twicing” flavor and allows re-fitting
the restored signal by adding back a local affine transformation of the residual term. We illustrate
the benefits of our method on numerical simulations for image restoration tasks.
Key words. inverse problems, image restoration, variational methods, re-fitting, twicing, boosting, debiasing
AMS subject classifications. 49N45, 65K10, 68U10.

1. 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 [44] or the total variation [38].
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 [19], Lasso-Gauss [36], OLS post-Lasso [1], Debiased Lasso (see [28, 1] 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
∗ IMB, CNRS, Université de Bordeaux, Bordeaux INP, F-33405 Talence, France ({charlesalban.deledalle, nicolas.papadakis}@math.u-bordeaux.fr).
† LTCI, CNRS, Télécom ParisTech, Université Paris-Saclay, 75013 Paris, France
(joseph.salmon@telecom-paristech.fr).
‡ IMB, CNRS, Université de Bourgogne, 21078 Dijon, France (samuel.vaiter@u-bourgogne.fr).

1

2

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

the classical post re-fitting solution described earlier. For instance, our methodology
successfully applies to the Tikhonov regularization [46], the isotropic total-variation
(iso-TV), the non-local means [3], the block matching 3D (BM3D) [10] and the Dual
Domain Image Denoising (DDID) [26]. In common variational contexts, e.g., `1 − `2
analysis [20] (encompassing the Lasso, the group Lasso [29, 54, 32], 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 [14] 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 [14] –
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 [47] and have recently been thoroughly
investigated: boosting [4], Bregman iterations and nonlinear inverse scale spaces
[33, 5, 52, 34], ideal spectral filtering in the analysis sense [23], SAIF-boosting [30, 43]
and SOS-boosting [37] 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.

CLEAR: Covariant LEAst-square Re-fitting

3

Interestingly, we have also realized that our scheme presents some similarities
with the classical shrinking estimators introduced in [41], especially as presented in
[22]. 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 [22, 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.
2. Background models and notation. We tackle the problem of estimating
an unknown vector x0 ∈ Rp from noisy and incomplete observations
(1)

y = Φx0 + w ,

where Φ is a linear operator from Rp to Rn and w ∈ Rn 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 x0 from y. We denote by x̂ : y 7→ x̂(y)
the procedure that provides such an estimation of x0 . 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 F (x, y) + G(x) .
x∈Rp

Another kind of estimator can be defined as the output of an iterative algorithm
(k, y) 7→ xk , e.g., a numerical solver. For a chosen iteration k, we define the final
estimator x̂(y) = xk . 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.
2.1. 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, vI ∈ R|I| is the sub-vector
whose elements are indexed by I ⊂ N and |I| is its cardinality. For any matrix M ,
MI 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 ∈ Rd , 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] : xi 6= 0} |.
1 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 7→ x̂(y) is differentiable almost everywhere.

4

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

y

200

Value

150

x̂(y)

100

50

Rx̂ (y)

Original: x0
Observation: y = x0 + w
TV: x̂(y)
Re-fitting: Rx̂ (y)

0

M x̂ (y)

-50
0

50

100

150

200

Position

(a)

(b)

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 Rinv
x̂ (y) as the orthogonal projection of y on Mx̂ (y).

2.2. 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 x0 belongs to the affine subspace
C = b + Im[A] with b ∈ Rp and A ∈ Rp×n , are a particular case of (2) where
F (x, y) = 21 ||y − Φx||22 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).
The Tikhonov regularization [46] (or Ridge regression [25]), used when ||Γx0 ||2
is small with Γ ∈ Rm×p , is an instance of (2) where F (x, y) = 21 ||y − Φx||22 and
G(x) = λ2 ||Γx||22 , 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 hard-thresholding [16], used when Φ = Id and x0 is sparse, is a solution
2
of (2) where F (x, y) = 21 ||y − x||22 and G(x) = λ2 ||x||0 , for some parameter λ > 0.
The hard-thresholding operation writes: x̂(y)I = yI and x̂(y)I c = 0, where I =
{i ∈ [p] : |yi | > λ} (see Appendix A.3).

The soft-thresholding [16], used when Φ = Id and x0 is sparse, is a solution of
(2) where F (x, y) = 21 ||y − x||22 and G(x) = λ||x||1 for some parameter λ > 0. The
soft-thresholding operation writes: x̂(y)I = yI − λ sign(yI ) and x̂(y)I c = 0, where I
is defined as for the hard-thresholding (see Appendix A.4).
The `1 synthesis (or Lasso [44, 16]), used when x0 is sparse, is a solution of (2)
where F (x, y) = 12 ||Φy − x||22 and G(x) = λ||x||1 for some parameter λ > 0. Provided
the solution is unique (see for instance [45]), it reads
(3)

x̂(y)I = (ΦI )+ yI − λ((ΦI )> ΦI )−1 sI

and x̂(y)I c = 0 ,

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

The `1 analysis, used when Γx0 is sparse with Γ ∈ Rm×p , is a solution of (2) where

(4)

F (x, y) =

1
||Φx − y||22
2

and G(x) = λ||Γx||1 ,

for some λ > 0 .

CLEAR: Covariant LEAst-square Re-fitting

5

Provided Ker Φ ∩ Ker Γ = {0}, there exists a solution given implicitly, see [48], as

(5)

x̂(y) = U (ΦU )+ y − λU (U > Φ> ΦU )−1 U > (Γ> )I sI ,

for almost all y and where I = supp(Γx̂(y)) = {i ∈ [m] : (Γx̂(y))i 6= 0} is called the
Γ-support of the solution, sI = 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) [38] is an instance of (4) where x0 ∈
p
p×b
Rp is identified to a b-dimensional signal,
is the
Ppfor which Γ = ∇ : R → R
discrete gradient operator and ||∇x||1 = 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 [29, 54, 32], used when Γx0 is block sparse with Γ ∈ Rp 7→ Rm×b ,
is a solution of (2) with
(6)

F (x, y) =

1
||Φx − y||22
2

and G(x) = λ||Γx||1,2 ,

Pm
for some parameter λ > 0 and where ||Γx||1,2 = i=1 ||(Γx)i ||2 . The isotropic TotalVariation (iso-TV) [38] is a particular instance of (6) where x0 ∈ Rp can
Ppbe identified
to a b-dimensional signal, for which Γ = ∇ : Rp → Rp×b and ||∇x||1,2 = 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., [6]), 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.
The (block-wise) non-local means [3], used when the image x0 ∈ Rp (with p =
p1P
× p2 ) is composed of many redundant patterns, is the solution of (2) for F (x, y) =
1
2
i,j wi,j ||Pi x − Pj y||2 and G(x) = 0, where Pi is the linear operator extracting a
2
patch (i.e., a small window) at pixel i of size (2b + 1)2 . The index i ∈ [p1 ] × [p2 ] spans
the image domain and j − i ∈ [−s, s]2 spans
 a limited search window. The weights
wi,j are defined as wi,j = ϕ ||Pi y − Pj y||22 where the kernel ϕ : R+ → [0, 1] is a
decreasing function. Assuming periodical boundary conditons, its solution is given as
P
X
j w̄i,j yj
with w̄i,j =
wi−k,j−k ,
(7)
x̂(y)i = P
j w̄i,j
k

2

where k ∈ [−d, d] spans the patch domain (see Appendix A.5).

2.3. 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 sI whose expression simplifies to λ sign(yI ) 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 sI . 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 [42]
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

6

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

(a)

(b)

(c)

(d)

(e)

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.

though TV has perfectly retrieved the support of ∇x0 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
(8)

argmin

x ; supp(x) ⊆ I

1
2 ||Φx

− y||22 .

In this section, we present a re-fitting procedure, discussed in [14], 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 7→ x̂(y)
from Rn to Rp 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.

7

CLEAR: Covariant LEAst-square Re-fitting

Definition 1. The model subspace associated to an a.e. differentiable estimator x̂
is defined at almost all points y ∈ Rn by the affine subspace of Rp
(9)

Mx̂ (y) = x̂(y) + Im [Jx̂ (y)] ,

where Jx̂ (y) ∈ Rp×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 7→ x̂(y) is given for almost all y ∈ Rn by
(10)

+
Rinv
x̂ (y) = x̂(y) + J(ΦJ) (y − Φx̂(y)) ∈ argmin

x∈Mx̂ (y)

1
2 ||Φx

− y||22 ,

where J = Jx̂ (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) = 21 ||Φx − y||22 , but extensions to a general
F (e.g., logistic regression) is straightforward and reads Rinv
x̂ (y) ∈ argmin F (x, y).
x∈Mx̂ (y)

Remark 4. When x̂(y) ∈ Im[J], then Mx̂ (y) = Im[J] and

Rinv
x̂ (y)

= J(ΦJ)+ (y).

3.2. Re-fitting examples. We now exemplify the previous definitions for the
various variational estimators introduced in Subsection 2.2.
The affine constrained least-squares have everywhere the same Jacobian matrix
>
J = A(ΦA)+ and its affine model subspace is Mx̂ (y) = b+Im[A(ΦA) ] (as Im[M + ] =
>
p
Im[M ]). Taking C = R , n = p, A = Id and b = 0, leads to x̂(y) = Φ+ y whose
model subspace is Mx̂ (y) = Im[Φ> ], reducing to Rp when Φ has full column rank. In
this case, the invariant re-fitting is Rinv
x̂ (y) = x̂(y).
The Tikhonov regularization has everywhere the same Jacobian matrix J =
(Φ> Φ + λΓ> Γ)−1 Φ> and everywhere the same affine model subspace Mx̂ (y) = Im[J].
+
If follows that Rinv
x̂ (y) = J(ΦJ) y. In particular, when Φ has full column rank,
p
inv
+
Mx̂ (y) = R and Rx̂ (y) = Φ y.

The soft- and hard-thresholding share a.e. the same Jacobian matrix given by
J = IdI ∈ Rp×|I| . Their model subspace reads as Mx̂ (y) = Im[IdI ] = Im[J] and the
invariant re-fitting is for both the hard-thresholding.
The `1 synthesis (or Lasso) has for Jacobian matrix J = IdI (ΦI )+ a.e. and when
ΦI has full column rank, it shares the same model subspace Mx̂ (y) = Im[IdI ] = Im[J]
as the hard- and soft-thresholding. Its invariant re-fitting is in this case Rinv
x̂ (y) =
IdI (ΦI )+ y. While the Lasso systematically underestimates the amplitude of the signal
by a shift λIdI ((ΦI )t ΦI )−1 sI , the re-fitting Rinv
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 sI are a.e. stable w.r.t. small perturbations [48]). It results that the
model subspace reads as Mx̂ (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):
(11)

+
Rinv
x̂ (y) = U (ΦU ) y .

8

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

(a) Noisy obs. y

(b) aniso-TV x̂(y)

(c) Ref. Rinv
x̂ (y) (= Rx̂ (y))

(d) iso-TV x̂(y)

(e) Inv. ref. Rinv
x̂ (y)

(f) Cov. ref. Rx̂ (y)

Fig. 3. (a) A noisy image y = x0 + w where x0 contains six overlapping squares. (b) The
solution x̂(y) of aniso-TV with λ = 10, and(c) its invariant re-fitting Rinv
x̂ (y) (which coincides in
this case with the covariant one Rx̂ (y)). (d) The solution x̂(y) of iso-TV with λ = 10, (e) its
invariant re-fitting Rinv
x̂ (y) and (f) its covariant one Rx̂ (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 Rinv
that performs a piece-wise average of
x̂ is the projector ΠIm[J] = U U
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 [49]
(12)

J = U (U > Φ> ΦU + λU > Γ> ΩΓU )−1 U > Φ> y

D
E

(
(Γx̂(y))i
(Γx̂(y))i
1
z
−
z
,
if i ∈ I,
i
i ||(Γx̂(y))i ||2
||(Γx̂(y))i ||2
where Ω : z ∈ Rm×b 7→ ||(Γx̂(y))i ||2
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 [49, 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., Mx̂ (y) = Ker[ΓI c ] = Im[J], and the
+
re-fitting is also given by Eq. (11), hence Rinv
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 [50, 18]. 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 [14] to compute the invariant re-fitting.
3.3. Results and limitations. Figure 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

9

CLEAR: Covariant LEAst-square Re-fitting

(a) Masked obs. y

(b) Tikhonov x̂(y)

(c) Inv. ref. Rinv
x̂ (y)

(d) Cov. ref. Rx̂ (y)

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

of all plateaus towards the data. Figure 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). Figure 3.(a)-(c) gives a similar illustration for a 2D aniso-TV denoising example.
While the invariant re-fitting acts properly for the `1 analysis estimator, it is
however less appealing in other scenarios. Figure 2.(c),(d) and Figure 3.(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 [6]), see Figure 2.(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 Figure 2.(d) and Figure 3.(e).
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 Rinv
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).
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.(e), Figure 3.(f) and Figure 4.(d).
4. 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

10

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

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 Φx0 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 Φx0 , 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 ∈ Rn of Φx0 , a new estimator Dx̂,z : Rn → Rp 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
Rx̂ . We are now equipped to introduce such a guess based re-fitting.

Definition 5. Let x̂ : Rn → Rp be differentiable at z ∈ Rn . An estimator Dx̂,z :
R → Rp is a guess based covariant least-square re-fitting of x̂ for z, if
n

(13)

Dx̂,z ∈ argmin ||Φh(z) − z||22 ,
h∈H

where H is the set of maps h : Rn
1. Affine map:
2. Covariant preserving:
3. Coherent map:

→ Rp satisfying, for all y ∈ Rn ,
h(y) = Ay + b
for some A ∈ Rp×n , b ∈ Rp ,
Jh (z) = ρJx̂ (z)
for some ρ ∈ R,
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 Dx̂,z (y).

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

Dx̂,z (y) = x̂(z) + ρJ(y − Φx̂(z))

where


 hΦJδ, δi
ρ=
||ΦJδ||22
 1

if

ΦJδ 6= 0 ,

otherwise ,

where J = Jx̂ (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 ∈ Rp . 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 ρ.

11

CLEAR: Covariant LEAst-square Re-fitting

(15)

Remark that Problem (13) can be recast as a one-dimensional problem
(
min
ρ∈R

||Φ(x̂(z) + ρJ(z − Φx̂(z))) −

z||22

= ||(Id − ρΦJ)(Φx̂(z) −

z)||22

)

whose unique solution, if ΦJ(Φx̂(z) − z) 6= 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δ||22 = hΦJδ, ΦJδi =
>
h(ΦJ) ΦJδ, δi = hΦJδ, δi. Injecting this in (14) gives ρ = 1.

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

Proposition 8. Let Y be a random vector in Rn such that E[Y ] = Φx0 , and
Cov[Y ] = Σ ∈ Rn×n . Then y 7→ Dx̂,z (y) satisfies
(16)

(17)
(18)

E[Dx̂,z (Y )] − x0 = (Id − ρJΦ)(x̂(z) − x0 ) ,

Cov[Dx̂,z (Y ), Y ] = ρJΣ ,

Cov[Dx̂,z (Y )] = ρ2 JΣJ > ,

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
(19)

E[(x̂(z) + ρJ(Y − Φx̂(z)))Y > ]−E[(x̂(z) + ρJ(Y − Φx̂(z)))]E[Y ]


>
= ρ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[Dx̂,z (Y )].
Proposition 8 provides a closed form expression for the bias, the cross-covariance
and the covariance of Dx̂,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̂.
(20)

Corollary 9. Let Tx̂,z (y) be the tangent estimator of x̂ at z ∈ Rn defined as
Tx̂,z (y) = x̂(z) + J(y − z) .

Let Y be a random vector in Rn such that E[Y ] = Φx0 and Cov[Y ] = Σ. Then
y 7→ Tx̂,z (y) and y 7→ Dx̂,z (y) satisfy
(21)

E[Tx̂,z (Y )] − x0 = (x̂(z) − x0 ) + J(Φx0 − z) ,

(22)

Cov[Dx̂,z (Y ), Y ] = ρCov[Tx̂,z (Y ), Y ] ,

(23)

Cov[Dx̂,z (Y )] = ρ2 Cov[Tx̂,z (Y )] ,

(24)

Corr[Dx̂,z (Y ), Y ] = Corr[Tx̂,z (Y ), Y ] ,

(25)

Corr[Dx̂,z (Y )] = Corr[Tx̂,z (Y )] ,

12

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

where the pcross correlation matrix is defined as Corr[X, Y ]i,j
=
Cov[X, Y ]i,j / Cov[X]i,i Cov[Y ]j,j , 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 Tx̂,z and that J does not
depend on y. It follows that Cov[Tx̂,z (Y ), Y ] = JΣ and Cov[Tx̂,z (Y )] = JΣJ > . These,
jointly with Proposition 8, conclude the proof.
Corollary 9 is essential in this work as it states that, by preserving the Jacobian
structure, Dx̂,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.
Proposition 10. Let Y be a random vector of Rn such that E[Y ] = Φx0 . Assume
ΦJ is an orthogonal projector, then y 7→ Dx̂,z (y) satisfies
(26)

||Φ(E[Dx̂,z (Y )] − x0 )||2 6 ||Φ(E[Tx̂,z (Y )] − x0 )||2 .

Proof. As ΦJ is an orthogonal projector, by virtue of Proposition 7, ρ = 1, then
(27) ||Φ(E[Tx̂,z (Y )] − x0 )||22 = ||Φ(x̂(z) + J(Φx0 − z) − x0 )||22

= ||Φ(x̂(z) + JΦ(x0 − x̂(z)) − x0 )||22 + ||ΦJ(Φx̂(z) − z)||22

which concludes the proof.

= ||Φ(E[Dx̂,z (Y )] − x0 )||22 + ||ΦJδ||22

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 Rn such that E[Y ] = Φx0 . Let
0 , ΦJδ0 i
ρ0 = hδ||ΦJδ
and δ0 = Φ(x0 − x̂(z)). Assume there exists α ∈ [0, 1] such that
2
0 ||
2

√
ρ − ρ0
6 1−α ,
ρ0

(28)
(29)

and

||ΦJ(δ − δ0 )||22 + 2hδ0 , ΦJ(δ − δ0 )i > −α

Then, y 7→ Dx̂,z (y) satisfies
(30)

hδ0 , ΦJδ0 i2
.
||ΦJδ0 ||22

||Φ(E[Dx̂,z (Y )] − x0 )||2 6 ||Φ(E[Tx̂,z (Y )] − x0 )||2 .

Proof. It follows from Proposition 8 and Corollary 9 that ||Φ(E[Dx̂,z (Y )]−x0 )||2 =
||(Id − ρΦJ)δ0 ||2 and ||Φ(E[Tx̂,z (Y )] − x0 )||2 = ||δ0 + ΦJ(δ − δ0 )||2 . Subsequently, we
get that Equation (30) holds true if
(31)
(32)

i.e.,

||(Id − ρΦJ)δ0 ||22 6 ||δ0 + ΦJ(δ − δ0 )||22 ,

ρ2 ||ΦJδ0 ||22 − 2ρhδ0 , ΦJδ0 i 6 ||ΦJ(δ − δ0 )||22 + 2hδ0 , ΦJ(δ − δ0 )i.

Using Assumption (29), a sufficient condition for Equation (30) to hold is
(33)

ρ2 ||ΦJδ0 ||22 − 2ρhδ0 , ΦJδ0 i2 + α

hδ0 , ΦJδ0 i2
60 .
||ΦJδ0 ||22

√
The roots of this polynomial are given by (1 ± 1 − α)ρ0 , which concludes the proof.

CLEAR: Covariant LEAst-square Re-fitting

13

Remark 12. Remark that requiring (28) is quite natural as it states that ρ should
be close enough to the optimal ρ0 minimizing the discrepancy with regards to Φx0
(i.e., minimizing ||Φh(z) − Φx0 ||22 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 = Φx0 , Assumption (29) holds true as it would read
0 > −αhδ0 , ΦJδ0 i2 /||ΦJδ0 ||22 (since δ = δ0 and ρ = ρ0 ). Another case of interest
is when ΦJ is an orthogonal projector for which (29) holds true as it would read
||ΦJδ||22 − ||ΦJδ0 ||22 > −||ΦJδ0 ||22 (using that ρ = ρ0 = 1, h·, ΦJ·i = ||ΦJ · ||22 , and
choosing α = 1). Hence, Proposition 11 recovers Proposition 10.
Remark 13. Using the same sketch of proof as Proposition 11, the condition
|(ρ − ρ0 )/ρ0 | 6 1 is sufficient to get ||Φ(E[Dx̂,z (Y )] − x0 )||2 6 ||Φ(x̂(z) − x0 )||2 . In other
words, even though ρ has a relative error of 100% w.r.t. ρ0 , the estimator y 7→ Dx̂,z (y)
still reduces the bias of the constant estimator y 7→ x̂(z). This result remains valid
when comparing y 7→ Dx̂,z (y) to the pseudo-oracle estimator y 7→ x̂(z) + J(y − Φx0 ),
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 Φx0 . In such cases, the estimator Dx̂,z can be
considered as a debiasing procedure of x̂, in the sense that it reduces the bias of Tx̂,z
while preserving its correlation structure (according to Corollary 9).
4.2. Definitions and properties. Using Dx̂,z defined in Theorem 6, we can
now give an explicit definition of CLEAR as Rx̂ (y) = Dx̂,y (y).

Definition 14 (CLEAR). The Covariant LEast-square Re-fitting associated to an
a.e. differentiable estimator y →
7 x̂(y) is, for almost all y ∈ Rn , given by

 hΦJδ, δi if ΦJδ 6= 0 ,
(34)
Rx̂ (y) = x̂(y) + ρJ(y − Φx̂(y)) with ρ =
||ΦJδ||22
 1
otherwise ,
where δ = y − Φx̂(y) and J = Jx̂ (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 Mx̂ (z) = Mx̂ (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 7→ Dx̂,z (y) for which ρ is constant. Note that the mapping
y 7→ Dx̂,z (y) is affine, but not the map y 7→ Rx̂ (y). Note that, as a consequence, the
statistical interpretations given in the previous section do not hold for Rx̂ (y) even
though they shed some light on its behavior.
Why not iterating the procedure? The re-fitting procedure computes x̃2 =
x̃1 + ρJ(y − Φx̃1 ), with x̃1 = x̂(y). One may wonder if it is beneficial to iterate the

14

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

Fig. 5. Geometrical illustration of the covariant re-fitting in a denoising problem of dimension
p = 3. We assume that Mx̂ (z) = Mx̂ (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.

process as x̃k+1 = x̃k + ρJ(y − Φx̃k ) (in the same vein as [47, 4, 33, 43, 37]). Consider
a denoising problem Φ = Id with Tikhonov or iso-TV, for which J is symmetrical
and x̂(y) ∈ Im[J] (see examples of Subsection 3.2). The sequence will converge if
and only if J(y − x̃k ) vanishes, i.e., x̃k must converge to J + Jy + ζ with ζ ∈ Ker[J].
By construction, x̃k ∈ 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 = Rinv
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 Rx̂ is the fact that it belongs to the model subspace
of x̂ as stated in the following proposition.
Proposition 17. Let y 7→ x̂(y) be an a.e. differentiable estimator. Then for
almost all y ∈ Rn , one has Rx̂ (y) ∈ Mx̂ (y).
Proof. As Mx̂ (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, Rx̂ (y) =
x̂(y) + J(y − Φx̂(y)), and, ΦRx̂ (y) = ΦRinv
x̂ (y).
Proof. By virtue of Proposition 7, ρ = 1 and then Rx̂ (y) = x̂(y) + J(y − Φx̂(y)).
+
The fact that ΦRx̂ (y) = ΦRinv
x̂ (y) comes from the fact that ΦJ(ΦJ) = ΦJ.

The next proposition provides, when JΦ satisfies a fixed point formulation, an
expression of Rx̂ (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 Rx̂ (y) = (1 − ρ)x̂(y) + ρJy.

15

CLEAR: Covariant LEAst-square Re-fitting

Proof. We have Rx̂ (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
(35)

x̂(y) = argmin F (y − Φx) + G(x) ,
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.
The affine constrained least-squares, the `1 synthesis, the `1 − `2 analysis, anisoTV 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 Rx̂ (y) = (1 − ρ)x̂(y) + ρJy.
4.3. 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 Rx̂ (y) = Rinv
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
Rx̂ (y) = (1 + ρ)x̂(y) − ρJΦx̂(y).

The soft- and hard-thresholding, used with Φ = Id, have the Jacobian matrix
J = IdI . As a consequence ΦJ = IdI 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 = IdI (ΦI )+ where ΦI has full column
rank. As for the thresholding, ΦJ = ΦU (ΦU )+ is an orthogonal projection and the
covariant re-fitting reads Rx̂ (y) = Rinv
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 Rx̂ (y) = Rinv
x̂ (y).
The `1 − `2 analysis has the Jacobian operator given in Eq. (12), which applied to
a vector d ∈ Rn is a solution of the following problem

(36)

Jd ∈

argmin

x ; supp(Γx) ⊆ I

where ω : z ∈ R

m×b

7→

X
i∈I

1
2 ||Φx

− d||22 + λ2 ω(Γx) ,

1
||(Γx̂(y))i ||2

||zi ||22



(Γx̂(y))i
− zi ,
||(Γ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 isoTV, 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 Rx̂ (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).

16

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

The non-local means has an intricate Jacobian matrix. Nevertheless, its directional
derivative has a simpler expression given, for any direction d ∈ Rn , by
P 0
P
P 0
X
w̄i,j dj − x̂(y)i j w̄i,j
j w̄i,j yj +
0
0
P
(37)
Jd =
with w̄i,j
=
wi−k,j−k
,
w̄
i,j
j
k

where x̂(y) is defined in Eq. (7) and
(38)

0
wi,j
= 2hPi y − Pj y, Pi d − Pj diϕ0 ||Pi y − Pj y||2



,

with ϕ0 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).
5. 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).
5.1. Algorithmic differentiation. Following [15], we consider restoration algorithms whose solutions x̂(y) = xk are obtained via an iterative scheme of the form
(
xk
= γ(ak ) ,
(39)
ak+1 = ψ(ak , y) .
Here, ak ∈ A is a sequence of auxiliary variables, ψ : A × Rn → A is a fixed point
operator in the sense that ak converges to a? , and γ : A → Rp is non-expansive
(i.e., ||γ(a1 ) − γ(a2 )|| 6 ||a1 − a2 ||, ∀a1 , a2 ∈ A) entailing xk converges to x? = γ(a? ).
As a result, for almost all y and for any direction d ∈ Rn , the directional derivatives Dxk = Jx̂k (y)d and Dak = Jak (y)d can be jointly obtained with xk and ak as
 k
x



ak+1

Dxk


 k+1
Da

(40)

where Γa =

∂γ(a)
∂a

ak

, Ψa =

∂ψ(a,y)
∂a

= γ(ak ) ,
= ψ(ak , y) ,
= Γa Dak ,
= Ψa Dak + Ψy d ,
ak

and Ψy =

∂ψ(ak ,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 (39). 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 [24, 31] for a comprehensive study.
5.2. Finite difference based differentiation. Another strategy is to approximate directional derivatives by finite differences, for any d ∈ Rn and ε > 0, as
(41)

Jx̂ (y)d ≈

x̂(y + εd) − x̂(y)
.
ε

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

17

CLEAR: Covariant LEAst-square Re-fitting

algorithm that provides x̂(y). For ε small enough, it performs as well as the approach
described in (40) (with x̂(y) = xk ) that requires the knowledge of the derivatives.
Indeed, if y 7→ x̂(y) is Lipschitz-continuous, then (41) converges to (40) when ε → 0
(by virtue of Rademacher’s theorem and [21, 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 Jx̂ (y)d. This finite difference approach has been used
in many fields, and notably for risk estimation, see e.g., [53, 40, 35]. We will apply
this black box strategy on state-of-the-art algorithms in Subsection 6.5.
5.3. Two-step computation for the general case. In the most general case,
the computation of the covariant re-fitting, given by
(42) Rx̂ (y) = x̂(y) + ρJ(y − Φx̂(y)) with

ρ=

hΦJδ, δi
||ΦJδ||22

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 (42). 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 (39), next by (40). 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.
5.4. One-step computation for specific cases. When x̂(y) fulfills the assumption JΦx̂(y) = x̂(y) of Proposition 19, the covariant re-fitting reads as
(43)

Rx̂ (y) = (1 − ρ)x̂(y) + ρJy

with ρ =

hΦ(Jy − x̂(y)), y − Φx̂(y)i
.
||Φ(Jy − x̂(y))||22

The computations of x̂(y) and Jy are then sufficient to compute the re-fitting Rx̂ (y).
As a result, in the case of algorithmic differentiation, (40) can be run once to get
Rx̂ (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 (39). 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.
5.5. Example on a primal-dual solver for `1 analysis. In this section we
instantiate Algorithm (40) to the case of the primal-dual sequence of [7]. By dualizing
the `1 analysis norm x 7→ λ||Γx||1 , the primal problem (4) can be reformulated, with
x? = x̂(y), as the following saddle-point problem
(44)

(z ? , x? ) ∈ arg max
minp
m
z∈R

x∈R

1
kΦx − yk22 + hΓx, zi − ιBλ (z) ,
2

where z ? ∈ Rm is the dual variable, and Bλ = {z ∈ Rm : ||z||∞ 6 λ} is the `∞ ball.
Problem (44) can be efficiently solved using the primal-dual algorithm of [7]. By
taking στ < 1/kΓk22 , θ ∈ [0, 1] and initializing (for instance,) x0 = v 0 = 0 ∈ Rp ,

18

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

z 0 = 0 ∈ Rm , the algorithm reads
 k+1
= ΠBλ (z k + σΓv k ) ,
 z

k+1
x
= (Id + τ Φ> Φ)−1 xk + τ (Φ> y − Γ> z k+1 ) ,
(45)
 k+1
v
= xk+1 + θ(xk+1 − xk ) ,

where the projection of z over Bλ is done component-wise as

zi
if |zi | 6 λ ,
(46)
ΠBλ (z)i =
λ sign(zi ) otherwise .

The sequence xk converges to a solution x? of the `1 analysis problem [7].
It is easy to check that the primal-dual sequence defined in (45) can be written
in the general form considered in (39), see for instance [15]. As a result, we can use
the algorithmic differentiation based strategy described by (40) as follows: for the
initialization x̃0 = ṽ 0 = 0 ∈ Rp , z̃ 0 = 0 ∈ Rm , and for β = 0, as
 k+1
z
= ΠBλ (z k + σΓv k ) ,




k+1

x
= (Id + τ Φ> Φ)−1 xk + τ (Φ> y − Γ> z k+1 ) ,


 k+1
v
= xk+1 + θ(xk+1 − xk ) ,
(47)
k+1
z̃
= Πzk +σΓvk (z̃ k + σΓṽ k ) ,




k+1

x̃
= (Id + τ Φ> Φ)−1 x̃k + τ (Φ> y − Γ> z̃ k+1 ) ,


 k+1
ṽ
= x̃k+1 + θ(x̃k+1 − x̃k ) ,

z̃i
if |zi | 6 λ + β ,
where Πz (z̃)i =
0
otherwise .
Recall that the re-fitting is Rxk (y) = x̃k , since JΦ is an orthogonal projector.
Remark that the algorithmic differentiation of (45) is exactly (47) for β = 0,
hence, x̃k = Rxk (y). However, if one wants to guarantee the convergence of the
sequence x̃k towards Rx̂ (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 rank2 . Let α > 0
be the minimum non zero value3 of |Γx? |i for all i ∈ [m]. Choose β such that ασ >
β > 0. Then, the sequence x̃k = Rxk (y) defined in (47) converges to the re-fitting
Rx̂ (y) of x̂(y) = x? .
The proof of this theorem is postponed to Appendix C.

A similar result was obtained in [13] when solving the `1 analysis problem (4)
with the Douglas-Rachford splitting algorithm described in [17, 8].
5.6. 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 ∈ Rm×b , as
(
zi
if ||zi ||2 6 λ ,
(48)
Πiso
Bλ (z)i =
zi
λ ||zi ||2 otherwise .
2 This
3 If

could be enforced as shown in [48].
|Γx? |i = 0 for all i ∈ [m], the result remains true for any α > 0.

19

CLEAR: Covariant LEAst-square Re-fitting

Algorithm Non-local means [3] and its directional derivative
Inputs:
Parameters:
Outputs:

noisy image y, direction d, noise standard-deviation σ
half search window width s, half patch width b, kernel function ϕ
x? = x̂(y) and x̃ = Jx̂ (y)d

Initialize W ← ϕ(2σ 2 (2b + 1)2 ) 1p1 ×p2
Initialize Wy ← ϕ(2σ 2 (2b + 1)2 ) y
Initialize W 0 ← 0p1 ×p2
Initialize Wy0 ← ϕ(2σ 2 (2b + 1)2 ) d
for k ∈ [−s, s]2 \ {0, 0} do
Compute e ← [ (y − Sk (y))2 ] ? κ
Compute w ← ϕ(e) ? κ
Update
W ←W +w
Update
Wy ← Wy + wSk (y)

(add weights for the central pixels [39])
(accumulators for the weighted sum)

(error between each k shifted patches [11, 12])
(contribution for each patch of its k shift)
(add weights at each position)
(add contribution of each k shifted patches)

Compute e0 ← [ 2(y − Sk (y))(d − Sk (d)) ] ? κ
Compute w0 ← [ e0 ϕ0 (e) ] ? κ
Update
W 0 ← W 0 + w0
Update
Wy0 ← Wy0 + w0 Sk (f ) + wSk (d)
end for
Compute x? ← Wy /W
Compute x̃ ← (Wy0 − W 0 x? )/W

(weighted mean)

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, Sk is the operator that shift all pixels in
the direction k, ? is the discrete convolution operator, and κ ∈ Rp1 ×p2 is such that κi,j = 1 if
(i, j) ∈ [−b, b]2 , 0 otherwise.

It follows that the algorithmic differentiation strategy reads as
 k+1
iso k
k
) ,


 z k+1 = ΠBλ (z +>σΓv

−1

x
=
(Id
+
τ
Φ
Φ)
xk + τ (Φ> y − Γ> (z k+1 )) ,

 k+1

v
= xk+1 + θ(xk+1 − xk ) ,
(49)
k+1
k
k
z̃
= Πiso

z k +σΓv k (z̃ + σΓṽ ) ,



k+1
>
−1
k


x̃
= (Id + τ Φ Φ)
x̃ + τ (Φ> y − Γ> z̃ k+1 ) ,

 k+1
ṽ
= x̃k+1 + θ(x̃k+1 − x̃k ) ,
(
z̃i 
iso
D
E
 if ||zi ||2 6 λ + β ,
where Πz (z̃)i =
zi
zi
λ
otherwise .
||zi ||2 z̃i − z̃i , ||zi ||2
||zi ||2

Unlike the `1 case, the re-fitted solution is not x̃k itself, but following Subsection 5.4,
it can be obtained at iteration k as
(50)

Rxk (y) = (1 − ρ)xk + ρx̃k

with ρ =

hΦ(x̃k − xk ), y − Φxk i
.
||Φ(x̃k − xk )||22

5.7. Example for the non-local means. In this section, we specify the update
rule (40) to an acceleration of the non-local means inspired from [11, 12]. We use the
procedure of [39] 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 (39),
where the fixed point solution is obtained directly at the first iteration.

20

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

PSNR: 23.49, SSIM: 0.218

(a)

(d)

PSNR: 29.98, SSIM: 0.900

PSNR: 37.64, SSIM: 0.985

(b)

(c)
Attenuated
structure

Re-tted
structure

Lost
structures

Left lost
structures

(e)

(f)

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

The pseudo-code obtained with the algorithmic differentiation scheme, as described in (40) is given in Figure 6. All variables with suffix 0 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(s2 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).
6. 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 PSNR4 and SSIM [51],
we numerically evaluate the re-fitting, discuss its benefit and limitations in several
scenarios and compare this method with popular approaches from the literature.

6.1. 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 [42], typically for thin detailed
structures, which are re-enhanced in our result.
4 PSNR

= 10 log10 2552 / p1 ||x̂(y) − x0 ||22 for an image ranging on [0, 255].

21

CLEAR: Covariant LEAst-square Re-fitting

PSNR: 9.51, SSIM: 0.261

(a)

PSNR: 25.21, SSIM: 0.734

(b)

PSNR: 27.28, SSIM: 0.831

(c)

Attenuated
structures
Residual
lost structures

(d)

(e)

(f)

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

The residuals x̂(y) − x0 and Rx̂ (y) − x0 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.
6.2. 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.
6.3. Denoising with non-local means. Figure 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

22

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

PSNR: 22.16, SSIM: 0.460

PSNR: 25.46, SSIM: 0.643

(a)

(b)

PSNR: 27.48, SSIM: 0.776

(c)

Attenuated
structures

Residual
lost structures

(d)

(e)

(f)

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

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.
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

23

CLEAR: Covariant LEAst-square Re-fitting

PSNR: 22.45, SSIM: 0.416

(a)

(c)

(b)

240

PSNR: 27.00, SSIM: 0.807

(d)

Original x̂(y)
Re-fitted Rx̂ (y)
Optimum original
Optimum re-fitted
Sub-optimum original
Sub-optimum re-fitted

220
200
Quadratic cost

PSNR: 24.80, SSIM: 0.545

180
160
140
120
100
80
0.2

0.4

0.6

0.8
1
1.2
1.4
Regularization parameter λ

1.6

1.8

2
PSNR: 28.89, SSIM: 0.809

(e)

(f)

PSNR: 22.84, SSIM: 0.312

(a)

PSNR: 35.99, SSIM: 0.938

(b)

(c)

PSNR: 28.28, SSIM: 0.823

(g)

PSNR: 23.96, SSIM: 0.694

(d)

Quadratic cost

10 2

10 1

10 0
1

2

3

4
5
6
7
Regularization parameter λ

(e)

8

9

10
PSNR: 38.22, SSIM: 0.935

(f)

PSNR: 46.42, SSIM: 0.986

(g)

Fig. 10. Experiment with aniso-TV: (top) poorly piece-wise constant case. (bottom) pure piecewise constant case. (a) Noise-free x0 . (b) Noisy y = x0 + w. (e) MSE of x̂(y) and its re-fitting
Rx̂ (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 λ.

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

24

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

PSNR: 22.18, SSIM: 0.397

(a)

(c)

(b)

PSNR: 26.62, SSIM: 0.724

(d)

Original x̂(y)
Re-fitted Rx̂ (y)
Optimum original
Optimum re-fitted
Sub-optimum original
Sub-optimum re-fitted

200

Quadratic cost

PSNR: 25.07, SSIM: 0.564

150

100

1

2

3

4

5

6
PSNR: 30.12, SSIM: 0.815

Filtering parameter h

(e)

(f)

PSNR: 22.18, SSIM: 0.847

(a)

PSNR: 24.68, SSIM: 0.914

(b)

(c)

PSNR: 29.20, SSIM: 0.823

(g)

PSNR: 24.64, SSIM: 0.910

(d)

240

Quadratic cost

220
200
180
160
140
120
1

2

3

Filtering parameter h

(e)

4

5

6
PSNR: 27.04, SSIM: 0.946

(f)

PSNR: 27.29, SSIM: 0.955

(g)

Fig. 11. Experiment with non-local means: (top) moderate patch redundancy case. (bottom)
high patch redundancy case. (a) Noise-free x0 . (b) Noisy y = x0 + w. (e) MSE of x̂(y) and its refitting Rx̂ (y) w.r.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 x0 is in agreement with the retrieved subspace model Mx̂ (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

25

CLEAR: Covariant LEAst-square Re-fitting

PSNR: 22.17, SSIM: 0.528

(a)

(b)

(c)

PSNR: 26.92, SSIM: 0.861

(d)

Original x̂(y)
Re-fitted R[x̂] (y)
Optimum original
Optimum re-fitted
Sub-optimum original
Sub-optimum re-fitted

200

Quadratic cost

PSNR: 25.00, SSIM: 0.7523

150

100

50
0

1

2

3

4
PSNR: 30.29, SSIM: 0.920

Filtering parameter log γ

(e)

(f)

PSNR: 22.16, SSIM: 0.452

(a)
200

(c)

(g)

PSNR: 26.60, SSIM: 0.721

(d)

Original x̂(y)
Re-fitted R[x̂] (y)
Optimum original
Optimum re-fitted
Sub-optimum original
Sub-optimum re-fitted

180

Quadratic cost

PSNR: 26.33, SSIM: 0.716

(b)

PSNR: 29.41, SSIM: 0.918

160
140
120
100
80
60
40
-2

-1

0

1

2

3

Filtering parameter log γ

(e)

4

5

6
PSNR: 31.02, SSIM: 0.858

(f)

PSNR: 29.91, SSIM: 0.845

(g)

Fig. 12. (top) Experiments with BM3D [10, 27] and with DDID [26] (bottom). (a) Noise-free
x0 . (b) Noisy y = x0 + w. (e) MSE of x̂(y) and its re-fitting Rx̂ (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 γ.

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 [10] and DDID [26] reaching state-of-the-art results in denoising. 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

26

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

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 [26]. For BM3D, we use the implementation of [27],
and choose to multiply by γ the inner parameters λhard
(in the hard thresholding
3D
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) 7→ xk 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” [2]. 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 Rx̂ (y).
Co-support identification based post re-fitting. Another solution referred to as
post re-fitting, and studied in, e.g., [19, 36, 1, 28, 1], consists in identifying the
(co-)support I = {i : (Γx̂(y))i 6= 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 xk , and unfortunately, supp(Γxk ) can
be far from supp(Γx̂(y)) even though xk can be made arbitrarily close to x̂(y). Such
erroneous support identifications lead to results that strongly deviates from Rx̂ (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

27

CLEAR: Covariant LEAst-square Re-fitting

PSNR: 22.20, SSIM: 0.428

(a)

PSNR: 26.04, SSIM: 0.747

(d)

(b)

PSNR: 26.68, SSIM: 0.712

(e)

PSNR: 28.15, SSIM: 0.830

(c)

PSNR: 28.87, SSIM: 0.819

(f)

Fig. 13. (a) Noise-free x0 . (b) Noisy y = x0 + 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.

PSNR: 20.11, SSIM: 0.266

(a)

PSNR: 23.31, SSIM: 0.728

(d)

(b)

PSNR: 22.12, SSIM: 0.698

(e)

PSNR: 23.38, SSIM: 0.730

(c)

PSNR: 23.77, SSIM: 0.742

(f)

Fig. 14. (a) Original x0 . (b) Blurred and noisy y = Φx0 + w. (c) Original aniso-TV x̂(y)
with λ = 0.7. Enhanced results by (d) iterative hard-thresholding, (e) post re-fitting with support
identification, and (f) our proposed covariant re-fitting.

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 Mx̂ (y). This is
appealing since the co-support of x̂(y) is difficult to identify, particularly in the analy-

28

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

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.
6.7. 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 [4] is a simple approach that consists in re-injecting to the current solution x̃k 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 x̃0 = 0, the iterations reads
(51)

x̃k+1 = x̃k + x̂(y − Φx̃k ) .

The first iterate is x̃1 = x̂(y), and x̃2 is known as the twicing estimate [47]. Such
approaches are popular in non-parametric statistics, e.g., for kernel smoothing [9].
When k increases, its bias tends to decreases while its variance increases, see
[43]. In denoising (i.e., when Φ = Id) with a linear estimator x̂(y) = W y (e.g., the
Tikhonov regularization), we get x̃k = (Id − (Id − W )k )y, for k > 0. In particular, the
twicing reads as x̃2 = (2W − W 2 )y. In this linear case, the covariant re-fitting reads
as Rx̂ (y) = (W + ρW − ρW 2 )y and coincides with the twicing when ρ = 1.

Iterative Bregman refinement. In [33], the authors proposed an iterative procedure, originally designed to improve iso-TV results, given by
!
k
X
k+1
k
(52)
x̃
= x̂ y +
(y − Φx̃ ) .
i=1

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 x̃k = (Id − (Id − W )k )y
and coincides with boosting (we refer to [5, 52, 23, 34] for related approaches).
SOS-boosting. In [37], 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
(53)

x̃k+1 = τ x̂(y + αΦx̃k ) − (τ α + τ − 1)x̃k .

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
x̃k = W y+α(W −Id)x̃k−1 , and in particular, for k = 2, we get x̃2 = (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 [37], 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.
SAIF-boosting. As described in [30], the diffusion of a filter
 consists in iteratively
re-applying the filter to the current estimate x̃k+1 = x̂ x̃k . The authors of [43]

29

CLEAR: Covariant LEAst-square Re-fitting

PSNR: 22.42, SSIM: 0.627

SOS-boosting Bregman

Boosting

(a)

(b)

PSNR: 24.89, SSIM: 0.773

(c)

PSNR: 26.81, SSIM: 0.829

(d)

PSNR: 25.20, SSIM: 0.785

PSNR: 25.44, SSIM: 0.794

PSNR: 25.99, SSIM: 0.809

PSNR: 26.71, SSIM: 0.793

PSNR: 27.01, SSIM: 0.828

PSNR: 24.68, SSIM: 0.718

PSNR: 22.70, SSIM: 0.638

PSNR: 22.42, SSIM: 0.627

PSNR: 24.85, SSIM: 0.768

PSNR: 24.86, SSIM: 0.770

PSNR: 24.85, SSIM: 0.769

PSNR: 24.86, SSIM: 0.769

(e) 1 iteration

(f) 2 iterations

(g) 5 iterations

(h) 20 iterations

Fig. 15. (a) Cartoon image x0 . (b) Noisy version y = x0 +w. (c) Iso-TV x̂(y) with λ = 1.5. (d)
Our covariant re-fitting Rx̂ (y). (e,f,g,h) From top to bottom, boosting [47, 4], Bregman iterations [5]
and SOS-boosting [37] at respectively 1, 2, 5 and 20 iterations.

noticed that, unlike the boosting method of [4], 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.
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

30

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

PSNR: 22.16, SSIM: 0.518

(a)

PSNR: 25.98, SSIM: 0.747

PSNR: 28.77, SSIM: 0.873

Boosting

(d)

PSNR: 26.69, SSIM: 0.779

PSNR: 27.05, SSIM: 0.794

PSNR: 27.35, SSIM: 0.798

PSNR: 25.15, SSIM: 0.669

Bregman

(c)

PSNR: 24.97, SSIM: 0.790

PSNR: 19.16, SSIM: 0.434

PSNR: 25.03, SSIM: 0.739

PSNR: 17.62, SSIM: 0.373

SOS-boosting

(b)

PSNR: 26.50, SSIM: 0.768

PSNR: 26.94, SSIM: 0.785

PSNR: 27.66, SSIM: 0.809

PSNR: 24.85, SSIM: 0.693

(e) 1 iteration

(f) 2 iterations

(g) 5 iterations

(h) 20 iterations

Fig. 16. (a) Image x0 with moderate patch redundancy. (b) Noisy version y = x0 +w. (c) Nonlocal means x̂(y) with h = 5.5. (d) Covariant re-fitting Rx̂ (y). (e,f,g,h) From top to bottom, boosting
[47, 4], Bregman iterations [5] and SOS-boosting [37] at respectively 1, 2, 5 and 20 iterations.

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 1n = 1n , (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

1 − |yλi | if |yi | > λ ,
Wii =
(54)
0
otherwise ,
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

1 if |yi | > λ ,
(55)
Jii =
0 otherwise ,

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
`1 analysis described in Equation (5), a candidate for W is
(56)

W = U (ΦU )+ − U (U > Φ> ΦU )−1 U > (Γ> )I R ,

CLEAR: Covariant LEAst-square Re-fitting

31

where R is a diagonal matrix with diagonal elements Rii = λ/|yi | for i ∈ I and 0 otherwise. One can check that if 1n ∈ 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 responsible for the systematic contraction of the `1 analysis regularization
(this simplifies to λ/|yi | 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 nonlocal means. Note that 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).
7. 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

32

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

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.
7.1. Acknowledgment. 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).
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||22 + ιC (x) where C = b + Im[A], b ∈ Rp and A ∈ Rp×n . The initial
problem can be recast as
(57)

argmin ||Φx − y||22 = b + A · argmin ||ΦAt − (y − Φb)||22 .
t∈Rn

x∈C

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 x = b + At, i.e., x = b + A(ΦA)+ (f − Φb) is a solution.

A.2. Retrieving the Tikhonov solution. We consider the minimization problem, defined for Γ ∈ Rm×p and λ > 0, of 21 ||Φx − y||22 + λ2 ||Γx||22 . The objective being
differentiable and convex, its first order optimality conditions give
(58)

Φ> (Φx − y) + λΓ> Γx = 0 ⇔ (Φ> Φ + λΓ> Γ)x = Φ> f .

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 miniλ2
1
2
mization problem, defined for λ > 0, of E(x,
 y) = 2 ||x − y||2 + 2 ||x||0 . The problem
is separable meaning that argmin E(x, y) i = argmin Ei (xi , yi ) with
x∈Rn

(59)

Ei (xi , yi ) =

1
2

(

xi ∈R

yi2
(xi − yi )2 + λ2

if xi = 0 ,
otherwise .

Since yi2 6 minxi [(xi − yi )2 + λ2 ] ⇔ |yi | 6 λ, we get
 2
1
yi
if |yi | 6 λ ,
(60)
min Ei (xi , yi ) =
λ2
otherwise ,
xi
2

which is reached by setting xi = 0 when |yi | 6 λ and xi = yi otherwise.

A.4. Retrieving the soft-thresholding solution. We consider the minimization problem, defined for λ > 0, of E(x, y) = 21 ||x − y||22 + λ||x||1 , which is, as for the
hard-thresholding, separable with Ei (xi , yi ) = 12 (xi − yi )2 + λ|xi |. By convexity, a
minimum is reached when zero belongs to its sub-differential. Hence, xi is solution if

sign(xi ) if |xi | > 0 ,
(61)
0 ∈ ∂Ei (xi , yi ) ⇔ xi ∈ yi − λ
[−1, 1]
otherwise ,
which holds true by setting xi = 0 when |yi | 6 λ and xi = yi − λ sign(yi ) otherwise.

CLEAR: Covariant LEAst-square Re-fitting

33

A.5. Retrieving the (block-wise) non-local means. We assume periodical
boundary conditions such that all quantities q indexed by i = (i1 , i2 ) ∈ Z2 satisfies
qi = q(i1 ,i2 ) = q(i1 +k1 p1 ,i2 +k2 p2 ) for all (k1 , k2 ) ∈ Z2 . This leads, for l ∈ [−b, b]2 , to the
following relationship
X
1X
1X
wi,j ||Pi x − Pj y||22 =
wi,j
F (x, y) =
(xi+l − yj+l )2
(62)
2 i,j
2 i,j
l
"
#
"
#
1X X
1X X
2
2
wi,j (xi+l − yj+l ) =
wi−l,j−l (xi − yj )
=
2 i,j
2 i,j
l
l
"
#
1X X
1X
=
wi−l,j−l (xi − yj )2 =
w̄i,j (xi − yj )2 .
2 i,j
2 i,j
l

For all i ∈ [p1 ] × [p2 ], studying the first optimality conditions gives
P
X
∂F (x, y)
j w̄i,j yj
w̄i,j (xi − yj ) = 0 ⇔ xi = P
=0⇔
.
(63)
∂xi
j w̄i,j
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

(64)

x̂(y) ∈ argmin F (y − Φx) + G(x) ,
x

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

(1 − ε)x̂(y) ∈ argmin F (y − εΦx̂(y) − Φx) + G(x) .
x

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
(66) F (y − εΦx̂(y) − Φv) + G(v) < F (y − εΦx̂(y) − (1−ε)Φx̂(y)) + G((1−ε)x̂(y)) ,
< F (y − Φx̂(y)) + (1−ε)G(x̂(y)) .

It follows that
(67)

F (y − εΦx̂(y) − Φv) + G(v) + εG(x̂(y)) < F (y − x̂(y)) + G(x̂(y)) .

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

F (y − Φw) + G(w) < F (y − Φx̂(y)) + G(x̂(y)) ,

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
x̂(y) − x̂(y − εΦx̂(y))
x̂(y) − (1 − ε)x̂(y)
= lim
= x̂(y) ,
ε→0
ε
ε
which concludes the proof.
(69)

JΦx̂(y) , lim

ε→0

34

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

Appendix C. Proof of Theorem 21.
Before turning to the proof of this theorem, let us introduce a first lemma.
Lemma 23. The re-fitting Rx̂ (y) of the `1 analysis regularization x? = x̂(y) is the
solution of the saddle-point problem
(70)

min max ||Φx̃ − y||22 + hΓx̃, z̃i − ιSI (z̃) ,

x̃∈Rp z̃∈Rm

where ιSI is the indicator function of the convex set SI = {p ∈ Rm : pI = 0} .

Proof. As ΦU has full column rank, the re-fitting of the solution (5) is the unique
solution of the constrained least-square problem (see Subsection 3.2)
(71)

Rx̂ (y) = U (ΦU )+ y = argmin ||Φx̃ − y||22 .
x̃∈Mx̂ (y)

Remark that x̃ ∈ Mx̂ (y) = Ker[IdtI cΓ] ⇔ (Γx̃)I c = 0 ⇔ ιS c (Γx̃) = 0, where SI c =
I
{p ∈ Rm : pI c = 0}.
Using Fenchel transform, ιS c (Γx̃) = maxz̃ hΓx̃, z̃i − ι∗S c (z̃), where ι∗S c is the
I
I
I
convex conjugate of ιS c . Observing that ιSI = ι∗S c concludes the proof.
I

I

Given Lemma 23, replacing Πzk +σΓvk in (47) by the projection onto SI , i.e.,

(72)

ΠSI (z̃)I c = z̃I c

and ΠSI (z̃)I = 0 ,

leads to the primal-dual algorithm of [7] applied to problem (70) which converges
to the re-fitted estimator Rx̂ (y). It remains to prove that the projection Πzk +σΓvk
defined in (47) converges to ΠSI in finite time.
Proof of Theorem 21. First consider i ∈ I, i.e., |Γx? |i > 0. By assumption on
α, |Γx? |i ≥ α > 0. Necessary zi? = λ sign(Γx? )i in order to maximize (44). Hence,
|z ? + σΓx? |i ≥ λ + σα. Using the triangle inequality shows that
(73)

λ + σα 6 |z ? + σΓx? |i 6 |z ? − z k |i + σ|Γx? − Γv k |i + |z k + σΓv k |i .

Choose ε > 0 sufficiently small such that σα − ε(1 + σ) > β. From the convergence
of the primal-dual algorithm of [7], the sequence (z k , xk , v k ) converges to (z ? , x? , x? ).
Therefore, for k large enough, |z ? − z k |i < ε, |Γx? − Γv k |i < ε, and
(74)

|z k + σΓv k |i > λ + σα − ε(1 + σ) > λ + β .

Next consider i ∈ I c , i.e., |Γx? |i = 0, where by definition |z ? |i 6 λ. Using again
the triangle inequality shows that
(75)

|z k + σΓv k |i 6 |z k − z ? |i + σ|Γv k − Γx? |i + |z ? |i .

Choose ε > 0 sufficiently small such that ε(1 + σ) < β. As (z k , xk , v k ) → (z ? , x? , x? ),
for k large enough, |z k − z ? |i < ε, |Γv k − Γx? |i < ε, and
(76)

|z k + σΓv k |i < λ + ε(1 + σ) 6 λ + β .

It follows that for k sufficiently large |z k + σΓv k |i 6 λ + β if and only if i ∈ I c , and
hence Πzk +σKvk (z̃) = ΠSI (z̃). As a result, all subsequent iterations of (47) will solve
(70), and hence from Lemma 23 this concludes the proof of the theorem.

CLEAR: Covariant LEAst-square Re-fitting

35

REFERENCES
[1] A. Belloni and V. Chernozhukov, Least squares after model selection in high-dimensional
sparse models, Bernoulli, 19 (2013), pp. 521–547.
[2] T. Blumensath and M. E. Davies, Iterative thresholding for sparse approximations, J. Fourier
Anal. Appl., 14 (2008), pp. 629–654.
[3] A. Buades, B. Coll, and J.-M. Morel, A review of image denoising algorithms, with a new
one, SIAM J. Multiscale Model. Simul., 4 (2005), pp. 490–530.
[4] P. Bühlmann and B. Yu, Boosting with the L2 loss: regression and classification, J. Am.
Statist. Assoc., 98 (2003), pp. 324–339.
[5] M. Burger, G. Gilboa, S. Osher, J. Xu, et al., Nonlinear inverse scale space methods,
Communications in Mathematical Sciences, 4 (2006), pp. 179–212.
[6] A. Chambolle, V. Duval, G. Peyré, and C. Poon, Geometric properties of solutions to the
total variation denoising problem, Preprint 1602.00087, Arxiv, 2016.
[7] A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with
applications to imaging, J. Math. Imaging Vis., 40 (2011), pp. 120–145.
[8] P. L. Combettes and J.-C. Pesquet, A Douglas-Rachford splitting approach to nonsmooth
convex variational signal recovery, IEEE J. Sel. Top. Sign, 1 (2007), pp. 564–574.
[9] P.-A. Cornillon, N. W. Hengartner, and E. Matzner-Løber, Recursive bias estimation for
multivariate regression smoothers, ESAIM: Probability and Statistics, 18 (2014), pp. 483–
502.
[10] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, Image denoising by sparse 3-D
transform-domain collaborative filtering, IEEE Trans. Image Process., 16 (2007), pp. 2080–
2095.
[11] J. Darbon, A. Cunha, T. F. Chan, S. Osher, and G. J. Jensen, Fast nonlocal filtering
applied to electron cryomicroscopy, in ISBI, 2008, pp. 1331–1334.
[12] C.-A. Deledalle, V. Duval, and J. Salmon, Non-local methods with shape-adaptive patches
(NLM-SAP), J. Math. Imaging Vis., 43 (2012), pp. 103–120.
[13] C.-A. Deledalle, N. Papadakis, and J. Salmon, Contrast re-enhancement of total-variation
regularization jointly with the Douglas-Rachford iterations, in SPARS, 2015.
[14] C.-A. Deledalle, N. Papadakis, and J. Salmon, On debiasing restoration algorithms: applications to total-variation and nonlocal-means, in SSVM, 2015, pp. 129–141.
[15] C.-A. Deledalle, S. Vaiter, G. Peyré, and J. M. Fadili, Stein unbiased gradient estimator
of the risk (SUGAR) for multiple parameter selection, SIAM J. Imaging Sci., 7 (2014),
pp. 2448–2487.
[16] D. L. Donoho and J. M. Johnstone, Ideal spatial adaptation by wavelet shrinkage,
Biometrika, 81 (1994), pp. 425–455.
[17] J. Douglas and H. Rachford, On the numerical solution of heat conduction problems in two
and three space variables, Trans. Amer. Math. Soc., (1956), pp. 421–439.
[18] V. Duval, J.-F. Aujol, and Y. Gousseau, A bias-variance approach for the nonlocal means,
SIAM J. Imaging Sci., 4 (2011), pp. 760–788.
[19] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, Least angle regression, Ann. Statist.,
32 (2004), pp. 407–499.
[20] M. Elad, P. Milanfar, and R. Rubinstein, Analysis versus synthesis in signal priors, Inverse
problems, 23 (2007), pp. 947–968.
[21] L. C. Evans and R. F. Gariepy, Measure theory and fine properties of functions, CRC Press,
1992.
[22] E. I. George, Minimax multiple shrinkage estimation, Ann. Statist., 14 (1986), pp. 188–205.
[23] G. Gilboa, A total variation spectral framework for scale and texture analysis, SIAM J. Imaging Sci., 7 (2014), pp. 1937–1961.
[24] A. Griewank and A. Walther, Evaluating derivatives: principles and techniques of algorithmic differentiation, Society for Industrial and Applied Mathematics (SIAM), 2008.
[25] A. E. Hoerl and R. W. Kennard, Ridge regression: Biased estimation for nonorthogonal
problems, Technometrics, 12 (1970), pp. 55–67.
[26] C. Knaus and M. Zwicker, Dual-domain image denoising, in ICIP, 2013, pp. 440–444.
[27] M. Lebrun, An analysis and implementation of the bm3d image denoising method, Image
Processing On Line, 2 (2012), pp. 175–213.
[28] J. Lederer, Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions,
arXiv preprint arXiv:1306.0113, (2013).
[29] Y. Lin and H. H. Zhang, Component selection and smoothing in multivariate nonparametric
regression, Ann. Statist., 34 (2006), pp. 2272–2297.
[30] P. Milanfar, A tour of modern image filtering: New insights and methods, both practical and

36

C-A. Deledalle, N. Papadakis, J. Salmon and S. Vaiter

theoretical, IEEE Signal Processing Magazine, 30 (2013), pp. 106–128.
[31] U. Naumann, Optimal Jacobian accumulation is NP-complete, Math. Program., 112 (2008),
pp. 427–441.
[32] G. R. Obozinski, M. J. Wainwright, and M. I. Jordan, High-dimensional support union
recovery in multivariate regression, in NIPS, 2008, pp. 1217–1224.
[33] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, An iterative regularization method
for total variation-based image restoration, SIAM J. Multiscale Model. Simul., 4 (2005),
pp. 460–489.
[34] S. Osher, F. Ruan, J. Xiong, Y. Yao, and W. Yin, Sparse recovery via differential inclusions,
Appl. Comput. Harmon. Anal., (2016).
[35] S. Ramani, T. Blu, and M. Unser, Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms, IEEE Trans. Image Process., 17
(2008), pp. 1540–1554.
[36] P. Rigollet and A. B. Tsybakov, Exponential screening and optimal rates of sparse estimation, Ann. Statist., 39 (2011), pp. 731–471.
[37] Y. Romano and M. Elad, Boosting of image denoising algorithms, SIAM J. Imaging Sci., 8
(2015), pp. 1187–1219.
[38] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259–268.
[39] J. Salmon, On two parameters for denoising with non-local means, Signal Processing Letters,
IEEE, 17 (2010), pp. 269–272.
[40] X. Shen and J. Ye, Adaptive model selection, J. Am. Statist. Assoc., 97 (2002), pp. 210–221.
[41] C. Stein, Inadmissibility of the usual estimator for the mean of a multivariate normal distribution, in Proc. 3th Berkeley Sympos. Math. Statist. and Prob., vol. 1, 1956, pp. 197–206.
[42] D. Strong and T. Chan, Edge-preserving and scale-dependent properties of total variation
regularization, Inverse problems, 19 (2003), pp. 165–187.
[43] H. Talebi, X. Zhu, and P. Milanfar, How to saif-ly boost denoising performance, IEEE
Trans. Image Process., 22 (2013), pp. 1470–1485.
[44] R. Tibshirani, Regression shrinkage and selection via the lasso, J. Roy. Statist. Soc. Ser. B,
(1996), pp. 267–288.
[45] R. J. Tibshirani, The lasso problem and uniqueness, Electron. J. Statist., 7 (2013), pp. 1456–
1490.
[46] A. N. Tikhonov, On the stability of inverse problems, Dokl. Akad. Nauk SSSR, 39 (1943),
pp. 176–179.
[47] J. W. Tukey, Exploratory data analysis, Reading, Mass., 1977.
[48] S. Vaiter, C.-A. Deledalle, G. Peyré, C. Dossal, and J. Fadili, Local behavior of sparse
analysis regularization: Applications to risk estimation, Appl. Comput. Harmon. Anal., 35
(2013), pp. 433–451.
[49] S. Vaiter, C.-A. Deledalle, G. Peyré, J. M. Fadili, and C. Dossal, The degrees of freedom
of partly smooth regularizers, to appear in Ann. Inst. Stat. Math., (2016).
[50] D. Van De Ville and M. Kocher, Nonlocal means with dimensionality reduction and surebased parameter selection, IEEE Trans. Image Process., 20 (2011), pp. 2683–2690.
[51] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image quality assessment: from
error visibility to structural similarity, IEEE Trans. Image Process., 13 (2004), pp. 600–612.
[52] J. Xu and S. Osher, Iterative regularization and nonlinear inverse scale space applied to
wavelet-based denoising, IEEE Trans. Image Process., 16 (2007), pp. 534–544.
[53] J. Ye, On measuring and correcting the effects of data mining and model selection, J. Am.
Statist. Assoc., (1998), pp. 120–131.
[54] M. Yuan and Y. Lin, Model selection and estimation in regression with grouped variables, J.
Roy. Statist. Soc. Ser. B, 68 (2006), pp. 49–67.

