Jean-François Aujol
email: jean-francois.aujol@cmla.ens-cachan.fr

Some first-order algorithms for total variation based image restoration

Keywords: Algorithms, duality, total variation regularization, image restoration AMS subject classifications: 68U10, 49M29, 65K10

Some first-order algorithms for total variation based image restoration

Introduction

During the last 15 years, total variation regularization has known a great success in image processing [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing[END_REF][START_REF] Chan | Image processing and analysis -Variational, PDE, wavelet, and stochastic methods[END_REF][START_REF] Andreu-Vaillo | Parabolic quasilinear equations minimizing linear growth functionals[END_REF]. It has been used in many applications such as image restoration, image deblurring, image zooming, image inpainting, . . . (see [START_REF] Chan | Image processing and analysis -Variational, PDE, wavelet, and stochastic methods[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing[END_REF] and references therein). In all these approaches, a total variation term |Du| is to be minimized in some way. The typical problem is the case of image restoration [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] with the minimization of a functional of the type:

Ω |Du| + 1 2µ f -u 2 (1)
|Du| stands for the total variation of u [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF], and if u is regular it is simply Ω |∇u| dx. Ω is the image domain, a convex Lipschitz open set in R 2 . f is the degraded image to restore. The minimizer u of functional [START_REF] Acar | Analysis of total variation penalty methods for ill-posed problems[END_REF] is the restored image we want to compute (see for instance [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF] for a thorough mathematical analysis of this problem). µ is a weighting parameter which controls the amount of denoising. In the case of zero mean Gaussian noise, µ can can be related to the standard deviation of the noise.

From a numerical point of view, total variation is not straightforward to minimize, since it is not differentiable in zero. A first approach is to regularize it, and instead to consider a term as β 2 + |∇u| 2 dx. We will refer to this choice as smoothed total variation regularization:

Ω β 2 + |∇u| 2 dx + 1 2µ f -u 2 (2)
The classical approach is then to use the associated Euler-Lagrange equation to compute the solution. Fixed step gradient descent [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], or later quasi-Newton methods [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Dobson | Convergence of an iterative method for total variation denoising[END_REF][START_REF] Acar | Analysis of total variation penalty methods for ill-posed problems[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Ng | On semismooth Newton methods for total variation minimization[END_REF][START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF] have been proposed for instance (see [START_REF] Chan | Image processing and analysis -Variational, PDE, wavelet, and stochastic methods[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing[END_REF] and references therein). Iterative methods have proved successful [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Bect | A l1-unified variational framework for image restoration[END_REF]. A projected-subgradient method can be found in [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF].

Ideas from duality have also been proposed: first by Chan and Golub [START_REF] Chan | A nonlinear primal-dual method for total variation-based image restoration[END_REF], later by A. Chambolle in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF], and then generalized in [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. Chambolle's projection algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] has grown very popular, since it is the first algorithm solving exactly problem [START_REF] Acar | Analysis of total variation penalty methods for ill-posed problems[END_REF] and not an approximation like [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF], with a complete proof of convergence. Moreover, it is straightforward to implement it. In [START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF], a very interesting combination of the primal and dual problems has been introduced. Second order cone programming ideas and interior point methods have proved interesting approaches [START_REF] Goldfarb | Second-order cone programming methods for total variation based image restoration[END_REF][START_REF] Fu | Efficient minimization methods of mixed l1-l1 and l2-l1 norms for image restoration[END_REF]. Recently, it has been shown that graph cuts based algorithms could also be used [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF][START_REF] Darbon | Image restoration with discrete constrained total variation part I: Fast and exact optimization[END_REF]. Finally, let us notice that it is shown in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] that Nesterov's schemes [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] provides fast algorithms both for minimizing functional (1) and [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF].

In this paper, we revisit Chambolle's projection algorithm. We show that a modification of Chambolle projection algorithm, recently suggested in [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF], can be seen as a particular instance of a more general algorithm proposed almost 30 years ago by Bermùdez and Moreno [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF]. It is in fact an adaptation of Uzawa algorithm [START_REF] Ciarlet | Introduction à l'analyse numérique matricielle et à l'optimisation[END_REF] to problem [START_REF] Acar | Analysis of total variation penalty methods for ill-posed problems[END_REF]. This is the first main contribution of the paper: shedding some new light on these projection based type algorithms. We then apply the approach of Bermùdez and Moreno to smoothed total variation regularization: this gives a new fast algorithm to minimize functionals such as [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF]. This is the second main contribution of the paper. We also prove the convergence of this new scheme. Notice that Bermùdez and Moreno algorithm has already been used for smoothed total variation based restoration in [START_REF] Almansa | A TV based restoration model with local constraints[END_REF], but with a different numerical scheme. To test the efficiency of these algorithms, inspired by [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF], we have decided to make some comparisons with a general class of efficient minimization algorithms introduced by Y. Nesterov in [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. It has been proved in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] that they are indeed very efficient for image restoration. We chose to use these type of algorithms, because as in the case of Bermùdez and Moreno approach, they consist in first order schemes and it is proved in [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] that they are optimal (in the sense that no algorithms, using only the values and gradients of the functional to minimize, has a better rate of convergence [START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF]). We also explain how a recent improvement of these algorithms in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] can be applied for image restoration. We give some numerical examples of all the schemes introduced in this paper: this is the third main contribution of the paper. Our experiments are in favor of Bermùdez-Moreno approach to get a fast approximation for smoothed total variation regularization, whereas Nesterov schemes seem to perform better for total variation regularization. Notice that to get a highly accurate solution for smoothed total variation regularization, Nesterov's schemes seem also to be the best choice. However, such an accuracy is not necessary for image restoration.

Before presenting the plan of the paper, let us emphasize once more the main contributions of the paper:

• Shedding some new light on the Chambolle projection algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], by seeing how it can be related to a particular instance of Bermùdez and Moreno algorithm [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF].

• Introducing a new and efficient scheme for smoothed total variation based image restoration.

• Presenting numerous numerical comparisons with a general class of algorithms recently introduced by Nesterov [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF].

The organization of the paper is the following. In Section 2, we recall Bermùdez-Moreno algorithm [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF]. We show how it can be applied to total variation regularization in Section 3. We also explain the relations between this scheme and Chambolle's projection algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], and we give some numerical examples. In Section 4, we detail how Bermùdez-Moreno algorithm can be applied to smoothed total variation based image restoration, providing a new algorithm to solve this type of problem. We then explain in Section 5 how these schemes can be used for image deblurring. Whereas Sections 2 to 5 are related to applications of Bermùdez and Moreno framework, Section 6 concerns a different type of approach (and thus this section can be read independently). In Section 6, we recall a general class of minimization algorithms introduced by Y. Nesterov in [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. These algorithms have proved very efficient in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] for solving image processing problems. We then explain how a recent improvement of these algorithms in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] can be applied for image restoration. In Section 7 we make some comparisons between the different schemes presented in this paper. Appendix A details the proof of convergence of Bermùdez-Moreno algorithm.

Bermùdez-Moreno algorithm

In this section, we present the algorithm proposed by Bermùdez and Moreno in [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF]. This is a general minimization algorithm. Surprisingly, this approach seems to have been ignored by the image processing community, although it provides efficient algorithms for solving classical image processing problems as we will see in the next sections of the paper. In particular, it gives an algorithm to solve problem (1) without resorting to some smooth approximation like problem [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF]. Notice that A. Chambolle's paper [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] with its projection algorithm to solve problem (1) (the first algorithm to solve exactly (1)) was published 12 years after the seminal work of Rudin at al [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. And yet, [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF], which was published 11 years before [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], already provided a similar algorithm with a proof of convergence. We follow here the presentation of [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF] and [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] (chapter II.3). The general minimization problem considered is the following (and we will see that this framework can be used in many image processing problems):

inf z∈V 1 2 Az, z -g, z + ψ(z) (3)
with V Hilbert space, ψ a proper convex lower semi continuous (l.s.c.) function defined on V :

ψ = φ o B * ( 4
)
where E is a Hilbert space, B a bounded linear operator, B :

E → V , B * : V → E, φ : E → R.
We recall that if H is a convex function, we say that it is proper if H(x) > -∞ for all x, and if there exists x 0 such that H(x 0 ) < +∞. We denote by dom H the set on which H(x) < +∞ [START_REF] Brezis | Analyse fonctionnelle. Théorie et applications[END_REF][START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF].

Assumptions on A: In all the paper, we will make the following assumptions on A: A is assumed to be a linear symmetric coercive operator, i.e. there exists α > 0 such that for all z in V :

Az, z V ≥ α z 2 V ( 5
)
Notice that it implies in particular that A is a monotone operator, i.e. Ay -Az, y -z ≥ 0 for all y, z ∈ V . We will also make the following assumptions on A:

A is continuous on the finite dimensional subspaces of V . There exists z 0 in dom ψ such that:

Az,z-z 0 +ψ(z) z → +∞ if z → +∞ (6)
Notice that in the next sections, all these assumptions will indeed be satisfied. In particular, since we will only consider operators A of the type A = γI for some γ > 0, the technical assumption (6) will be trivially verified.

Notations:

We use the following notations [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF]. If H is a maximal monotone operator, we denote by H λ its Yosida approximation (L λ is the resolvent of λH):

H λ = I -L λ λ where L λ = (I + λH) -1 (7)
Bermùdez and Moreno derive their results for H = ∂φ -ωI. Here we choose ω = 0, and we take the operator H as (φ being defined in (4)):

H = ∂φ (8)
Notice that since φ is assumed to be a convex proper lower semi continuous function, its subdifferential ∂φ is a maximal monotone operator [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF][START_REF] Brezis | Analyse fonctionnelle. Théorie et applications[END_REF][START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF][START_REF] Andreu-Vaillo | Parabolic quasilinear equations minimizing linear growth functionals[END_REF]. We had to recall the notions of subdifferential of a convex function, maximal monotone operator, and Yosida approximations, because Bermùdez and Moreno approach is based on convex analysis. The algorithm they propose to compute the solution of (3) relies on the associated Euler-Lagrange equation (which in this case happens to be a subdifferential inclusion). See Appendix A for further details.

Algorithm: In [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF], Bermùdez and Moreno propose to use the following algorithm to minimize (3). y 0 being arbitrary, consider the iterative scheme:

u m = A -1 (g -By m ) y m+1 = H λ (B * u m + λy m ) (9)
They prove the following convergence result (proposition 3.1 in [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF]):

Theorem 1 Assume that A is a linear symmetric coercive operator satisfying ( 5) and ( 6), and that φ is a convex proper lower semi continuous function. Assume furthermore that:

0 < 1 λ < 2α B * 2 (10)
Then the sequence (u m ) defined by ( 9) is such that: lim m→+∞ u m = u (for the strong topology of V ) with u solution of: g -Au ∈ B∂φ(B * u), i.e. u unique solution of (3). Moreover, y m ⇀ y in E weak, with: y ∈ ∂φ(B * u).

The proof of Theorem 1 is detailed in Appendix A.

Relation with Forward-Backward Splitting : It was pointed out to the author by one of the anonymous reviewer that Bermùdez-Moreno algorithm can be seen as a particular instance of Forward-Backward Splitting applied to the dual problem of (3); see also Remark 3.2 in [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF]. Theorem 1 can then be deduced as a consequence of results from e.g. [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

One of the main interest of Theorem 1 is that it is not restricted to the case when φ is a support function [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF]. However, due to the importance of total variation regularization in image processing, we first consider the case of problem (1) in Sections 3 and 3.4. We will consider the case of problem (2) in Section 4, where φ is no longer a support function.

Application to total variation regularization

In this section, we show how Bermùdez-Moreno algorithm (9) can be used for total variation regularization. In Section 3.1, we first consider the continuous setting to derive the link with Bermùdez-Moreno's work. In Section 3.2, we then consider the discrete case and we show that Bermùdez and Moreno algorithm consists in solving the dual problem of (1) with a projected gradient algorithm, whose convergence is guaranteed thanks to Theorem 1. We give some numerical examples in Section 3.3. In Section 3.4, we explain the connection between Bermùdez-Moreno framework and other existing approaches.

Continuous setting

Let us consider the celebrated ROF model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]:

inf u∈L 2 (Ω) J(u) + 1 2µ f -u 2 L 2 (Ω) (11)
Here J(u) is the total variation of u extended to L 2 (Ω) (since in dimension 2, we have BV (Ω) ⊂ L 2 (Ω) [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF]):

J(u) = Ω |Du| if u ∈ BV (Ω) +∞ otherwise (12)
In fact, (11) is a particular case of (3). Indeed, take

V = L 2 (Ω), E = (L 2 (Ω)) 2 , A = 1 µ I, g = 1 µ f . A is of course coercive with coercivity constant α = 1 µ . J(u) = ψ(u) = φ(B * (u)), and
J(u) = sup v∈K u, div v (13)
Hence φ is the support function of K (closed convex set in (L 2 (Ω)) 2 ):

K = v ∈ (L 2 (Ω)) 2 / div v ∈ L 2 (Ω) , v ∞ ≤ 1 with |v| = v 2 1 + v 2 2 ( 14
)
We have:

φ(w) = sup v∈K w, v (L 2 (Ω)) 2 , B = -div = ∇ * , and B * = ∇ (15)
where we have used the fact that K is symmetric to 0. We recall that w, v (L 2 (Ω)

) 2 = w 1 , v 1 L 2 (Ω) + w 2 , v 2 L 2 (Ω) .
Moreover, since φ is the support function of K, then H λ (v) is the orthogonal projection of v λ onto K [44, 10] , i.e.: H λ (v) = P K v λ , where if x = (x 1 , x 2 ),

P K (x) = x 1 max{1, |x|} , x 2 max{1, |x|} (16)
Bermùdez-Moreno algorithm [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF] in this case is: u 0 arbitrary, and:

u m = f + µdiv y m y m+1 = P K y m + 1 λ ∇u m (17)
Applying Theorem 1, we get the following result:

Proposition 1 If λ > µ 2 B * 2
, then the sequence (u m , y m ) defined by scheme [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] is such that u m → u (in L 2 (Ω) strong) and y m ⇀ y (in L 2 (Ω) × L 2 (Ω) weak), with u solution of (11).

Discrete setting

From now on, and until the end of the paper, we will restrict our attention to the discrete setting. We take here the same notations as in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. The image is a two dimension vector of size N × N . We denote by X the Euclidean space R N ×N , and Y = X × X. The space X will be endowed with the inner product (u, v) = 1≤i,j≤N u i,j v i,j and the norm u = (u, u). To define a discrete total variation, we introduce a discrete version of the gradient operator. If u ∈ X, the gradient ∇u is a vector in Y given by: (∇u) i,j = ((∇u) 1 i,j , (∇u) 2 i,j ). with (∇u)

1 i,j = u i+1,j -u i,j if i < N 0 if i = N and (∇u) 2 i,j = u i,j+1 -u i,j if j < N 0 if j = N
The discrete total variation of u is then defined by:

J(u) = 1≤i,j≤N |(∇u) i,j | (18)
We also introduce a discrete version of the divergence operator. We define it by analogy with the continuous setting by div = -∇ * where ∇ * is the adjoint of ∇: that is, for every p ∈ Y and u ∈ X, (-div p, u) X = (p, ∇u) Y . It is easy to check that:

(div (p)) i,j =    p 1 i,j -p 1 i-1,j if 1 < i < N p 1 i,j if i = 1 -p 1 i-1,j if i = N +    p 2 i,j -p 2 i,j-1 if 1 < j < N p 2 i,j if j = 1 -p 2 i,j-1 if j = N (19)
From now on, we will use these discrete operators. Notice that in all the rest of the paper (except in the appendix), we place ourself in the discrete setting. We will sometimes use continuous notations; however, the reader has to keep in mind that only the discrete case is considered.

We will use Meyer G space for oscillating patterns [START_REF] Meyer | Oscillating patterns in image processing and in some nonlinear evolution equations[END_REF][START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF]:

G = {v ∈ X / ∃g ∈ Y such that v = div (g)} ( 20
)
and if v ∈ G: v G = inf g ∞ / v = div (g), g = (g 1 , g 2 ) ∈ Y, |g i,j | = (g 1 i,j ) 2 + (g 2 i,j ) 2 (21)
where g ∞ = max i,j |g i,j |. Moreover, we will use the notation:

G µ = {v ∈ G / v G ≤ µ} (22)
With these classical finite differences, we have: ∇u 2 ≤ 8 u 2 . Hence ∇ 2 = ∇ * 2 ≤ 8. And in fact it is possible to show [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] that

∇ 2 = ∇ * 2 = 8.
Let us consider new variables: Then we can rewrite [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] into: p 0 arbitrary, and

v m = u m µ , p m = y m , τ = µ λ (23)
v m = f µ + div p m p m+1 = P K (p m + τ ∇v m ) (24)
Applying Theorem 1, we get the following result:

Proposition 2 Let X the Euclidean space R N ×N , and Y = X × X. If τ < 1 4 , then the sequence (v m , p m ) defined by scheme [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] is such that v m → v in X and p m → p in Y with µv solution of [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF].

Notice that (24) can be written in a more compact way:

p m+1 = P K p m + τ ∇ f µ + div p m (25)

Numerical examples

We show here some numerical experiments with scheme [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. We will make some comparisons with other existing algorithms in Section 7 and study their numerical accuracy. On Figure 1, we display the classical images Lenna and cameraman that we use in this paper to illustrate our study. We also show their noisy versions (degraded by additive zero mean Gaussian noise with standard deviation σ = 20, the dynamic range of the gray values of the image being [0,255]). On Figure 2, we show the restoration we get with [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. These results have the classical behavior of total variation based image restoration.

Figure 2: Total variation restoration of the noisy images presented on the bottom row of Figure 1 with scheme [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. In both cases, the Lagrange multiplier is µ = 30.

Relation with Chambolle projection algorithm

In [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], A. Chambolle proposes a nonlinear projection algorithm to minimize the ROF model ( 11). This algorithm is based on the remark that the solution of ( 11) is given by u = f -P Gµ (f ), where P Gµ is the orthogonal projector onto G µ (defined by ( 22)). [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] gives an algorithm to compute P Gµ (f ). It indeed amounts to finding:

min µdiv (p) -f 2 X : p / |p i,j | ≤ 1 ∀i, j = 1, . . . , N (26)
This problem can be solved by a fixed point method: v 0 = 0, p 0 = 0, and

v m = f µ + div p m p m+1 i,j = p m i,j +τ (∇v m )i,j 1+τ |(∇v m )i,j | (27)
It is shown in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] that if τ < 1/8 in ( 27), then µv m converges to the solution of [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]. In practice, convergence of ( 27) is generally observed as long as τ < 1/4. An extension of this algorithm to color images has been proposed in [START_REF] Bresson | Fast minimization of the vectorial total variation norm and applications to color image processing[END_REF]. The case of more general Hilbert space has been considered in [START_REF] Aujol | Constrained and SNR-based solutions for TV-Hilbert space image denoising[END_REF].

In [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF], A. Chambolle has proposed a modification of his projection algorithm. Instead of using [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], he suggests in [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] to use a simple projected gradient method to compute the projection P Gµ :

   v m = f µ + div p m p m+1 i,j = p m i,j +τ (∇v m )i,j max{1,|p n i,j +τ (∇v m )i,j |} (28)
And this last equation is exactly scheme [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. In [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF], A. Chambolle has proved the stability of (28). However, since the functional is not elliptic [START_REF] Polyak | Introduction to optimization[END_REF], the convergence of the projection algorithm is not straightforward. In this paper, the convergence of (28) provided τ < 1/4 is a consequence of Proposition 2. See also [START_REF] Duval | Projected gradient based color image decomposition[END_REF] where a direct proof of convergence (inspired by this work)) of a similar projection algorithm is proposed. Notice that a partial proof of convergence of the projection algorithm has independently been proposed in [START_REF] Zhu | Duality-based algorithms for total variation image restoration[END_REF]: the authors get the same type of result as the one of Proposition 2 here, but with only the convergence of v m in (24) (in their result, the sequence p m is not guaranteed to converge). Moreover, as we will see in the next section, the general algorithm [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF] proposed by Bermùdez and Moreno [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF] can be of interest to other image restoration problems, such as smoothed total variation regularization based ones [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF]. Numerical comparisons of all these schemes (( 24), ( 27)) will be discussed in Section 7.2. The fact that in the case of the ROF problem (11), Bermùdez-Moreno algorithm (scheme (( 24)) is just a projected gradient algorithm on the dual problem has many implications:

1. Let us notice that problem [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF] is of the type:

inf u∈Q E(u) ( 29
)
where E is a convex Lipschitz non differentiable function, and Q a convex closed set. For this type of problem, it can be shown [START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF] (Theorem 3.2.1 page 138) that no algorithm (only using the values and gradients of

E) has a better rate of convergence than O 1 √ k
(in term of objective function) uniformly on all problems of the form [START_REF] Durand | Multiplicative noise removal using L 1 fidelity on frame coefficients[END_REF], with k the number of iterations of the algorithm. Nevertheless, it is also proved in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] (see also [START_REF] Weiss | Algorithmes rapides d'optimisation convexe. Applications à la restauration d'images et à la détection de changement[END_REF] Theorem 3.12 page 36) that the projected gradient method for minimizing a convex Lipschitz differentiable functional on a closed convex set is of order O 1 k . Scheme ( 24) is therefore an algorithm of order O 1 k for solving [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]. 2. It is well-known that the projected gradient algorithm is a particular instance of the proximal forwardbackward algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. This provides a general framework for minimizing the sum of two convex functions.

The convergence result of Proposition 2 could also be derived from [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. This approach has been used for instance in [START_REF] Yuan | Convex Hodge decomposition and regularization of image flows[END_REF] to prove the convergence of a similar algorithm to [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. The connection between the projected gradient algorithm and the proximal forward-backward algorithm is emphasized in [START_REF] Durand | Multiplicative noise removal using L 1 fidelity on frame coefficients[END_REF][START_REF] Fadili | Monotone operator splitting for fast sparse solutions of inverse problems[END_REF]. Of course, as explained at the end of Section 2, the relation (in the particular case of total variation regularization) between Bermùdez-Moreno algorithm and the proximal forward-backward algorithm comes from the fact that in general Bermùdez-Moreno framework is a particular instance of Forward Backward Splitting. See also [START_REF] Setzer | Split Bregman Algorithm, Douglas-Rachford Splitting and frame shrinkage[END_REF] where the connection is made between Forward-Backward Splitting and the Split Bregman algorithm recently proposed in [START_REF] Goldstein | The split Bregman algorithm for L1 regularized problems[END_REF] 4 Smoothed total variation regularization

In this section, we consider the following problem:

inf u Ω β 2 + |∇u| 2 dx + 1 2µ f -u 2 2 (30)
We refer to this problem as the smoothed total variation based regularization problem. For small values of β it can be seen as an approximation of [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]. This type of regularization is very common in image processing (see [START_REF] Aubert | Mathematical Problems in Image Processing[END_REF][START_REF] Chan | Image processing and analysis -Variational, PDE, wavelet, and stochastic methods[END_REF] and references therein). Compared to total variation regularization, it has the advantage of being a smooth regularization. And compared to stronger regularization such as ∇u 2 , it has the advantage of not eroding too much the edges of the image.

In Section 4.1, we explain how Bermùdez-Moreno algorithm (9) can be used to solve this problem. The new algorithm we propose has a fixed point iteration step. We show the convergence of this fixed point iteration in Section 4.3. We will show some numerical examples with this new scheme in Section 7.1.

Presentation of the scheme

Let us denote by

φ β (ξ) = Ω β 2 + |ξ| 2 dx (31)
We have

∂φ β (ξ) = ξ β 2 + |ξ| 2 (32)
Let us consider the following scheme:

u m = f + µdiv y m y m+1 = I-(I+λ∂φ β ) -1 λ (∇u m + λy m ) ( 33
)
Applying Theorem 1, we get:

Proposition 3 Let X the Euclidean space R N ×N , and Y = X × X. If λ > 4µ, then the sequence (u m , y m ) defined by scheme [START_REF] Fadili | Monotone operator splitting for fast sparse solutions of inverse problems[END_REF] is such that u m → u in X and y m → y in Y with u solution of [START_REF] Duval | Projected gradient based color image decomposition[END_REF].

The second equation of (33) implies:

λy m+1 = ∇u m + λy m -(I + λ∂φ β ) -1 (∇u m + λy m ) (34)
As in the total variation case, let us set:

v m = u m µ and τ = µ λ and y m = p m (35)
(33) becomes:

v m = f µ + div p m (I + λ∂φ β ) λ(τ ∇v m + p m -p m+1 ) = λ(τ ∇v m + p m ) ( 36
)
Let us set:

w m+1 = τ ∇v m + p m -p m+1 (37)
From the second line of (36), we get:

w m+1 + ∂φ β λ(w m+1 ) = τ ∇v m + p m (38)
But

∂φ β λw m+1 = λw m+1 β 2 + |λw m+1 | 2 = w m+1 β 2 λ 2 + |w m+1 | 2 (39)
We thus get from ( 38)

w m+1 + w m+1 β 2 τ 2 µ 2 + |w m+1 | 2 = τ ∇v m + p m (40)
Using the notations γ = βτ µ , and C m = τ ∇v m + p m , the previous equation becomes:

w m+1 1 + 1 γ 2 + |w m+1 | 2 = C m (41)
(41) is easily solved with a fixed point iteration. Indeed we have the following result:

Proposition 4 Let X the Euclidean space R N ×N , and Y = X × X. Consider the sequence x 0 = w m :

x k+1 = C m γ 2 + |x k | 2 1 + γ 2 + |x k | 2 (42)
Then

x k → w m+1 in Y as k → +∞.
The proof of this result will be detailed in Section 4.3. Bermùdez and Moreno algorithm has already been used for smoothed total variation based restoration in [START_REF] Almansa | A TV based restoration model with local constraints[END_REF]. The authors of [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] use a different approach than in this paper.

To solve [START_REF] Ng | On semismooth Newton methods for total variation minimization[END_REF], they take the square of both sides of (41), and they use a Newton method to compute |w m+1 |. They then compute w m+1 with [START_REF] Ng | On semismooth Newton methods for total variation minimization[END_REF]. But with such an approach, the authors of [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] report poor numerical results. We also tried this approach, and we have seen the same poor results as in [START_REF] Almansa | A TV based restoration model with local constraints[END_REF]. We therefore advocate the use of the fixed point algorithm proposed here to solve [START_REF] Ng | On semismooth Newton methods for total variation minimization[END_REF], which we prove to converge without further assumption (notice that another alternative would be to solve directly [START_REF] Ng | On semismooth Newton methods for total variation minimization[END_REF] with Newton method). The final scheme to solve (30) is thus:

           v m = f µ + div p m w m+1 =   1 + 1 r β 2 τ 2 µ 2 +|w m+1 | 2   -1 (τ ∇v m + p m ) p m+1 = τ ∇v m + p m -w m+1 (43)
The second equation is solved with a fixed point iteration [START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF]. We will see that in practice, a single iteration is enough, and thus the second line of (43) reduces to:

w m+1 =   1 + 1 β 2 τ 2 µ 2 + |w m | 2   -1 (τ ∇v m + p m ) (44)
Applying Theorem 1, we have the following convergence result:

Proposition 5 Let X the Euclidean space R N ×N , and Y = X × X. If τ < 1 4 , then the sequence (v m , w m , p m ) defined by scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] is such that v m → v in X, w m → w in Y , and p m → p in Y with µv solution of (30).

Interpretation of scheme (43)

One first needs to remember that we are interested in solving problem [START_REF] Duval | Projected gradient based color image decomposition[END_REF]. Using the change of notation v = u/µ, solving (30) is equivalent to solving:

inf v β 2 µ 2 + |∇v| 2 dx + 1 2 f µ -v 2 (45)
The associated Euler-equation is:

0 = v - f µ -div   ∇v β 2 µ 2 + |∇v| 2   (46)
The most classical methods to solve this equation are the fixed step gradient descent as in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and the quasi-Newton method (which can be seen also as semi-quadratic regularization) as for instance in [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Dobson | Convergence of an iterative method for total variation denoising[END_REF][START_REF] Acar | Analysis of total variation penalty methods for ill-posed problems[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Chan | On the convergence of the lagged diffusity fixed point method in total variation image restoration[END_REF][START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF]. The idea of the quasi-Newton method is to linearize the non-linear term in the above equation, and to consider an iterative scheme of the type:

0 = v m+1 - f µ -div   ∇v m+1 β 2 µ 2 + |∇v m | 2   (47)
Here, we propose a different iterative scheme to solve [START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF].

0 = v m - f µ -div p m (48)
with

p m = z m β 2 µ 2 + |z m | 2 (49)
In the limit, we would like to have z m → ∇v. To update p m , we use the following equation:

p m+1 = p m + τ ∇v m -z m+1 (50)
If (p m ) converges, then v m → v with [START_REF] Setzer | Split Bregman Algorithm, Douglas-Rachford Splitting and frame shrinkage[END_REF], and z m → ∇v with (50) as m → +∞. The system of equations ( 48)-( 49)-( 50) can be rewritten into:

           v m = f µ + div p m z m+1   τ + 1 r β 2 µ 2 +|z m+1 | 2   = τ ∇v m + p m p m+1 = p m + τ (∇v m -z m+1 ) (51)
If we make the change of variable z m = w m /τ , then scheme (51) is exactly [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF], i.e. Bermùdez-Moreno algorithm for solving problem [START_REF] Duval | Projected gradient based color image decomposition[END_REF].

Convergence of the fixed point iteration

In this section, we detail the proof of Proposition 4. The proof relies on Weizfeld method [START_REF] Weisfeld | Sur le point pour lequel la somme des distances de points donnés est minimum[END_REF][START_REF] Chan | On the convergence of the lagged diffusity fixed point method in total variation image restoration[END_REF][START_REF] Facciolo | Irregular to regular sampling, denoising and deconvolution[END_REF]. We adopt here the presentation of [START_REF] Chan | On the convergence of the lagged diffusity fixed point method in total variation image restoration[END_REF] for Weizfeld method. Let us first introduce some notations. We consider the following functional:

F (u) = 1 2 u -C 2 + (γ 2 + |u| 2 ) 1/4 2 (52)
We have:

∇F (u) = u -C + u γ 2 + |u| 2 (53)
Let us define:

A(u) = I + I γ 2 + |u| 2 (54)
Notice that u → A(u) is continuous, and that λ min (A(u)) ≥ 1, where λ min (M ) stands for the smallest eigenvalue of M . Let us finally define:

G(v, u) = F (u) + v -u, ∇F (u) + 1 2 v -u, A(u)(v -u) (55)
Notice that G consists in a linearization of F . In fact, G defines a general Weizfeld method for the problem:

inf u F (u) (56)
Notice that since F is strictly convex and coercive, there exists a unique u solution of (56), and u is the solution of:

∇F (u) = u 1 + 1 γ 2 + |u| 2 -C = 0 (57)
We now define the iteration of Weizfeld method:

u m+1 = argmin v G(v, u m ) ( 58
)
Since G is strictly convex and coercive, there exists a unique u m+1 solution of (58). It satisfies the Euler-Lagrange equation:

∇F (u m ) + A(u m )(u m+1 -u m ) = 0 (59)
i.e.:

u m+1 1 + 1 γ 2 + |u m | 2 = C (60)
which is precisely iteration [START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF].

Proposition 6 If u is fixed, then for all v we have:

G(v, u) -F (v) ≥ 0.

Proof.

A standard computation leads to:

G(v, u) -F (v) = u -v, -1 2
u + v γ 2 + |u| 2 + γ 2 + |u| 2 -γ 2 + |v| 2 dx = -1 2 |u| 2 -|v| 2 γ 2 + |u| 2 dx + γ 2 + |u| 2 -γ 2 + |v| 2 dx
Using the notation a = γ 2 + |u| 2 and b = γ 2 + |v| 2 , we get:

G(v, u) -F (v) = -1 2 a -b a + a -b dx = (a -b) 2 2a dx ≥ 0 (61)
The following lemma holds:

Lemma 1
We have for all m:

F (u m+1 ) ≤ F (u m ) ( 62
)
and lim m→+∞ u m+1 -u m = 0 (63)

Proof. From Proposition 6, we have

F (u m+1 ) ≤ G(u m+1 , u m ). But from (58), we get G(u m+1 , u m ) ≤ G(u m , u m ) = F (u m
). We thus deduce inequality (62).

We now concentrate on proving (63). From Proposition 6, we have:

F (u m+1 ) ≤ G(u m+1 , u m ) = F (u m ) + u m+1 -u m , ∇F (u m ) + 1 2 u m+1 -u m , A(u m )(u m+1 -u m ) = F (u m ) - 1 2 u m+1 -u m , A(u m )(u m+1 -u m )
where we have used Equation (59). We thus deduce that (since λ min (A(u)) ≥ 1):

1 2 u m+1 -u m 2 ≤ 1 2 u m+1 -u m , A(u m )(u m+1 -u m ) ≤ F (u m ) -F (u m+1 )
We finally get that:

u m+1 -u m ≤ 2(F (u m ) -F (u m+1 )) (64)
We have just seen before that F (u m ) is a positive, monotone decreasing sequence. Hence F (u m ) is a convergent sequence, and in particular F (u m ) -F (u m+1 ) → 0, which concludes the proof.

We are now in position to prove the convergence of the fixed point iteration as stated in Proposition 4: Proof. From (60), one sees that u m is uniformly bounded. Therefore, up to a subsequence, u m converges to some v. Moreover, from Lemma 1, we see that u m+1 also converges to v. Passing to the limit in (60), we see that v = u where u is the unique minimizer of (56). We conclude that the whole sequence u m goes to u.

We end this section by stating a result about the convergence rate of the fixed point algorithm [START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF]. We denote by ũ the solution of Problem (56). We use the following notations:

γ m = G(ũ, u m ) -F (ũ) 1 2 ũ -u m , A(u m )(ũ -u m ) (65)
and

η = 1 -λ min (A(ũ) -1 ∇ 2 F (ũ)) (66) Proposition 7 1. F (u m+1 ) -F (ũ) ≤ γ m (F (u m ) -F (ũ)).
2. η < 1 and 0 ≤ γ m ≤ η, for m sufficiently large. In particular, F (u m ) has a linear convergence rate of at most η.

3. u m is r-linearly convergent with a convergent rate of at most √ η.

Proof. We refer the interested reader to the proof of Theorem 6.1 in [START_REF] Chan | On the convergence of the lagged diffusity fixed point method in total variation image restoration[END_REF].

Image deconvolution

In this section, we consider the problem of image deconvolution. We explain how Bermùdez-Moreno algorithm [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF] can be applied to this problem, by using the iterative approach of [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Bect | A l1-unified variational framework for image restoration[END_REF]. In all the previous sections, we have considered the denoising problem:

1 2µ u -f 2 + φ β (u) (67)
with the convention that φ 0 (u) = Ω |Du|. As probably noticed by the reader, Bermùdez-Moreno scheme can be applied to functional of the type:

1 2µ Au -f 2 + φ β (u) (68)
provided that A is an easily invertible operator. However, in the case of image deblurring, the operator A is ill-posed, and we can therefore not apply Bermùdez-Moreno scheme directly. A possible alternative is to use an iterative approach as proposed in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] or [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF]. This type of approach has grown very popular and is now widely used to handle sparsity constraints [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. Here we use the presentation of [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF]. The trick of the method lies in the following result:

Proposition 8 Let B a linear positive symmetric invertible operator with B < 1.

Let C = B(I -B) -1 .
Then, for all u we have:

Bu, u = inf w u -w 2 + Cw, w (69)
Moreover, the minimum is reached for

w = (I + C) -1 (u) = (I -B)(u) (70)
Here, we choose ν > 0 such that νA * A < 1, and we set B = νA * A. Let us set:

H(u, w) = 1 2µν u -w 2 + Cw, w + 1 2µ f 2 -2 Au, f (71)
Using Proposition 8, it is easy to see that

1 2µ f -Au 2 = inf w H(u, w) (72)
Let us now define

F (u, w) = H(u, w) + φ β (u) (73)
Let us consider the following algorithm:

w n = (I -A * A)(u n ) u n+1 = argmin u 1 2µν w n + νA * f -u 2 + φ β (u) (74)
Setting v n = w n + νA * f , it can be written:

v n = u n + νA * (f -Au n ) u n+1 = argmin u 1 2µν v n -u 2 + φ β (u) (75)
The following convergence result is shown in [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF]:

Proposition 9 Let X the Euclidean space R N ×N . The sequence (u n , v n ) defined by scheme (75) is such that (u n , v n ) → (u, v) in X × X with (u, v -νA * f ) minimizer of (73).
In practice, to solve the second line of (75) we use scheme [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] if β = 0 and scheme (33) if β > 0. Notice that (75) can also be interpreted as a Forward-Backward splitting algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] to solve problem (68). Numerical results: (75) is very easy to implement. Notice that in such an iterative approach, one of the key point is to be able to solve each iteration efficiently, which is the case with scheme ( 24) or [START_REF] Fadili | Monotone operator splitting for fast sparse solutions of inverse problems[END_REF]. We show some numerical results on Figure 3. As expected, total variation regularization deconvolution gives sharper edges, whereas smoothed total variation based deconvolution preserves better the textures.

Nesterov algorithms

In the previous section, we have introduced first order numerical schemes [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] or [START_REF] Fadili | Monotone operator splitting for fast sparse solutions of inverse problems[END_REF] to solve image restoration problems. To see how efficient they are, we have decided to compare them with state of the art first order numerical schemes. It has been shown in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] that Nesterov schemes are very efficient to solve image restoration problems: they beat all the other existing first order algorithms. These schemes were recently introduced in [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF][START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF], and they have proved to be a significant improvement in convex optimization. Nevertheless, except in the work by P. Weiss et al. [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF], Nesterov schemes have not been applied yet in the image processing community. Notice that contrary to the previous sections, Section 6 is not related to Bermùdez and Moreno framework, and thus it can be read independently. Let us emphasize again that we present Nesterov algorithms because they are optimal first order schemes, and because we want to compare Bermùdez and Moreno approach with state of the art first order schemes.

We first recall Nesterov schemes in Section 6.1. Motivated by [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] and our first numerical results, we have decided to implement some improvements of Nesterov schemes recently introduced by Y. Nesterov in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. A first variant is presented in Section 6.2 and a second one in Section 6.3.

Nesterov schemes

In [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF][START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF], Y. Nesterov proposes efficient schemes to minimize functionals such as [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF] or [START_REF] Duval | Projected gradient based color image decomposition[END_REF]. We follow here the presentation of [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF]. We consider the following minimization problem:

inf u∈Q E(u) (76)
where E is a convex Lipschitz differentiable function, and Q a convex closed set. We denote by ũ a solution of (76). For this type of problem, it can be shown [START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF] (Theorem 2.1.7 page 61) that no algorithm (only using the values and gradients of E) has a better rate of convergence than O 1 k 2 uniformly on all problems of the form (76) (k is the number of iterations of the algorithm).

In the framework developped by Nesterov, the convergence rate is in term of objective function, and not in term of distance to the minimizer. For instance, a convergence rate of O 1

k 2 for problem (76) means that |E(u k ) -E(ũ)| ≤ C k 2
, where ũ is the solution of problem (76), and u k the approximation of ũ at iteration k. Of course, without further assumption, it gives no information on the convergence rate of u k to ũ. To get such a piece of information, a coercivity hypothesis is needed for the functional E. Nevertheless, as shown in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] and in the present paper, the convergence rate obtained with Nesterov's results are in accordance with the behavior of the algorithms (in the sense that if a scheme is supposed to converge as O 1 k 2 and a second one as O 1 k , it is indeed numerically observed that the first scheme converges faster to the solution).

Let us notice that the constant hidden in the convergence rate in Nesterov's theory is always proportional to L × u 0 -ũ 2 , where L is the Lipschitz constant of ∇E, and u 0 the intial guess for the minimizer ũ.

In [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] is given an O 1 k 2 algorithm for solving problem (76) which we detail here-after (it is thus optimal in the sense of Nesterov).

Let . be a norm and d a convex function such that there exists σ > 0 and x 0 in Q satisfying for all x the inequality:

d(x) ≥ σ 2 x -x 0 2 . 1. Set k = 0, v 0 = 0, x 0 ∈ Q, L Lipschitz constant of ∇E.
2. Set k = k + 1, and compute η k = ∇E(x k ).

Set y

k = argmin y∈Q η k , y -x k + 1 2 L y -x k 2 . 4. Set v k = v k-1 + k+1 2 η k . 5. Set z k = argmin y∈Q L σ d(x) + v k , z . 6. Set x k+1 = 2 k+3 z k + k+1 k+3 y k .
Proposition 10 [39] The previous algorithm ensures that :

0 ≤ E(y k ) -E(ũ) ≤ 4Ld(ũ) σ(k + 1)(k + 2) ( 77
)
The idea behind Nesterov's scheme is similar to the one of the conjugate gradient algorithm [START_REF] Ciarlet | Introduction à l'analyse numérique matricielle et à l'optimisation[END_REF]: the direction of descent is at step k + 1 is computed by taken into account the information of the complete sequence of gradient (∇E(x 0 ), . . . , ∇E(x k )), and not only ∇E(x k ).

Primal Nesterov algorithm: For β > 0, we remind the reader that we set φ β (u) = β 2 + |∇u| 2 dx. Nesterov algorithm can be used to solve the following problem:

inf u∈Kα φ β (f + u) (78)
where

K α = x ∈ L 2 / x 2 ≤ α .
This problem is equivalent to problem (30) (see [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF] for a complete analysis). The advantage of formulation (78) is that Nesterov's scheme can directly be applied. See [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] (Algorithm 2 page 12) for a detailed implementation of this algorithm. We will refer to it as the primal Nesterov algorithm. We just give here the sketch of the algorithm (P Kα is the orthogonal projection onto K α ):

1. Set k = 0, v 0 = 0, x 0 = 0, L = div 2 /β = 8/β. 2. Set k = k + 1, and compute η k = -div ∇(x k +f ) √ β 2 +|∇(x k +f )| 2 . 3. Set y k = P Kα (x k -η k /L), with K α = x ∈ L 2 / x 2 ≤ α . 4. Set v k = v k-1 + k+1 2 η k . 5. Set z k = P Kα (-v k /L). 6. Set x k+1 = 2 k+3 z k + k+1 k+3 y k .
7. The output of the algorithm is: u = y lim + f .

Dual Nesterov algorithm:

Of course, due to the non-differentiability in zero of the total variation, Nesterov scheme cannot be applied directly to problem [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]. The basic idea is to apply Nesterov's scheme to the dual version of [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF], that is to: inf f -u∈Gµ

1 2 u 2
, where G µ is given by ( 22), i.e.:

inf q∈K E(q) ( 79
)
where

E(q) = 1 2 f -µdiv q 2 and K = x ∈ L 2 × L 2 / x ≤ 1 .
If we denote by ũ the solution of [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF], and by q the solution of (79), we have ũ = f -µdiv q.

See [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] (Algorithm 3 page 20) for a detailed implementation of this algorithm. We will refer to it as the dual Nesterov algorithm. We just give here the sketch of the algorithm (P K is the orthogonal projection onto K):

1. Set k = 0, v 0 = 0, x 0 = 0, L = µ div 2 = 8µ. 2. Set k = k + 1, and compute η k = ∇ f -µdiv (x k ) . 3. Set y k = P K (x k -η k /L), with K = x ∈ L 2 × L 2 / x ≤ 1 . 4. Set v k = v k-1 + k+1 2 η k . 5. Set z k = P K (-v k /L). 6. Set x k+1 = 2 k+3 z k + k+1 k+3 y k .
7. The output of the algorithm is: u = f -µdiv (y lim ).

Notice that in the dual Nesterov algorithm, the set K is included in L 2 × L 2 ; whereas in the case of the primal Nesterov algorithm, the set K α is embedded in L 2 .

In [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF], very good numerical results are reported both for the primal and the dual Nesterov algorithms (much better than steepest gradient descent for instance). We have therefore decided to use them as reference in the comparisons presented here-after. We remind the reader that for the projected gradient scheme [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] for solving [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF], we mention in Subsection 3.4

Accelerated Nesterov algorithm

In [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF], Y. Nesterov proposes a way to speed up the minimization algorithms introduced in [START_REF] Nesterov | Introductory lectures on convex optimization: a basic course[END_REF]. The idea is to improve the estimation of the Lipschitz constant of the functional to minimize (in view of equation ( 77)). In this subsection, we show how it can be used for image restoration. Consider the general minimization problem

inf u E(u) + ψ(u) (80)
We set φ(u) = E(u) + ψ(u), and:

ψ(u) = χ Q (u) = 0 if u ∈ Q +∞ otherwise. (81)
Problem ( 80) is therefore the same as (76). As previously, E is a convex Lipschitz differentiable function, and Q a convex closed set. We denote by ũ a solution of (80). We set:

T L (y) = argmin x∈Q m L (y, x) (82) with m L (y, x) = E(y) + ∇E(y), x -y + L 2 x -y 2 + ψ(x) (83)
Moreover, it is shown in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] that

φ ′ (T L (y)) = L(y -T L (y)) + ∇E(T L (y)) -∇E(y) (84)
In [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] is given an efficient algorithm for solving problem (80):

• Set k = 0, A 0 = 0, v 0 = 0, x 0 ∈ Q, L 0 = L Lipschitz constant of ∇E, ψ 0 (x) = 1 2 x -x 0 2 . Set γ u > 1 and γ d ≥ 1. • Set L = L k . REPEAT: Set a = 1+ √ 1+4A k L 2L . Set y = A k x k +av k
A k +a , and compute T L (y).

If: φ ′ (T L (y)), y -T L (y) < 1 2L φ ′ (T L (y)) 2 2 , then L = γ u L. UNTIL: φ ′ (T L (y)), y -T L (y) ≥ 1 2L φ ′ (T L (y)) 2 DEFINE y k = y, M k = L, a k+1 = a, A k+1 = A k + a k+1 , L k+1 = M k /γ d , x k+1 = T M k (y k ), ψ k+1 (x) = ψ k (x) + a k+1 E(x k+1 ) + ∇E(x k+1 ), x -x k+1 + ψ(x) , v k+1 = argmin x ψ k+1 (x).
Output: the output of the algorithm is u = x lim .

The following convergence result is shown in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]:

Proposition 11 [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] Let L E the Lipschitz constant of ∇E. Assume that 0 < L 0 ≤ L E . Then the previous algorithm ensures that :

0 ≤ φ(x k ) -φ(ũ) ≤ 4γ u L E ũ -x 0 2 k 2 (85)
where we recall that φ(u) = E(u) + ψ(u).

To apply this new algorithm, the only points to check are how to solve (82) and how to compute v k . This is explained by the two following lemmas.

Lemma 2

The solution of problem ( 82) is given by:

T L (y) = P Q y - 1 L ∇E(y) (86)
with P Q orthogonal projection onto Q.

Proof. It is easy to see that:

m L (y, x) = C(y) + L 2 x -y - 1 L ∇E(y) 2 2 + ψ(x) (87)
where C(y) is a function depending only on y. The result of the lemma follows from the fact that ψ = χ Q .

Lemma 3 v k = argmin x ψ k (x) is given by

v k = P Q x 0 - k p=1 a p ∇E(x p ) ( 88
)
with P Q orthogonal projection onto Q.

Proof. Remembering that ψ 0 (x) = 1 2 x -x 0 2 , it is easy to see that:

ψ k (x) = C(k) + k p=1 a p ψ(x) + 1 2 x -x 0 + k p=1 a p ∇E(x p ) 2 2 (89)
where C(k) is a function depending only on k. The result of the lemma follows from the fact that ψ = χ Q .

In practice, as proposed in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF], we use

γ u = γ d = 2.
Application to problem (78): The above algorithm can directly be applied to (78), with

E(u) = φ β (u + f ) = β 2 + |∇(f + u)| 2 dx and ψ(u) = χ Kα (u) (90)
where

K α = x ∈ L 2 / x 2 ≤ α . Of course, one has: ∇E(u) = -div ∇(f +u) √ β 2 +|∇(f +u)| 2 .
One just has to set x 0 = 0, L = div 2 /β = 8/β. The solution is given by f + x lim . Notice that here, the projection onto Q = K α is straightforward: P Kα (x) = αx max{α, x 2} . We will refer to this algorithm as the accelerated primal Nesterov algorithm.

Application to problem [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]: The basic idea is to apply the accelerated Nesterov scheme to the dual version of [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF], that is to (79)), i.e.: inf q E(q) + ψ(q) (91)

with E(q) = 1 2 f -µdiv q 2 2 and ψ(q) = χ K (q) with K = {g ∈ L 2 × L 2 , g 2 1 + g 2 2 ≤ 1}. We therefore have:

∇E(q) = ∇(f -µdiv q).
One just has to set u 0 = 0, L = µ div 2 = 8µ. The solution is given by f -µdiv x lim . Notice that here, the projection onto Q = K is straightforward:

P K (x 1 , x 2 ) = 1 max{1, x } (x 1 , x 2 )), with x = (x 1 , x 2 ) and x = x 2 1 + x 2
2 . We will refer to this algorithm as the accelerated dual Nesterov algorithm.

Variant for the accelerated Nesterov algorithm

In [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF], Y. Nesterov proposes in fact a more general algorithm than the one we have presented in Section 6.2. We show here how it can be used to solve image restoration problems. We still consider the general minimization problem inf

u E(u) + ψ(u) (92)
But this time ψ is assumed to be a strongly convex function with parameter µ ψ > 0: in the case when ψ is C 2 , it means that the smallest eigenvalue of ∇ 2 ψ is µ ψ > 0.

We set φ(u) = E(u) + ψ(u). As previously, E is a convex Lipschitz differentiable function We denote by ũ a solution of (92). Moreover, in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] is given an efficient algorithm for solving problem (92): this is exactly the algorithm presented in Section 6.2, the only difference being that in the step REPEAT, instead of setting

a = 1+ √ 1+4A k L 2L
, we set:

a = b + √ b 2 + 4A k b 2 with b = 1 + µ ψ A k L (93)
The following convergence result is shown in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]:

Proposition 12 [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] Let L E the Lipschitz constant of E, and µ ψ the convexity parameter of ψ. Assume that 0 < L 0 ≤ L E . Then the previous algorithm ensures that (85) still holds. Moreover, we also have:

0 ≤ φ(x k ) -φ(ũ) ≤ γ u L E ũ -x 0 2 1 + µ ψ 8γ u L E -2(k-1) (94)
Notice that (84) still holds in this case. To apply this new algorithm, the only points to check are how to solve (82) and how to compute v k . We particularize the problem, and we consider the restoration problem [START_REF] Duval | Projected gradient based color image decomposition[END_REF], i.e. in (92) we take:

E(u) = φ β (u + f ) = β 2 + |∇(f + u)| 2 dx and ψ(u) = 1 2µ u 2 (95)
Notice that we have:

L E = div 2 /β = 8/β and µ ψ = 1 µ (96)
The two following lemmas hold.

Lemma 4

The solution of problem (82) is given by:

T L (y) = Ly -∇E(y) L + 1 µ (97)
Proof. It is easy to see that:

∇ x (m L (y, x)) = ∇E(y) + L(x -y) + x µ (98)
Lemma 5 v k = argmin x ψ k (x) is given by:

v k = 1 1 + P k p=1 a p µ x 0 - k p=1 a p ∇E(x p ) (99)
Proof. Remembering that ψ 0 (x) = 1 2 x -x 0 2 , it is easy to see that:

ψ k (x) = 1 2 x -x 0 2 + k p=1 a p ψ(x) + k p=1 a p (E(x p ) + ∇E(x p ), x -x p ) (100)
In the next section, we will refer to this algorithm as the variant of the accelerated primal Nesterov algorithm. In practice, we take x 0 = 0, γ u = 2 and γ d = 2.

Numerical examples

In this section, we present some numerical examples with the schemes introduced in this paper. See also [START_REF] Aujol | Some algorithms for total variation based image restoration[END_REF] for more numerical results. Notice that all the experiments presented in this paper were run with Matlab, on a laptop with a processor at 2 GHz and 2 Gb of RAM. In all the presented algorithms, the cost of one iteration of the algorithm is proportional to the size of the image. This cost is around 0.03 second for a 256×256 image with either the fixed point algorithm [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF], the projected gradient algorithm [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF], or Chambolle projection algorithm [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. The primal or dual Nesterov algorithms (Section 6.1) have a cost per iteration which is twice higher. This cost per iteration is between 8 and 10 times higher with the different variants of Nesterov algorithms (Section 6.2 and 6.3).

In Section 7.1, we consider the case of smoothed total variation regularization, and in Section 7.2 we are interested in total variation regularization.

Smoothed total variation regularization

Restoration results obtained with the new scheme (43): We illustrate here the efficiency of scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] (based on Bermùdez and Moreno framework) to solve problem [START_REF] Duval | Projected gradient based color image decomposition[END_REF]. This new scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] has the advantages of being simple and stable. Moreover, it seems quite fast (less then 2 seconds for a 256×256 image to get a normalized L 2 error smaller than 0.5 with β = 1, the images having their values in the range [0; 255]). On Figure 4, we show the restoration results we get on the noisy images of Figure 1. The curvature parameter β of ( 43) is fixed to 10. As expected, the textures are better preserved with this model than with total variation regularization (compare with Figure 2), but the edges are not as sharp.

Influence of the number of iterations in the fixed point loop [START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF]: We now want to see the speed of convergence of [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF], and how it depends on the number of iterations in loop [START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF], and on the parameter β. For different values of β, we compute an ideal image by running 10 000 iterations of [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] with 500 iterations for the fixed point [START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF]. We can then compute at each iteration the L 2 error between a computed image with [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] and the target ideal image. On Figure 5, we show the behavior of the algorithm with respect to the number of iterations for the fixed point iteration, for different values of β. Clearly, it shows that 1 iteration is a very good choice: this will be our choice until the end of the paper. It is also clear that the convergence of ( 43) is much faster for large values of β.

Notice that there exist some theoretical convergence results about iterative schemes using an inner fixed-point loop. For instance, it is shown in [START_REF] Bonesky | A generalized conditional gradient method for nonlinear operator equations with sparsity constraints[END_REF] that one fixed-point iteration of the method of [START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF] is enough to get convergence (for the algorithm of [START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF]). This confirms the numerical observation made in this paper that in the inner fixed-point loop, one iteration may be sufficient to get convergence.

Comparisons with Nesterov schemes: On Figure 6, we compare our new algorithm [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] with the primal Nesterov algorithm (Section 6.1), the accelerated primal Nesterov algorithm (Section 6.2), and the variant of the accelerated primal Nesterov algorithm (Section 6.3). Notice that since algorithm [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] uses a fixed point iteration, On graphs (b) and (d), the primal Nesterov algorithm is not shown. Notice that the time for 1 iteration of the primal Nesterov algorithm is around twice the time for 1 iteration of the fixed point based algorithm [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF]. The accelerated primal Nesterov algorithm and its variant are between 4 and 5 times slower per iteration than the primal Nesterov algorithm. To get a fast approximation, the fixed point based algorithm [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] seems to be the best choice (the accuracy is good enough for image restoration). To get a highly accurate solution, the variant of the accelerated Nesterov scheme seems to be the most efficient.

we refer to it as fixed point algorithm in the caption of Figure 6. The convergence speed of these three last algorithms depends on the Lipschitz constant of the energy to minimize: the smaller this constant, the faster the method. It thus means here the larger β, the faster the method. Notice that here the images we consider have their range in [0, 255] (while for instance in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] the images are normalized in [0, 1]: this has some impact on the values β proposed here).

To make comparisons, we compute the L 2 -norm of the difference between the original image and the ideal image (obtained by running [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] with 10 000 iterations). We then set this L 2 -norm as the constraint in the primal Nesterov algorithm and the accelerated primal Nesterov algorithm . It is to be noticed that such a choice makes a small bias in favor of our scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF]. However, the obtained results are sufficiently convincing to forget this bias.

It can be seen that, the larger β, the faster the algorithms. For large values of β, all the algorithms are fast. However, when β decreases to zero, then scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] seems to bring a significant increase in speed of convergence towards a good approximation. It seems indeed that [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] can lead to a good approximation of the minimizer with few iterations. However, when one is interested in getting a very accurate solution, then the variant of the accelerated primal Nesterov algorithm seems to be the best choice. This is in accordance with the result of Proposition 12. Notice that the cost per iteration of scheme ( 43) is twice lower than for the primal Nesterov algorithm , while the accelerated primal Nesterov algorithm and its variant are between 4 and 5 times slower per iteration than the primal Nesterov algorithm.

Notice that the quality of the restored image obtained with scheme (43) after a few iterations (10 iterations for β = 25, 20 iterations for β = 10, 80 iterations for β = 1, 200 iterations for β = 0.1) is visually very good: the normalized L 2 error is then smaller than 0.3. For a restoration purpose, there is no need for the accuracy of the variant of the accelerated primal Nesterov algorithm. It is more important to have a fast approximation than a slow and very accurate solution.

Total variation regularization

In this section, we consider problem [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]. We want to compare five different algorithms. The first one is the projection algorithm of [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]: we refer to it as Chambolle projection algorithm. We use τ = 0.249 in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. The second one is the modification of this algorithm as proposed in [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF], and which we proved to be Bermùdez-Moreno algorithm [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] in the case of problem [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF]: since it is an adaptation of Uzawa method [START_REF] Ciarlet | Introduction à l'analyse numérique matricielle et à l'optimisation[END_REF] to problem [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF], we refer to it as Uzawa algorithm. We use τ = 0.249 in [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. The third algorithm we use here is the dual Nesterov algorithm presented in Section 6.1, as proposed in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF]. Motivated by the results of [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] and [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF], we use it as the reference algorithm. The fourth algorithm we use here is the accelerated dual Nesterov algorithm of [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] presented in Section 6.2. The fifth algorithm we use is our new scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF]. Since it uses a fixed point algorithm, we refer to it as fixed point method.

For a given image and a given regularization parameter µ, a reference ideal image is computed by running 10 000 iterations with the dual Nesterov algorithm. Here, the bias will therefore be in favor of the dual Nesterov algorithm. However, we think that the results are convincing enough to forget this bias.

A convergence speed result is presented on Figure 7: we give the L 2 -norm of ũ -u n , where u n is the computed image at iteration n, and ũ the ideal image to obtain. As can be seen on Figure 7, the dual Nesterov algorithm is faster then Uzawa algorithm, which is itself faster than Chambolle projection algorithm. The accelerated dual Nesterov algorithm seems to be the best choice to get a highly accurate solution (50 iterations to get a normalized L 2 error of 0.3). However, 1 iteration with the accelerated dual Nesterov algorithm is around 4 times slower than with the dual Nesterov algorithm: the dual Nesterov algorithm seems thus a good compromise when one is interested in getting a very good approximation. Nevertheless, 1 iteration with the dual Nesterov algorithm is around 2 times slower than with Uzawa, Chambolle, or scheme (43) (while all three have the same computation time per iteration). For typical image restoration problems (with Gaussian noise), ( 24) seems 30% faster then [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] (for instance, it takes 70 iterations for [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] to get a normalized L 2 error of 1 while it takes 110 iterations for [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] to get the same accuracy). Algorithm [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] seems to be a good alternative when one is only interested in getting an approximation with a small number of iterations; for instance, scheme [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] is the fastest (in term of computation time) to get a normalized L 2 error of 2 (around 30 iterations).

In [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF], the authors explain that the dual Nesterov algorithm is much faster then the projected gradient method (24) (Uzawa algorithm) for total variation regularization. We confirm that it is indeed much faster when one is interested in computing an accurate solution. Notice also that in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF], the comparison criterion used is the is itself around 2 times slower than with Uzawa, Chambolle, or scheme (43) (while all three have the same computation time per iteration). To get a highly accurate solution, the accelerated dual Nesterov algorithm seems to be the best choice. However, the dual Nesterov algorithm seems to be the best compromise when one is only interested in getting a good approximation (which is the case for image restoration).

value of the total variation of the computed image. This is indeed the quantity which is controlled in Nesterov's approach for solving [START_REF] Bioucas-Dias | Thresholding algorithms for image restoration[END_REF] (see Proposition 10). Here, the criterion is the L 2 difference of the computed solution for some iteration with the ideal solution. Figure 7 is surely in favor of the approach developed in [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF]. However, the difference during the first iterations is not that large, and thus the projected gradient algorithm (24) (Uzawa algorithm) can still be considered as a good method when one is only interested in getting a fast approximation of the solution.

Dual Nesterov algorithm for solving [START_REF] Duval | Projected gradient based color image decomposition[END_REF]: In view of Figure 7, one should be tempted to use the dual Nesterov algorithm for solving [START_REF] Duval | Projected gradient based color image decomposition[END_REF]. It is easy to compute the dual problem. If we denote by ũ the solution of ( 30), then we have ũ = f -µdiv p with p solution of:

inf p∈K 1 2µ µdiv p -f 2 -β 1 -|p| 2 dx ( 101
)
where K = p ∈ L 2 × L 2 / p ∞ ≤ 1 . However, the gradient of the functional in (101) is not Lipschitz, and we therefore cannot use directly the dual Nesterov algorithm.

A Proof of convergence of Bermùdez-Moreno algorithm

In this section, we follow [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF] and [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] (chapter II.3). Our goal is to give the reader some intuition on why the result of Theorem 1 holds. We remind the reader that we use the notations: H λ = I-L λ λ , with L λ = (I + λH) -1 and H = ∂φ with φ proper convex lower semi continuous function. We will use the next lemma:

Lemma 6 1 λ 2 L λ (v 1 ) -L λ (v 2 ) 2 + H λ (v 1 ) -H λ (v 2 ) 2 ≤ 1 λ 2 v 1 -v 2 2 (102)
Proof. This is an immediate consequence of definitions [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF]. Problem ( 3) is related to: ∀z , Au, z -u + ψ(z) -ψ(u) ≥ g, z -u

The relation is given by the next lemma (whose proof is straightforward (see [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] proposition 2.2 page 37)):

Lemma 7 u is solution of (103) if and only if u is solution of (3).

We remind the reader that B∂φ(B * u) = ∂ψ(u) (see [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] proposition 5.7 page 27). Problem (103) is related to the subdifferential inclusion:

g -Au ∈ B∂φ(B * u)

The relation is given by the next proposition:

Proposition 13 u is solution of (104) if and only if u is solution of (103).

Proof. The fact that u solution of (104) implies that u solution of (103) is a direct consequence of the definition of the subdifferential of a convex function [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF]. The reciprocal result is more complicated, and we refer the reader to chapter II.3 of [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] for a detailed proof.

We will make use of the next lemma (Lemma 2.1 in [START_REF] Bermùdez | Duality methods for solving variational inequalities[END_REF]):

Lemma 8 H maximal monotone operator. Then the two following conditions are equivalent:

(i) y ∈ H(v) (ii) y = H λ (v + λy)
An immediate consequence of the previous lemma is the following result:

Proposition 14 u is a solution of (104) if and only if (u, y) is a solution of:

Au = g -By y = H λ (B * u + λy) (105)
We are now in position to prove Theorem 1.

Proof. From (102), we get: We now deduce from (106) that:

1 λ 2 L λ (B * u + λy) -L λ (B * u m + λy m ) 2 E + y -y m+1 2 E ≤ 1 λ 1 λ -2α B * 2 B * (u -u m ) 2 + y -y m 2 E ( 109
)
We eventually get that, since 0 < 1 λ < 2α B * 2 , as long as u m = u: y -y m+1 E < y -y m E . We deduce that y -y m 2 E is a convergent sequence in R. Thus passing to the limit in (109), we get: lim m→+∞ B * (u-u m ) E = 0. Using (108), we eventually get that u m → u.

There remains to prove that y m also converges. We first remark that now, passing to the limit in (109), we get: L λ (B * u m + λy m ) → L λ (B * u + λy). But since L λ = I -λH λ , we get with the second line of (105) that: L λ (B * u + λy) = B * u. From the second line of (9), we get:

y m+1 = H λ (B * u m + λy m ) = y m + 1 λ (B * u m -L λ (B * u m + λy m )) (110)
Passing to the limit, we eventually get that: lim m→+∞ {y m+1 -y m } = 0. Now we can conclude that y m ⇀ y in E weak, since the application

v ∈ E → H λ (B * u(v) + λv) (111)
with u(v) solution of: Au = g -Bv, is non expansive (see [START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] Corollary 4 p.199).

Figure 1 :

 1 Figure 1: The classical Lenna and cameraman image, and their noisy version (additive zero mean Gaussian noise with standard deviation σ = 20).

Figure 3 :

 3 Figure 3: (a) and (b): degraded images (convolved with a Gaussian kernel with standard deviation η = 5, and then degraded by a zero mean Gaussian noise with standard deviation σ = 10), the original images are on top row of Figure 1; (c) and (d): total variation restoration with schemes (75) and (24), with λ = 5; (e) and (f ): smoothed total variation based restoration with schemes (75) and (43), with λ = 5 and β = 10.

Figure 4 :

 4 Figure 4: Smoothed total variation based restoration of the noisy images presented on the bottom row of Figure 1 with scheme (43) with β = 10. In both cases, the Lagrange multiplier is µ = 30.

Figure 5 :Figure 6 :

 56 Figure5: Comparisons of the number of iterations for the fixed point loop[START_REF] Ng | On semismooth Newton methods for total variation minimization[END_REF] in algorithm[START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF]: 1 or 200. The L 2 error is given with respect to the number of iterations of (43) (vertical logarithmic scale). Graph (a) is with β = 0.1: after 60 iterations of (43), both errors are the same. Graph (b) is with β = 10: after 10 iterations of (43), both errors are the same. We thus advocate to use only 1 iteration for the fixed point iteration[START_REF] Nikolova | The equivalence of half-quadratic minimization and the gradient linearization iteration[END_REF].

Figure 7 :

 7 Figure7: Speed of convergence (total variation regularization): the L 2 norm of the error is given at each iteration (vertical logarithmic scale). (a) gives the speed of convergence for iterations 1 to 600, and (b) for iterations 100 to 600. On graph (a), from top to bottom are the speed of convergence of the fixed point based algorithm[START_REF] Pazy | On the asymptotic behavior of iterates of nonexpansive mappings in Hilbert space[END_REF] with β = 0.1, the speed of convergence of Chambolle projection algorithm[START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] with τ = 0.249, the speed of Uzawa scheme[START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] with τ = 0.249, the speed of the dual Nesterov algorithm, and the speed of the accelerated dual Nesterov algorithm. On graph (b) are only shown the dual Nesterov algorithm, and the accelerated dual Nesterov algorithm. 1 iteration with the accelerated dual Nesterov algorithm is around 4 times slower than with the dual Nesterov algorithm. But 1 iteration with the dual Nesterov algorithm is itself around 2 times slower than with Uzawa, Chambolle, or scheme (43) (while all three have the same computation time per iteration). To get a highly accurate solution, the accelerated dual Nesterov algorithm seems to be the best choice. However, the dual Nesterov algorithm seems to be the best compromise when one is only interested in getting a good approximation (which is the case for image restoration).

1 λ 2 2 E + y -y m+1 2 E ≤ 1 λ 2 B 2 E 2 E ( 106 ) 2 E

 122212221062 L λ (B * u + λy) -L λ (B * u m + λy m ) * (u -u m ) + λ(y -y m ) = y -y m 2 E + 2 λ B * (u -u m ), y -y m E + 1 λ 2 B * (u -u m )But if we subtract the first line of (9) to the first line of (105), we have:A(u -u m ) = B(y m -y).Taking the inner product with (u -u m ), we deduce:A(u -u m ), u -u m = B(y m -y), u -u m = y m -y, B * (u -u m )(107)Hence:y -y m , B * (u -u m ) = -A(u -u m ), u -u m ≤ -α u -u m 2 E ≤ -α B * 2 B * (u -u m )

  that the convergence rate is O 1 k . This was already an improvement over the O 1

	√	k	bound for non differentiable function [38] (Theorem 3.2.1 page 138). With the dual Nesterov algorithm,
	we have now a O 1 k 2 algorithm for solving problem (11).

Acknowledgements

The author would like to thank Vicent Caselles, Antonin Chambolle and Pierre Weiss for fruitful discussions about this work. The author also would like to thank the anonymous reviewers for their useful comments. This work has been done with the support of the French "Agence Nationale de la Recherche" (ANR), under grant NATIMAGES (ANR-08-EMER-009), "Adaptivity for natural images and texture representations".