Aurélie Fischer
email: aurelie.fischer@parisdescartes.fr

Selecting the length of a principal curve within a Gaussian model

Keywords: Principal curves, model selection, oracle inequality, slope heuristics AMS 2010 Mathematics Subject Classification: 62G08, 62G05

Principal curves are parameterized curves passing "through the middle" of a data cloud. These objects constitute a way of generalization of the notion of first principal component in Principal Component Analysis. Several definitions of principal curve have been proposed, one of which can be expressed as a least-square minimization problem. In the present paper, adopting this definition, we study a Gaussian model selection method for choosing the length of the principal curve, in order to avoid interpolation, and obtain a related oracle-type inequality. The proposed method is practically implemented and illustrated on cartography problems.

Introduction

Principal curves can be thought of as a nonlinear generalization of Principal Component Analysis. Instead of searching for the first principal component of a data cloud, the purpose is to design a curve passing "through the middle" of the observations, as illustrated in Figure 1. Principal curves have many applications in various areas, such as physics (Hastie and Stuetzle [START_REF] Hastie | Principal curves[END_REF], Friedsam and Oren [START_REF] Friedsam | The application of the principal curve analysis technique to smooth beamlines[END_REF]), character and speech recognition (Kégl and Krzyżak [START_REF] Kégl | Piecewise linear skeletonization using principal curves[END_REF], Reinhard and Niranjan [START_REF] Reinhard | Parametric subspace modeling of speech transitions[END_REF]), but also mapping and geology (Brunsdon [10], Stanford and Raftery [START_REF] Stanford | Finding curvilinear features in spatial point patterns: principal curve clustering with noise[END_REF], Banfield and Raftery [START_REF] Banfield | Ice floe identification in satellite images using mathematical morphology and clustering about principal curves[END_REF], Einbeck, Tutz and Evers [START_REF] Einbeck | Exploring multivariate data structures with local principal curves[END_REF][START_REF] Einbeck | Local principal curves[END_REF]), natural sciences (De'ath [START_REF] De'ath | Principal curves: a new technique for indirect and direct gradient analysis[END_REF], Corkeron, Anthony and Martin [START_REF] Corkeron | Ranging and diving behaviour of two 'offshore' bottlenose dolphins, Tursiops sp., off eastern Australia[END_REF], Einbeck, Tutz and Evers [START_REF] Einbeck | Exploring multivariate data structures with local principal curves[END_REF]) and medicine (Wong and Chung [START_REF] Wong | Principal curves to extract vessels in 3D angiograms[END_REF], Caffo, Crainiceanu, Deng and Hendrix [START_REF] Caffo | A case study in pharmacologic colon imaging using principal curves in single photon emission computed tomography[END_REF]). where I = [a, b] is a closed interval of the real line. The original definition of a principal curve, due to Hastie and Stuetzle [START_REF] Hastie | Principal curves[END_REF], relies on the self-consistency property of principal components. A smooth (infinitely differentiable) parameterized curve f (t) = (f 1 (t), . . . , f d (t)) is a principal curve for X if f does not intersect itself, has finite length inside any bounded subset of R d , and is self-consistent, which means that

f (t) = E[X|t f (X) = t]. (1)
Here, the so-called projection index t f (x) is defined by

t f (x) = sup t ∈ I : x -f (t) = inf t ′ x -f (t ′ ) ,
where • denotes the standard Euclidean norm of R d . So, t f (x) is the largest real number t minimizing the Euclidean distance between x and f (t), as shown in Figure 2.

The self-consistency property may be interpreted by saying that each point of the curve f is the mean of the observations projecting on f around this point.

A number of other points of view, more or less related to this original definition, have been proposed thereafter. Tibshirani [START_REF] Tibshirani | Principal curves revisited[END_REF], keeping the self-consistency property, adopts a semiparametric approach and defines principal curves in terms of a mixture model, whereas Delicado [START_REF] Delicado | Another look at principal curves and surfaces[END_REF] generalizes another property of principal components, leading to the notion of "principal curves of oriented points". The definitions of Verbeek, Vlassis and Kröse [START_REF] Verbeek | A soft k-segments algorithm for principal curves[END_REF] and Einbeck, Tutz and Evers [START_REF] Einbeck | Local principal curves[END_REF] are based on local principal components put together and the "locally defined principal curves" of Ozertem and Erdogmus [START_REF] Ozertem | Locally defined principal curves and surfaces[END_REF] correspond to the ridge of a density function. Recently, Genovese, Perone-Pacifico, Verdinelli and Wasserman [START_REF] Genovese | The geometry of nonparametric filament estimation[END_REF] discussed a closely related problem, called nonparametric filament estimation.

x 1 f (t 1 ) x 2 f (t 2 ) x 3 f (t 3 ) x 4 f (t 4 ) x 5 f (t 5 ) x 6 f (t 6 ) x 7 f (t 7 )
Figure 2: The projection index t f . For all i, t i stands for t f (x i ).

In this paper, we will adopt the principal curve definition of Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF], which is closely related to the original one, but presents the advantage of avoiding the implicit formulation. Instead, this definition takes the form of an empirical risk minimization problem, which is easier to handle than (1). In the definition of Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF], a principal curve of length L for X is a parameterized curve minimizing the least-square type criterion

∆(f ) = E inf t∈I X -f (t) 2
over all curves of length not larger than L > 0. Such a principal curve always exists provided E X 2 < ∞, but it may not necessarily be unique. Note that Sandilya and Kulkarni [START_REF] Sandilya | Principal curves with bounded turn[END_REF] have proposed a similar definition, using a constraint on the turn instead of the length of the curve.

In practice, the distribution of X is unknown, and we have at hand a sample X 1 , . . . , X n of independent random variables distributed as X. In this situation, ∆(f ) is replaced by its empirical counterpart

∆ n (f ) = 1 n n i=1 inf t∈I X i -f (t) 2 .
In order to construct a satisfactory principal curve, a good choice of the length is crucial. Indeed, a principal curve constrained to have a too small length will not be able to capture the shape of the data, whereas a too long curve may lead to interpolation problems, as illustrated in Figure 3. In the present contribution, we propose to study the length selection problem, using the approach of non-asymptotic model selection by penalization introduced by Birgé and Massart [START_REF] Birgé | From model selection to adaptive estimation[END_REF] and Barron, Birgé and Massart [START_REF] Barron | Risk bounds for model selection via penalization[END_REF]. To this end, we will consider a Gaussian framework. A related point of view in the context of almost surely bounded random variables is discussed in Biau and Fischer [START_REF] Biau | Parameter selection for principal curves[END_REF]. The rest of the paper is organized as follows. In Section 2, we consider the problem of choosing the length of a principal curve using a Gaussian model selection approach, and show that the curve obtained by minimizing some appropriate penalized criterion satisfies an oracle-type inequality. Section 3 presents some experimental results in the context of cartography. Proofs are collected in Section 4 for the sake of clarity.

[A] [B] [C]

Length selection

We investigate a Gaussian model selection method in order to choose the length of a principal curve. Our context is similar to that of Caillerie and Michel [START_REF] Caillerie | Model selection for simplicial approximation[END_REF], who tackle model selection questions for graphs called simplicial complexes. In the sequel, the Euclidean space R d is equipped with the inner product defined by

u, v = 1 d d j=1 u j v j , (2)
and • denotes the associated Euclidean norm.

We assume that we observe random vectors X 1 , . . . , X n with values in R d following the model

X i = x ⋆ i + σξ i , i = 1, . . . , n, (3)
where the x ⋆ i are unknown, the ξ i are independent standard Gaussian vectors of R d and σ > 0 stands for the noise level, which is supposed known. Let us denote by -→ X = t ( t X 1 , . . . , t X n ) the (column) vector made of all coordinates of the random vectors X i , i = 1, . . . , n. Defining -→ x ⋆ and -→ ξ in the same way, the model (3) can be rewritten under the form

-→ X = -→ x ⋆ + σ -→ ξ .
Let F and G be two fixed points of R d and L a countable subset of ]0, +∞[. We introduce a countable collection {F ℓ } ℓ∈L , where each set F ℓ is a class of parameterized curves f : I → R d with length ℓ and endpoints F and G. Our aim is to select the length ℓ. To do this, we consider the criterion ∆ ′ n defined by

∆ ′ n (f ) = 1 n n i=1 inf t∈I X i -f (t) 2 = 1 n n i=1 inf x i ∈Γ f X i -x i 2 ,
where Γ f denotes the range of the curve f . Due to the definition of the norm • chosen above [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF], this is the empirical criterion ∆ n (f ) normalized by the dimension d. Suppose that, for all ℓ ∈ L, -→ x

(n) ℓ := (x 1ℓ , . . . , xnℓ ) minimizes 1 n n i=1 X i -x i 2 among all -→ x ∈ C ℓ := f ∈F ℓ (Γ f ) n .
In order to determine the length ℓ, our purpose is to minimize in ℓ a criterion of the type

crit(ℓ) = 1 n n i=1 X i -xiℓ 2 + pen(ℓ),
where pen : L → R + is a penalty function, which should avoid the selection of a too large ℓ. Our goal is to design an appropriate penalty. Observe that the classical asymptotic model selection criteria AIC (Akaike [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]), BIC (Schwarz [START_REF] Schwarz | Estimating the dimension of a model[END_REF]) or Mallows'C p (Mallows [START_REF] Mallows | Some comments on C p[END_REF]), which involve the "number of parameters" to be estimated, are not suitable in this framework. Therefore, our approach will rely on the non-asymptotic model selection theory developed by Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF] and Barron, Birgé and Massart [START_REF] Barron | Risk bounds for model selection via penalization[END_REF].

When the considered models are linear subspaces, the penalty can be chosen proportional to the dimension of the model, according to Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF]. Here, the models C ℓ are not linear subspaces of R nd and the dimension must be replaced by another quantity.

In order to measure the complexity of these nonlinear models, we will use metric entropy.

The metric entropy of a set S is given by

H(S, • , ε) = ln N (S, • , ε),
were the covering number N (S, • , ε) is the minimal number of balls with radius ε for the norm • needed to cover S.

Our approach is based on a general model selection theorem for nonlinear Gaussian models (Massart [27]). Let us denote by

• nd the normalized norm of R nd , defined by the inner product -→ u , - → v nd = 1 nd nd i=1 u i v i .
For every ℓ ∈ L, let ϕ ℓ be a function such that ϕ ℓ ≥ φ ℓ , where φ ℓ is given by

φ ℓ (u) = κ u 0 H(C ℓ , • nd , ε)dε, ( 4
)
with κ an absolute constant. We define d ℓ by the equation

ϕ ℓ 2σ √ d ℓ √ nd = σd ℓ √ nd .
Assume that there exists a family of weights {w ℓ } ℓ∈L satisfying ℓ∈L e -w ℓ = Σ < ∞.

Under these assumptions and with this notation, Theorem 4.18 in Massart [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF] can be written in the following manner:

Theorem 2.1. Let η > 1 and pen(ℓ) ≥ η σ 2 nd d ℓ + √ 2w ℓ 2 .
Then, almost surely, there exists a minimizer l of the penalized criterion

crit(ℓ) = 1 n n i=1 X i -xiℓ 2 + pen(ℓ).
Moreover, writing xi := xi l for all i=1,. . . ,n, we have

1 n n i=1 E xi -x ⋆ i 2 ≤ c(η) inf ℓ∈L (d 2 ( -→ x ⋆ , C ℓ ) + pen(ℓ)) + σ 2 nd (Σ + 1)
,

where d 2 ( -→ x ⋆ , C ℓ ) = inf - → y ∈C ℓ 1 n n i=1 y i -x ⋆ i 2 .
This result establishes, for a penalty pen(ℓ) which is large enough, an oracle-type inequality in expectation for the xi , i = 1, . . . , n. Provided a control of the Dudley integral (4) (Dudley [START_REF] Dudley | The sizes of compact subsets of Hilbert space and continuity of Gaussian processes[END_REF]), this theorem will apply in our context and allow us to select the length ℓ of the curve. To assess this integral, we will need some technical lemmas, which are proved in Section 4.

The first step consists in controlling the metric entropy of the classes C ℓ , ℓ ∈ L. Note that, for all ℓ ∈ L, f ∈F ℓ Γ f corresponds to an ellipsoid of R d , as stated in the next lemma. In the sequel, this ellipsoid will be denoted by E ℓ .

Lemma 2.1. Every parameterized curve of R d with endpoints F and G and length ℓ (ℓ > F G), is included in an ellipsoid E ℓ with first principal axis of length ℓ, the other axes having length

λ = √ ℓ 2 -F G 2 .
In particular, in R 2 , E ℓ is an ellipse with foci F and G (see Figure 4), and in R 3 , it is a ellipsoid of revolution around the axis passing through these two points.

• F • G λ = √ ℓ 2 -F G 2 ℓ Figure 4:
In the plane R 2 , ellipse E ℓ with foci F and G and axes ℓ and λ.

We obtain then the following upper bound for

N (C ℓ , • nd , ε), ℓ ∈ L. Lemma 2.2. Suppose that ℓ ≥ λ ≥ ε. The covering number of C ℓ for the normalized norm • nd of R nd satisfies N (C ℓ , • nd , ε) ≤ 2 ε nd (ℓλ d-1 ) n .
Bounding the integral

φ ℓ (u) = κ u 0 H(C ℓ , • nd , ε)dε
for all ℓ ∈ L, we can then define an adequate function ϕ ℓ .

Lemma 2.3. The function ϕ ℓ given by

ϕ ℓ (r) =        κr √ nd   ln 2ℓ 1/d λ 1-1/d r + √ π   if r ≤ λ ϕ ℓ (λ) + (r -λ)ϕ ′ ℓ (λ) if r ≥ λ satisfies, for all r, ϕ ℓ (r) ≥ φ ℓ (r).
Finally, in order to apply Theorem 2.1, we have to assess d ℓ , defined by the equation

ϕ ℓ 2σ √ d ℓ √ nd = σd ℓ √ nd ,
which is the purpose of the next lemma.

Lemma 2.4. Let ϕ ℓ be given by Lemma 2.3. Suppose that

σ ≤ λ 4κ   ln 2 + 1 d ln ℓ λ + √ π   -1
.

Then, equation

ϕ ℓ 2σ √ d ℓ √ nd = σd ℓ √ nd admits a solution d ℓ satisfying d ℓ ≤ 8κ 2 nd ln ℓ 1/d λ 1-1/d 2σκ √ π + π .
We are now in a position to state the main result of this section.

Theorem 2.2. Assume that there exists a family of weights {w ℓ } ℓ∈L such that ℓ∈L e -w ℓ = Σ < ∞, and that, for every ℓ ∈ L,

σ ≤ λ 4κ   ln 2 + 1 d ln ℓ λ + √ π   -1 . ( 5
)
Then, there exist constants c 1 and c 2 such that, for all η > 1, if

pen(ℓ) ≥ ησ 2 c 1 ln ℓ 1/d λ 1-1/d σ + c 2 + 4w ℓ nd , ( 6
)
then, almost surely, there exists a minimizer l of the penalized criterion

crit(ℓ) = 1 n n i=1 X i -xiℓ 2 + pen(ℓ).
Moreover, if xi = xi l for all i = 1, . . . , n, we have

1 n n i=1 E xi -x ⋆ i 2 ≤ c(η) inf ℓ∈L {d 2 ( -→ x ⋆ , C ℓ ) + pen(ℓ)} + σ 2 nd (Σ + 1)
,

where d 2 ( -→ x ⋆ , C ℓ ) = inf - → y ∈C ℓ 1 n n i=1 y i -x ⋆ i 2 .
Let us now comment on the theorem.

The first remark is about the fact that Theorem 2.2 involves unknown constants. Lemma 2.4 indicates that c 1 = 16κ 2 and c 2 = πln(2κ √ π) could be chosen. However, these values are (likely too large) upper bounds. Furthermore, the variance noise σ has been supposed known and is involved in the penalty. Nevertheless, the noise level is generally unknown in practice. In fact, the expression (6) does not provide a penalty function directly, but gives its shape instead. Note that it is possible to estimate σ separately and then proceed by plug-in. However, there is another solution to assess c 1 , c 2 and σ, relying on the slope heuristics. This penalty calibration method introduced by Birgé and Massart [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] (see also Arlot and Massart [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF], Lerasle [START_REF] Lerasle | Optimal model selection in density estimation[END_REF] and Saumard [START_REF] Saumard | The slope heuristics in heteroscedastic regression[END_REF]) precisely allows to tune a penalty known up to a multiplicative constant.

According to the formula binding ℓ and λ, the quantity ln(ℓ 1/d λ 1-1/d ) in the penalty characterizes each model of curves with length ℓ. The other elements varying over the collection of models are the weights {w ℓ } ℓ∈L . For linear models C ℓ with dimension D ℓ , a possible choice for w ℓ is w ℓ = w(D ℓ ) where w(D) = cD + ln |{k ∈ L, D k = D}| and c > 0 (see Massart [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF]). If there is no redundancy in the models dimension, this strategy amounts to choosing w ℓ proportional to D ℓ . By analogy, w ℓ may here be chosen proportional to ln(ℓ 1/d λ 1-1/d ). More formally, we set w ℓ = c ln ℓ 1/d λ 1-1/d , where the constant c > 0 is such that ℓ∈L

1 ℓ c/d λ c(1-1/d) = Σ < +∞.
This choice finally yields a penalty proportional to ln(ℓ 1/d λ 1-1/d ), which may be calibrated in practice thanks to the slope heuristics.

In addition, observe that condition [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] says that the noise level σ should not be too large with respect to λ. In other words, if λ = √ ℓ 2 -F G 2 is of the same order as σ, it is not possible to obtain a suitable principal curve with length ℓ.

Finally, let us point out that due to the exponent n in the covering number in Lemma 2.2-a comment is given in Section 4 after the proof of the lemma (Remark 4.1)-, the penalty shape obtained does not tend to 0 as n tends to infinity. This point is intrinsically related to the geometry of the problem. Indeed, its resolution is not made easier by increasing the size of the sample, since nothing has been specified about the repartition of the x ⋆ i 's. A possible direction for future research could consist in dealing with the framework in which these points are distributed along the curve following a uniform distribution.

Experimental results

In this section, we propose to illustrate the length selection method practically. The experiments presented here are carried out with the software MATLAB. As announced in Section 2, the penalty given in Theorem 2.2 will be calibrated thanks to the slope heuristics. Two strategies may be used : the dimension jump method consists in identifying an abrupt jump in the models complexity, whereas the other solution is to observe that the empirical contrast is proportional to the penalty shape for complex models and use the slope of this line to assess the constant.

In this practical implementation, we considered polygonal lines, which present the advantage that projecting on the curve reduces to projection on a segment. However, the method described below could probably be replaced by a more sophisticated technique dealing with smooth curves. In the sequel, the maximal number k of segments is taken large enough, to ensure that the only parameter reflecting the complexity of the curve is the length. Then, the length ℓ of the principal curve is chosen as follows:

1. For a range of values of the length ℓ, compute fℓ by minimizing the empirical criterion ∆ n (f ) and record

∆ n ( fℓ ) = 1 n n i=1 ∆( fℓ , X i ).
2. Let w ℓ be proportional to ln ℓ 1/d λ 1-1/d and consider a penalty of the form

pen(k, ℓ) = c ln(ℓ 1/d λ 1-1/d ).
3. Select the constant ĉ using the slope heuristics.

4. Retain the curve fl obtained by minimizing the penalized criterion

crit(ℓ) = 1 n n i=1 X i -xiℓ 2 -2ĉ ln(ℓ 1/d λ 1-1/d ).
In step 1 of the algorithm, the criterion ∆ n (f ) is minimized by using a MATLAB optimization routine.

The weights w ℓ were chosen as suggested in the discussion after Theorem 2.2. This is a convenient choice, which does not modify the penalty shape.

To apply the slope heuristics in step 3, we employ the MATLAB package CAPUSHE, implemented by Baudry, Maugis and Michel [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF]. We tried both the dimension jump and the slope estimation methods, which results were found to be very similar.

Finally, recall that the endpoints F and G of the principal curve have been assumed to be fixed. From a practical point of view, several methods can be employed to choose these two points from the observations. A possible solution is to define F and G with the aid of the points that are farthest from each other in the minimum spanning tree of the data (or of some subset of the data). Figure 5 gives some examples of such trees, which can be constructed using Kruskal's algorithm [START_REF]On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF] or Prim's one [START_REF] Prim | Shortest connection networks and some generalizations[END_REF].

Figure 5: Some examples of minimum spanning trees.

We will now present two examples of applications to mapping. Indeed, Brunsdon [START_REF] Brunsdon | Path estimation from GPS tracks[END_REF] has shown that principal curves may be useful in that area, in order to estimate paths from GPS tracks. More specifically, principal curves can serve as a means to compute an average path from GPS data registered by several people moving on a given street.

The results obtained for some hairpin street located in France and the Labyrinth of the Jardin des Plantes in Paris are visible in Figure 6 and 7 respectively. Each figure gives first an air photography of the place and the corresponding GPS tracks data points. Then, the resulting principal curve is shown both on the data cloud and as an overlay on the photography, which allows to assess the performance of the method. Moreover, the principal curves fitted using our model selection approach (denoted by LS in the sequel) can be compared for both data sets to those obtained with a benchmark algorithm. Indeed, Figure 8 gives the outputs of the Polygonal Line Algorithm, which is based on a local control of the curvature and was proposed by Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF] (PL hereafter). On the whole, the considered streets are correctly recovered. For the hairpin road, we observe that the PL output is smoother, but the right-hand end of the curve looks better on the LS result. Note that the direction taken by this part of the PL curve is wrong.

Regarding the Labyrinth, LS perform quite well and yields this time the smoothest curve. The PL principal curve is not at all as smooth as before: this can be explained by the fact that we used for both experiments the default parameters of this algorithm, which may be more or less suitable depending on the different characteristics of the data set, such as curvature or sampling density.

Noting that for the first example, there is a part of the street somewhat hidden by trees, a particularly interesting application could be to use principal curves to draw a map of paths in a forest.

Proofs

Proof of Lemma 2.1

Let c = F G/2. Note that ℓ > 2c. In a well-chosen orthonormal coordinate system of R d , F has coordinates (-c, 0, . . . , 0) and G (c, 0, . . . , 0). A curve with length ℓ and endpoints F and G is included in the set delimited by the points M(x 1 , . . . , x d ) such that

MF + MG = ℓ.
Let us show that this equation defines an ellipsoid with first principal axis ℓ, the other axes having length λ. Let M(x 1 , . . . , x d ) be such that MF + MG = ℓ. Then,

MF 2 = (x 1 + c) 2 + d j=2 x 2 j and MG 2 = (x 1 -c) 2 + d j=2 x 2 j .
Therefore,

MF -MG = MF 2 -MG 2 MF + MG = (x 1 + c) 2 -(x 1 -c) 2 ℓ = 4x 1 c ℓ .
As a result, on the one hand,

(MF + MG) 2 + (MF -MG) 2 = ℓ 2 + 16x 2 1 c 2 ℓ 2 ,
and on the other hand,

(MF + MG) 2 + (MF -MG) 2 = 2(MF 2 + MG 2 ) = 4 d j=1 x 2 j + 4c 2 .
Hence,

ℓ 2 + 16x 2 1 c 2 ℓ 2 = 4 d j=1 x 2 j + 4c 2 ,
which may be rewritten

x 2 1 1 - 4c 2 ℓ 2 + d j=2 x 2 j = ℓ 2 4 -c 2 , ( 7
)
or, equivalently,

x 2 1 ℓ 2 /4 + d j=2 x 2 j ℓ 2 /4 -c 2 = 1, (8)
where ℓ 2 /4c 2 > 0 since ℓ > 2c. In other words, the point M belongs to an ellipsoid with one axis of length ℓ and d -1 axes of length 2

ℓ 2 /4 -c 2 = √ ℓ 2 -F G 2 = λ.
Reciprocally, if M(x 1 , . . . , x d ) satisfies equation [START_REF] Birgé | From model selection to adaptive estimation[END_REF], with ℓ 2 4 c 2 > 0, then

MF 2 = (x 1 + c) 2 + d j=2 x 2 j = (x 1 + c) 2 + ℓ 2 4 -c 2 -x 2 1 + 4x 2 1 c 2 ℓ 2 = 2x 1 c + ℓ 2 4 + 4x 2 1 c 2 ℓ 2 = ℓ 2 + 2x 1 c ℓ 2 . Hence, MF = ℓ 2 + 2x 1 c ℓ . Similarly, MG = ℓ 2 - 2x 1 c ℓ . Now, |x 1 | ≤ ℓ 2 4c
, since otherwise

x 2 1 ℓ 2 /4 > ℓ 2 4c 2 > 1, contradicting equation (8). Finally, MF +MG = ℓ 2 + 2x 1 c ℓ + ℓ 2 - 2x 1 c ℓ = ℓ.

Proof of Lemma 2.2

First, we shall compute the covering number of an ellipsoid E ℓ of R d . The next lemma is a particular case of Proposition 5 in von Luxburg, Bousquet and Schölkopf [START_REF] Luxburg | A compression approach to support vector model selection[END_REF].

Lemma 4.1. Assume that ℓ ≥ λ ≥ ε. The number of balls of radius ε needed to cover E ℓ , ellipsoid in dimension d with principal axes ℓ, λ, . . . , λ, satisfies

N (E ℓ , • , ε) ≤ 2 ε d ℓλ d-1 .
Proof of Lemma 4.1. The number of balls of radius ε needed to cover E ℓ satisfies

N (E ℓ , • , ε) ≤ ℓ ε + 1 λ ε + 1 d-1
, where ⌊y⌋ denotes the floor of y, i.e., the largest integer less than or equal to y. Indeed, the ellipsoid E ℓ is inscribed in a parallelepiped with sides of lengths ℓ, λ, . . . , λ, and the number of balls of radius ε needed to cover a parallelepiped with sides of length c 1 , . . . , c d is

d j=1 c j ε + 1 .
By assumption, ℓ ≥ λ ≥ ε, so that ℓ ε + 1 ≤ 2ℓ ε and λ ε + 1 ≤ 2λ ε . Hence,

N (E ℓ , • , ε) ≤ 2 ε d ℓλ d-1 .
Thus, according to Lemma 4.1,

N (E ℓ , • , ε) ≤ 2 ε d ℓλ d-1 .
Let U be a collection of at most 2

ε d ℓλ d-1 centers of balls in (R d , • ), corresponding to an ε-covering of E ℓ . For each vector -→ u = t ( t u 1 , . . . , t u n ) ∈ R nd
, where the u i 's are elements of U, we have n i=1 B(u i , ε) ⊂ B( -→ u , ε) (balls for the normalized norm of R d and R nd respectively). Consequently,

N (C ℓ , • nd , ε) ≤ 2 ε nd (ℓλ d-1 ) n .
Remark 4.1. It is not possible to get rid of the exponent n in this covering number computation. Indeed, if we consider only one curve f of length ℓ, covered by N balls of R d , N n balls of R nd are needed to cover (Γ f ) n , since every point among the n considered points along the curve can be in any one of the N balls. However, our upper bound for the number of balls in R nd needed to recover C ℓ is probably too large. Indeed, since the n points are constrained to be located on the same curve of length ℓ, we would not need to use all balls of R nd obtained by combining the centers of the balls of R d covering E ℓ . We could undoubtedly get a better upper bound by solving the combinatorial problem consisting in counting all acceptable combinations of balls in R d for a class C ℓ .

Proof of Lemma 2.3

We begin by a technical lemma which will be useful for bounding integrals.

Lemma 4.2. For x ∈ ]0, 1],

x 0 ln 1 t dt ≤ x ln 1 x + √ π .
Proof of Lemma 4.2. We have

x 0 ln 1 t dt = t ln 1/t x 0 + x 0 1 2 ln 1/t dt = x ln 1/x + 1 √ 2 +∞ √ 2 ln 1/x e -u 2 /2 du ≤ x( ln 1/x + √ π).
This inequality is due to the fact that, for a ≥ 0,

1 √ 2π +∞ a e -u 2 /2 du ≤ e -a 2 /2 . ( 9
)
Indeed, if g denotes the function defined by g(a) = e -a 2 /2 -1 √ 2π

+∞ a e -u 2 /2 du, then g ′ (a) = e -a 2 /2 ( 1 √ 2πa). This function g is increasing on [0, 1/ √ 2π] and decreasing on [1/ √ 2π, +∞]. Since g(0) = 1/2 and lim +∞ g = 0, we obtain g(a) ≥ 0 for all a ≥ 0. Hence, inequality (9) is proved.

Back to the proof of Lemma 2.3, note that according to Lemma 2.2, the metric entropy

H(C ℓ , • nd , ε) satisfies H(C ℓ , • nd , ε) ≤ nd ln 2ℓ 1/d λ 1-1/d ε . If r ≤ λ, φ ℓ (r) = κ r 0 H(C ℓ , • nd , ε)dε ≤ κ √ nd r 0 ln 2ℓ 1/d λ 1-1/d ε dε. The change of variables t = ε 2ℓ 1/d λ 1-1/d gives φ ℓ (r) ≤ 2κ √ ndℓ 1/d λ 1-1/d r 2ℓ 1/d λ 1-1/d 0 ln 1 t dt. Now, Lemma 4.2 indicates that for x ∈ ]0, 1], x 0 ln 1 t dt ≤ x ln 1 x + √ π .
Thus,

φ ℓ (r) ≤ κr √ nd   ln 2ℓ 1/d λ 1-1/d r + √ π   .
For r ≤ λ, let ln 2 ℓ λ

1/d ≥ √ π - 1 2   ln 2 ℓ λ 1/d   -1/2 ≥ √ π - 1 2 √ ln 2 ≥ 0,
which shows that ϕ ′ ℓ (λ) ≥ H(λ). Then, let ϕ ℓ (r) = ϕ ℓ (λ) + (rλ)ϕ ′ ℓ (λ) for r ≥ λ, so that, finally, φ ℓ (r) ≤ ϕ ℓ (r) for all r.

Proof of Lemma 2.4

Notice first that ϕ ℓ is concave. Indeed, the second derivative of the restriction ϕ ℓ | ]0,λ] of ϕ ℓ to ]0, λ] is equal to observe that looking for a solution d ℓ amounts to solving, for r > 0,

- κ √ nd 2r
ϕ ℓ (r) = √ nd 4σ r 2 .
This equation admits a unique solution r ℓ = 2σ d ℓ nd , since ϕ ℓ is concave and r → r 2 convex. Moreover, the solution r ℓ satisfies r ℓ ≤ λ if, and only if,

ϕ ℓ (λ) ≤ √ nd 4σ λ 2 , that is κλ √ nd   ln 2 + 1 d ln ℓ λ + √ π   ≤ √ nd 4σ λ 2 ,
which means that

σ ≤ λ 4κ   ln 2 + 1 d ln ℓ λ + √ π   -1
.

If this condition is satisfied, the equation becomes

κr ℓ √ nd   ln 2ℓ 1/d λ 1-1/d r ℓ + √ π   = √ nd 4σ r 2 ℓ ,
which is equivalent to 4σκ

  ln 2ℓ 1/d λ 1-1/d r ℓ + √ π   = r ℓ .
So, 4σκ √ π ≤ r ℓ , and then,

r ℓ ≤ 4σκ   ln ℓ 1/d λ 1-1/d 2σκ √ π + √ π   .
Since r ℓ = 2σ d ℓ nd , we obtain

d ℓ ≤ 8κ 2 nd ln ℓ 1/d λ 1-1/d 2σκ √ π + π .

Figure 1 :

 1 Figure 1: An example of principal curve.

Figure 3 :

 3 Figure 3: Principal curves fitted with [A] a too small length, [B] a too large length and [C]an appropriate one.

Figure 6 :

 6 Figure 6: Principal curve fitted with the LS algorithm for the hairspin data.

Figure 7 :

 7 Figure 7: Principal curve fitted with the LS algorithm for the Labyrinth data.

Figure 8 :

 8 Figure 8: Principal curves fitted with the PL Algorithm. [A] Hairpin street data. [B]Labyrinth data.

κ λ 0 H

 0 ℓ , • nd , ε)dε = (C ℓ , • nd , ε)dε + κ r λ H(C ℓ , • nd , ε)dε ≤ φ ℓ (λ) + (rλ)H(λ) ≤ ϕ ℓ (λ) + (rλ)ϕ ′ ℓ (λ),whereH(λ) = κ H(C ℓ , • nd , λ). Indeed, ϕ ℓ (λ) ≥ φ ℓ (λ)by definition of ϕ ℓ , and on the other hand ϕ ′ ℓ (λ) -H(λ

  ϕ ℓ | ]0,λ] is concave. As ϕ ℓ is obtained by extending ϕ ℓ | ]0,λ] using the tangent to this function in λ, ϕ ℓ is also concave.

Figure 9 :

 9 Figure 9: Function ϕ ℓ κ √ nd for d = 2, ℓ = 6 and λ = 3.


Multiplicative Noise Removal Using L1 Fidelity on Frame Coefficients

Keywords: Curvelets, Douglas-Rachford splitting, ℓ 1 data-fidelity, Multiplicative noise removal, Nonsmooth optimization, Proximal calculus, Tight frames, TV regularization, Variational methods

We address the denoising of images contaminated with multiplicative noise, e.g. speckle noise. Classical ways to solve such problems are filtering, statistical (Bayesian) methods, variational methods, and methods that convert the multiplicative noise into additive noise (using a logarithmic function), apply a variational method on the log data or shrink their coefficients in a frame (e.g. a wavelet basis), and transform back the result using an exponential function.

We propose a method composed of several stages: we use the log-image data and apply a reasonable under-optimal hard-thresholding on its curvelet transform; then we apply a variational method where we minimize a specialized hybrid criterion composed of an ℓ 1 data-fidelity to the thresholded coefficients and a Total Variation regularization (TV) term in the log-image domain; the restored image is an exponential of the obtained minimizer, weighted in a such way that the mean of the original image is preserved. Our restored images combine the advantages of shrinkage and variational methods and avoid their main drawbacks. Theoretical results on our hybrid criterion are presented. For the minimization stage, we propose a properly adapted fast scheme based on Douglas-Rachford splitting. The existence of a minimizer of our specialized criterion being proven, we demonstrate the convergence of the minimization scheme. The obtained numerical results clearly outperform the main alternative methods especially for images containing tricky geometrical structures.

Introduction

In various active imaging systems, such as synthetic aperture radar (SAR), laser or ultrasound imaging, the data representing the underlying (unknown image) S o : Ω → R + , Ω ⊂ R 2 , are corrupted with multiplicative noise. Such a noise severely degrades the image (see Fig. 2(a)-(b) as well as the noisy images in section 6). In order to increase the chance to restore a better image, several independent measurements for the same scene should be realized, thus yielding a set of data:

Σ k = S o η k + n k , ∀k ∈ {1, • • • , K}, K ≥ 1, (1)
where η k : Ω → R + , and n k represent the multiplicative and the additive noise relevant to each measurement k. The additive noise n k has usually a very weak effect and is systematically neglected in

0 5 -6 0 2 -1 0 1 1 2 -1 0 1 -1 0 1 (a) pdf(η = η k ) (b) pdf(n) (c) Hist(W n) (d) pdf(η) (e) pdf(n) (f) Hist(W n)
η k : pdf(η k ) = µ e -µη k 1l R+ (η k ) where 1l R+ (t) def = 1 if t ≥ 0, 0 otherwise. ( 2
)
The data used for denoising is the average of the set of all K measurements that could be realized (see e.g. Fig. 2(b)). In classical SAR modeling one takes1 µ = 1 in (2), so the usual data production model reads (see e.g. [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Tur | When is speckle noise multiplicative?[END_REF][START_REF] Ulaby | Handbook of Radar Scattering Statistics for Terrain[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] among many other references):

Σ = 1 K K k=1 Σ k = S o 1 K K k=1 η k = S o η . (3)
Since all η k are independent, the resultant mean of the multiplicative noise η in (3) follows a Gamma distribution, see Fig. 1(d) for K = 10:

η = 1 K K k=1 η k : pdf(η) = K K η K-1 Γ(K) e -Kη , (4)
where Γ is the usual Gamma-function and pdf stands for probability density function.

Conversion of the multiplicative noise into additive noise. A large variety of methods-references are given in § 1.1-rely on the conversion of the multiplicative noise into additive noise using

v = log Σ = log S o + log η = u o + n. (5)
In this case the probability density function of n reads (see Fig. 1(b)-(e)):

n = log η : pdf(n) = K K e K(n-e n ) Γ(K) . (6)
In our experiments, we will consider noisy images for K = 1 (see e.g. Fig. 2(a)), in which case pdf(η) and pdf(n) are trivial to derive from ( 4) and [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF].

Preexisting restoration methods

Various adaptive filters for multiplicative noise removal have been proposed, see e.g. [START_REF] Krissian | Oriented speckle reducing anisotropic diffusion[END_REF][START_REF] Yu | Speckle reducing anisotropic diffusion[END_REF] and references therein. Experiments have shown that filtering methods work well when the multiplicative noise is weak, i.e. when K is large. However, in practice it is seldom possible to get a large K.

Bayesian MAP, regulrization and variational methods

Even though the methods belonging to this class are inspired by different motivations, it is well known that given v, the general approach define the restored function û as

û = argmin u F v (u),
where u is defined on a continuous or a finite domain Ω ⊂ R 2 , F v is a criterion combining data-fidelity Ψ with priors (regularization) Φ balanced via a parameter ρ > 0:

F v (u) = ρΨ(u) + Φ(u) for      Ψ(u) = Ω ψ u(ξ), v(ξ) dξ, Φ(u) = Ω φ(|∇u(ξ)|) dξ, ξ = (ξ 1 , ξ 2 ) ∈ Ω. (7)
In what follows, the overscript " . " denotes a minimizer of a given criterion and when there is no ambiguity, the corresponding restored image. In [START_REF] Antoniadis | Wavelet thresholding for some classes of non-gaussian noise[END_REF], ψ : R + → R + assesses closeness to data, ∇ stands for gradient (possibly in a distributional sense or a discrete approximation when Ω is finite), φ : R + → R + is an increasing function modeling the prior on |∇u| and | • | is a norm on the gradient field. Rudin, Osher

and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] proved that when Ω is finite, φ(t) = t yields solutions û that preserve edges. The resultant (prior) regularization term, known as Total Variation (TV), is very popular. It reads

u TV def = Ω |∇u(ξ)| dξ,
where |∇u| = (∂u/∂ξ 1 ) 2 + (∂u/∂ξ 2 ) 2 .

Next we sketch the main modern MAP or variational methods that were used for multiplicative noise removal. Such methods were applied to raw-data ( 3)-( 4) as well as to log-data ( 5)- [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF]. In [START_REF] Achim | SAR image filtering based on the heavy-tailed rayleigh model[END_REF], the authors develop a Bayesian MAP estimator on the log-image where the data-fitting term Ψ is derived from [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF] and a specially designed prior term. The method of Aubert-Aujol (AA) [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF] uses the raw data (3)-( 4) and Ψ is the log-likelihood derived from [START_REF] Achim | SAR image filtering based on the heavy-tailed rayleigh model[END_REF]. The criterion reads (AA) F Σ (S) = ρ log S(ξ) + Σ(ξ) S(ξ) dξ + S TV .

The method proposed in [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF] combines the log-likelihood associated to [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF] and a smoothed TV regularization. A simple and fast method using graph-cut minimization is proposed in [START_REF] Denis | SAR Image Regularization With Fast Approximate Discrete Minimization[END_REF].

The methods in [START_REF] Rudin | Total variation based image restoration with free local constraints[END_REF][START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF] consider white Gaussian multiplicative noise with mean 1 and a very small variance. In [START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF], Rudin Lions and Osher (RLO) minimize a criterion of the form

F Σ (S) = ρ Σ S -1 2 2 + S TV . (9)
The methods proposed by Shi and Osher (SO) in [START_REF] Shi | A nonlinear inverse scale space method for a convex multiplicative noise model[END_REF] deal with various noises and provide iterative TV regularization using relaxed inverse scale space (RISS) flows. The first one uses the log-data v and is based on the criterion

F v (u) = ρ u -v 2 2 + u TV .
The corresponding RISS flow reads (SO)

u t = div ∇u |∇u| + ρ(v -u + v), v t = α(v -u), α > 0, with v 0 = 0 and u 0 = mean(v). ( 10
)
The second one generalizes several multiplicative noise models in the image domain [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF][START_REF] Darbon | The use of levelable regularization functions for MRF restoration of SAR images while preserving reflectivity[END_REF][START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF] and the RISS flow uses iterative TV regularization on exp (S) for convergence reasons. According to the authors, both methods provide similar solution qualities. We focus on the first one mainly because of its simplicity.

Whatever smooth data-fidelity is chosen, it was proven [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF][START_REF] Nikolova | Weakly constrained minimization. Application to the estimation of images and signals involving constant regions[END_REF] that TV regularization yields images containing numerous constant regions (the well known stair-casing effect), hence textures and fine details are removed; this is clearly visible in the restorations obtained using these methods as reported in section 6. The results of [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] initiated a flood of papers to construct smooth edge-preserving convex functions φ, see e.g. [START_REF] Acar | Analysis of bounded variation penalty methods for ill-posed problems[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Vogel | Iterative method for total variation denoising[END_REF], and [START_REF] Aubert | Mathematical problems in image processing[END_REF] for a recent overview. Even though smoothness at zero alleviates stair-casing, it was proven in [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF] that these functions φ lead to images with underestimated edges. This is particularly annoying if the sought-after image has neat edges or spiky features since the latter are eroded. Nevertheless, the good point is that these methods enable the introduction of priors in the image domain.

Multiscale shrinkage for the log-data

Many authors-e.g. [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Achim | SAR image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Fukuda | Suppression of speckle in synthetic aperture radar images using wavelet[END_REF][START_REF] Pizurica | A review of wavelet denoising in MRI and ultrasound brain imaging[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] and references therein-focus on restoring the log-data v in (5) using decompositions into some multiscale frame for L 2 (R 2 ). Let us remind that the analysis operator W of a frame W = w i : i ∈ I , where I is a set of indices, i.e. (W z)

[i] = z, w i , ∀i ∈ I, satisfies the generalized Parseval condition c 1 z 2 2 ≤ i∈I | z, w i | 2 ≤ c 2 z 2 2 , ∀z ∈ L 2 (R 2 )
, where c 2 ≥ c 1 > 0 are the frame bounds. When c 1 = c 2 , the frame is said to be tight. In particular W can be a wavelet basis. The frame is chosen so that u o admits a sparse representation, i.e.

u o ≈ i∈J ζ i w i for some ζ i ∈ R with #J ≪ #I.
The usual strategy is to decompose the log-data v into an well-adapted W :

y = W v = W u o + W n . ( 11
)
The rationale is that the noise W n in y approaches a Gaussian distribution-see Figs. 1(c) and (f)according to the Central Limit Theorem. The tendency to normality gets better as K increases. Under different frameworks, coefficients are denoised using shrinkage estimators:

y T [i] = T (W v)[i] , ∀i ∈ I, (12)
where T : R → R is a symmetric function satisfying 0 ≤ T (t) ≤ t, ∀t ≥ 0. The most widely used examples are soft and hard thresholding inaugurated in [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], and are given in ( 16) and (36) later on. Since then, various shrinkage functions T have been explored, e.g. [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF][START_REF] Belge | Wavelet domain image restoration with adaptive edge-preserving regularization[END_REF][START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Simoncelli | Noise removal via Bayesian wavelet coding[END_REF][START_REF] Wang | Solution of inverse problems in image processing by wavelet expansion[END_REF] to name only a few. See also those derived in a Bayesian framework and applied to multiplicative noise removal e.g. in [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Achim | SAR image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF].

Let W = { w i : i ∈ I} be a dual frame. The associated dual operator W : y → i∈I y[i] w i is a left inverse of W , that is v = i∈I (W v)[i] w i .
Then a denoised log-image v T is generated by expanding the shrunk coefficients y T in the dual frame reads

v T = i∈I T ((W v)[i]) w i = i∈I T (y[i]) w i , (13)
and the sought-after image reads S T = exp(v T ).

The major problems with shrinkage denoising methods, as sketched in ( 12)- [START_REF] Bobichon | Regularized multiresolution methods for astronomical image enhancement[END_REF], is that shrinking large coefficients entails an erosion of the spiky image features, while shrinking small coefficients towards zero yields Gibbs-like oscillations in the vicinity of edges and a loss of texture information. On the other hand, if shrinkage is not strong enough, some coefficients bearing mainly noise will remain almost unchanged-we call such coefficients outliers-and (13) suggests they generate artifacts with the shape of the functions w i . An illustration can be seen in Fig. 2(b-h). Several improvements, such as translation invariant thresholding [START_REF] Coifman | Translation-invariant de-noising[END_REF] and block-Stein thresholding (BS) [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF], were brought in order to alleviate these artifacts. Results obtained using the BS method are presented in section 6: the above mentioned artifacts remain visible. Another inherent difficulty comes from the fact that coefficients between different scales are not independent, as usually assumed, see e.g. [START_REF] Antoniadis | Wavelet thresholding for some classes of non-gaussian noise[END_REF][START_REF] Belge | Wavelet domain image restoration with adaptive edge-preserving regularization[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Simoncelli | Bayesian denoising of visual images in the wavelet domain[END_REF]. In summary, as shrinkage-based estimation relies on sparsity of the representation, it is able to capture efficiently faint structures in the image. But this comes at the price of an intricate choice of the shrinkage function and the associated parameters (e.g. threshold).

Our approach is hybrid

We initially restore the log-data [START_REF] Achim | SAR image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF] and then derive the restored image Ŝ. Our objective is to avoid the main drawbacks of variational and sparsity-based shrinkage methods and to take benefit of the best of both worlds. A way to achieve such a goal is to combine the information brought by the coefficients of the frame-transformed data along with pertinent regularization in the domain of the log-image. This idea for the purpose of additive Gaussian noise removal have been investigated in several papers [START_REF] Bobichon | Regularized multiresolution methods for astronomical image enhancement[END_REF][START_REF] Candès | New multiscale transforms, minimum total variation synthesis. Applications to edge-preserving image reconstruction[END_REF][START_REF] Chan | Total variation improved wavelet thresholding in image compression[END_REF][START_REF] Coifman | Combining the calculus of variations and wavelets for image enhancement[END_REF][START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF][START_REF] Froment | Artifact free signal denoising with wavelets[END_REF][START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF][START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF].

Although guided by different arguments, hybrid methods amount to define the restored function û as

û ∈ argmin u Φ(u) subject to |(W (u -v)) [i]| ≤ µ i , ∀i ∈ I .
If the use of an edge-preserving regularization, such as TV for Φ is a pertinent choice, the strategy for the selection of parameters (µ i ) i∈I is more tricky since it must take into account the magnitude of the data coefficients (y[i]) i∈I . However, deciding on the value of µ i based solely on y[i], as done in these papers, is too rigid since there are either correct data coefficients that incur smoothing (µ i > 0), or noisy coefficients that are left unchanged (µ i = 0). A way to deal with this situation is to determine (µ i ) i∈I based both on the data and on the prior term Φ. According to the theoretical results derived in [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF][START_REF] Nikolova | A variational approach and impulse noise[END_REF],

this objective can be carried out by defining a non-smooth data-fitting term for the coefficients, as done by some of the authors of this paper in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF] which gave rise to very successful numerical results.

To the best of our knowledge, hybrid methods have never been applied to multiplicative noise removal before, whereas the latter is a challenging problem: it arises in important applications but up to now, there is no entirely satisfactory methods to solve it.

We propose a method where the restored log-image û is defined as the minimizer of a criterion composed of an ℓ 1 -fitting to the (suboptimally) hard-thresholded frame coefficients of the log-data and a TV regularization in the log-image domain (section 2). The rationale behind this choice and several theoretical properties are presented as well. This method uses some ideas from a previous work of some of the authors [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. The minimization scheme to compute the log-restored image uses a Douglas-Rachford splitting scheme specially adapted to our criterion (section 3). It involves original derivations and proofs.

The sought-after image is of the form B exp (û) where B is a bias-correction term derived from the noise distribution (section 4). The full algorithm to remove multiplicative noise is summarized in section 5. A large variety of experiments and comparisons to other contemporary multiplicative noise removal methods are presented in section 6. Concluding remarks and open questions are discussed in section 7.

Restoration of the frame coefficients of the log-data

In this section we consider how to restore the frame coefficients y = W v of the log-data image v obtained according to [START_REF] Achim | SAR image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF]. We focus on methods which, for a given preprocessed data set, lead to convex optimization problems. The denoised coefficients are denoted by x.

We assume that v ∈ L 2 (Ω) which ensures that y = W v ∈ ℓ 2 (I).

Specific requirements to restore the coefficients

Given the log-data v obtained according to [START_REF] Achim | SAR image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF], we first apply a frame transform as in [START_REF] Belge | Wavelet domain image restoration with adaptive edge-preserving regularization[END_REF] to get

y[i] = w i , v = w i , u o + w i , n , ∀i ∈ I, (14)
where u o denotes the unknown original log-image. The noise contained in the i-th datum reads n, w i ; its distribution is of the form displayed in Fig. 1(c) or (f). However the signal to noise ratio (SNR) of the coefficients is ill-assorted. When u o has a sparse representation in the frame, many coefficients contain only noise. For this reason, we apply a hard-thresholding to all coefficients

y TH [i] def = T H y[i] , ∀i ∈ I, (15)
where the hard-thresholding operator T H , with a threshold T , reads [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] T

H (t) = 0 if |t| ≤ T, t otherwise. ( 16
)
The resultant set of coefficients is systematically denoted by y TH . We choose an underoptimal threshold T in order to preserve as much as possible the information relevant to edges and to textures, an important part of which is contained in small coefficients. Let's point out that with hard thresholding, all kept coefficients are unaltered and thus contain the original information on the sought-after image.

Let us consider

v TH = i∈I y TH [i] w i = i∈I1 y[i] w i , (17)
where

I 1 = {i ∈ I : |y[i]| > T }. (18)
The image v TH contains artifacts with the shape of the w i 's, for all y[i], that are dominated by noise and above the threshold T , as well as a lot of information about the fine details in the original (unknown) log-image u o . In all cases, whatever the choice of T , the image v TH is unsatisfactory-see Fig. 2

(c)-(h).

Given the under-thresholded data y TH , we focus on hybrid methods where the restored coefficients x minimize a function where Ψ is a data-fidelity term in the domain of the frame coefficients and Φ is an edge-preserving regularization term in the log-image domain. The restored log-image û is of the form

F y : ℓ 2 (I) → R ∪ {∞} of the form:      x = argmin x∈ℓ 2 (I) F y (x) F y (x) = Ψ(x, y TH ) + Φ( W x), (19)
û = W x . (20)
Let us denote

I 0 = I \ I 1 = {i ∈ I : |y[i]| ≤ T }, (21)
where I 1 is given in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF]. In order to specify the shape of Ψ and Φ, we analyze the information borne by the coefficients y TH [i] relevant to I 0 and to I 1 .

(I 0 ) The coefficients y[i] for i ∈ I 0 usually correspond to high-frequency components which can be of the two types described below.

[i] = y[i], by (17) we get v TH = j =i x[j] w j + y[i] w i which
shows that û would contain an artifact with the shape of w i . For this reason, x[i] must be restored according to the prior Φ.

A specialized hybrid criterion

This analysis clearly defines the goals that the minimizer x of F y is expected to achieve. In particular, x must involve an implicit classification between coefficients that fit to y TH exactly and coefficients that are restored according to the prior term Φ. In short, restored coefficients have to fit y TH exactly if they are in accordance with the regularization term Φ and have to be restored via the latter otherwise. Since [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF][START_REF] Nikolova | A variational approach and impulse noise[END_REF] it is known that criteria F y where Ψ is non-smooth at the origin (e.g. ℓ 1 ) can satisfy x[i] = y TH [i] for coefficients that are in accordance with the prior Φ, while the other coefficients are restored according to Φ, see also [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. For these reasons, we focus on a criterion of the form [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] where

Ψ(x) = i∈I1 λ i |(x -y)[i]| + i∈I0 λ i |x[i]| = Λ(x -y TH ) 1 , (22) for
Λ def = diag{λ i : i ∈ I} , (23)
Φ(x) = Ω ∇ W x dξ = W x TV . ( 24
)
Remark 1 We should emphasize that the TV regularization term Φ in our criterion will not favor minimizers û = W x that involve constant regions unlike the usual variational methods, as discussed in § 1.1.1. The reason is that the non differentiability of Ψ at y TH ensures that x[i] = y TH [i] for a certain number of coefficients [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF]. Then W x keeps some fixed structures of the form y TH [i] w i which can prevent from stair-casing in û, provided that there is no large sets of outliers corresponding to the same region of the image all of which are restored according to the TV term Φ. Such an important benefit clearly depends on the level of the noise, the threshold T and the choice of (λ i ) i∈I .

Well-posedness of the minimization problem

The theorem below ensures the existence of a minimizer of the criterion F y defined above. Its proof can be found in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF].

Theorem 1 [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF] For y ∈ ℓ 2 (I) and T > 0 given, consider F y as defined by (15), ( 19) and ( 22)- [START_REF] Coifman | Combining the calculus of variations and wavelets for image enhancement[END_REF], where Ω ∈ R 2 is open, bounded and its boundary ∂Ω is Lipschitz, and

W = w i : i ∈ I is a frame of L 2 (Ω). Suppose that 1. W is the pseudo inverse of W , i.e. W = (W * W ) -1 W * where W * is the adjoint operator; 2. λ min = min i∈I λ i > 0.
Then F y has a minimizer x ∈ ℓ 2 (I).

Remark 2 Since W is a multiscale frame for L 2 (R 2 ) and x ∈ ℓ 2 (I), the restored log-image û as given in (20) satisfies û ∈ L 2 (Ω). Note that assumption 1 holds true in many cases, e.g. for tight frames this choice of W is standard since it leads to W = c -1 W * , where c ∈ (0, ∞) is a constant.

Let us notice that the minimizer of F y is not necessarily unique. Given y, we denote by G y the set of all minimizers of F y :

G y def = x ∈ ℓ 2 (I) : F y (x) = min x F y (x) . (25)
For every y TH , the set G y is convex and corresponds to images û = W x which are visually very similar, as stated in the theorem below whose proof can be found in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. Some guidelines for the choice of λ i were investigated in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. The conclusions can be summarized as follows.

(a) If i ∈ I 1 , the parameter λ i should be close to, but below the upper bound w i TV , since above this bound, outliers y[i] cannot be restored.

(b) For i ∈ I 0 , a reasonable choice is

λ i = max k =i Ω (∇ w i (ξ)) T ∇ w k (ξ) |∇ w k (ξ)| dξ ,
where . T stands for transposed. If λ i is below this bound, some neighboring coefficients that are erroneously set to zero might not be restored correctly, even though Gibbs-like oscillations are reduced.

Another important remark is that, for some multiscale transforms, the bounds discussed above are constant. This is for instance the case for the wavelet transform.

For the frame W = w i : i ∈ I we here focus on the second generation curvelet transform because of the following facts.

• This transform is known to provide near-optimal non-adaptive sparse representation of piecewise smooth images away from smooth edges. Thus it is a very good candidate to capture efficiently the geometrical content of the log-data.

• The bounds on (λ i ) i∈I are nearly constant, so we use only two values for λ i ,

λ i = λ0 > 0, ∀i ∈ I 0 and λ i = λ1 > 0, ∀i ∈ I 1 .
Then the diagonal matrix Λ in (23

) satisfies Λ[i, i] ∈ { λ0 , λ1 }.
• This curvelet transform corresponds to a tight frame which will turn out to be helpful for the subsequent optimization scheme, see in particular § 3.2.2.

Minimization scheme

Let us rewrite the minimization problem defined by ( 15), ( 19) and ( 22)-( 24) in a more compact form:

find x such that F y (x) = min x F y (x) for F y (x) = Ψ(x) + Φ(x) . (26)
Clearly, Ψ in [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF] and Φ in [START_REF] Coifman | Combining the calculus of variations and wavelets for image enhancement[END_REF] are proper lower-semicontinuous convex functions, hence the same holds true for F y . The set G y introduced in ( 25) is non-empty by Theorem 1 and can be rewritten as

G y = {x ∈ ℓ 2 (I) x ∈ (∂F y ) -1 (0)},
where ∂F y stands for subdifferential operator. Minimizing F y amounts to solving the inclusion

0 ∈ ∂F y (x) ,
or equivalently, to finding a solution to the fixed point equation

x = (Id + γ∂F y ) -1 (x) , ∀γ > 0, (27)
where (Id + γ∂F y ) -1 is the resolvent operator associated to ∂F y and Id is the identity map on the Hilbert space ℓ 2 (I). The schematic algorithm resulting from ( 27), namely

x (k+1) = (Id + γ∂F y ) -1 (x (k) ), (28)
is a fundamental tool for finding the root of any maximal monotone operator [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF], such as e.g. the subdifferential of a convex function. By [START_REF] Darbon | The use of levelable regularization functions for MRF restoration of SAR images while preserving reflectivity[END_REF], γ > 0 can be seen as the stepsize of the algorithm.

Unfortunately, the resolvent operator (Id + γ∂F y ) -1 for F y in ( 26) cannot be calculated in closed-form;

we focus on iterative splitting methods.

Generalities on splitting methods

Splitting methods do not attempt to evaluate (28) directly; instead, they are based on separate evaluations of the resolvent operators (Id + γ∂Ψ) -1 and (Id + γ∂Φ) -1 . The latter are usually easier and this turns out to be true for our functionals Ψ and Φ in [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF].

Splitting methods for monotone operators have numerous applications for convex optimization and monotone variational inequalities. Even though the literature is abundant, these can basically be systematized into three main classes: the forward-backward [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF][START_REF] Tseng | Applications of a splitting algorithm to decomposition in convex programming and variational inequalities[END_REF][START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF], the Douglas/Peaceman-Rachford [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF],

and the little-used double-backward [START_REF] Lions | Une méthode itérative de resolution d'une inéquation variationnelle[END_REF][START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in hilbert space[END_REF]. A recent theoretical overview of all these methods can be found in [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF][START_REF] Eckstein | A family of projective splitting methods for the sum of two maximal monotone operators[END_REF]. Forward-backward can be seen as a generalization of the classical gradient projection method for constrained convex optimization, hence it inherits all its restrictions. Typically, one must assume that either Ψ or Φ is differentiable with Lipschitz continuous gradient, and the stepsizes γ must fall in a range dictated by the gradient modulus of continuity; see [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] for an excellent account. Since both Ψ and Φ are non differentiable, forward-backward splitting is not adapted to our criterion [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF].

We will focus on a Douglas/Peaceman-Rachford approach since differentiability of neither of the functions Ψ or Φ is required. The derivation of our algorithm relies on the calculation of the proximity operators as defined next.

Proximity operators

Proximity operators were inaugurated in [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] as a generalization of convex projection operators.

Definition 1 (Moreau [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]) Let ϕ : H → R be a proper, lower-semicontinuous and convex function where H is a Hilbert space. Then, for every z ∈ H, the function h → ϕ(h) + zh 2 2 /2, for h ∈ H, achieves its infimum at a unique point denoted by prox ϕ z, i.e.

prox ϕ z = argmin h∈H ϕ(h) + 1 2 z -h 2 2 . ( 29
)
The operator prox ϕ : H → H thus defined is called the proximity operator of ϕ.

It will be convenient in the sequel to introduce the reflection operator

rprox ϕ def = 2prox ϕ -Id. ( 30
)
By the minimality condition for [START_REF] Denis | SAR Image Regularization With Fast Approximate Discrete Minimization[END_REF], it is straightforward that ∀z, p ∈ H we have

p = prox ϕ z ⇐⇒ z -p ∈ ∂ϕ(p). (31)
Then, for any ϕ, we have (Id + ∂ϕ) -1 = prox ϕ , and we can write for all γ > 0, (Id + γ∂Ψ) -1 = prox γΨ and (Id + γ∂Φ) -1 = prox γΦ .

One can note that (28) can also be written as x (k+1) = prox γFy (x (k) ) for γ > 0.

Douglas-Rachford splitting algorithms

The Douglas/Peaceman-Rachford family is the most general preexisting class of maximal monotone operator splitting methods. Given a fixed scalar γ > 0 and a sequence µ k ∈ (0, 2), this class of methods can be expressed via the following recursion written in the compact form

x (k+1) = 1 - µ k 2 Id + µ k 2 (2prox γΨ -Id) • (2prox γΦ -Id) x (k) , = 1 - µ k 2 Id + µ k 2 rprox γΨ • rprox γΦ x (k) . (33)
Since our problem (26) admits solutions, the following result ensures that iteration (33) converges for our criterion F y .

Theorem 3 Let γ > 0 and µ k ∈ (0, 2) be such that k∈N µ k (2µ k ) = +∞. Take x (0) ∈ ℓ 2 (I) and consider the sequence of iterates defined by [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF]. Then, (x (k) ) k∈N converges weakly to some point x ∈ ℓ 2 (I)

and prox γΦ x ∈ G y .
This theorem is a straightforward consequence of [25, Corollary 5.2]. For instance, the sequence µ k = 1, ∀k ∈ N, satisfies the requirement of the theorem.

According to this theorem, the restored log-image x, as defined by ( 19) and ( 22)-( 24), reads

x = prox γΦ x, ( 34
)
where x is the output of the Douglas-Rachford algorithm (33).

Douglas-Rachford splitting algorithm for our criterion

To implement the Douglas-Rachford recursion (33) that solves our minimization problem [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF], we need to compute the proximity operators prox γΨ and prox γΦ . This is detailed in this subsection.

Proximity operator of Ψ

The proximity operator of γΨ is established in the lemma stated below.

Lemma 1 Let x ∈ ℓ 2 (I). Then

prox γΨ (x)[i] = y TH [i] + T S γλi (x[i] -y TH [i]) , ∀i ∈ I , (35)
where T S θ is the soft-thresholding operator with threshold θ

T S θ (t) = 0 if |t| ≤ θ, t -θsign(t) otherwise, t ∈ R. (36)
Proof. Ψ as given in ( 22) is an additive separable function in each coordinate i ∈ I. Thus, solving the proximal minimization problem of Definition 1 is also separable and amounts to solving independently #I scalar minimization problems.

For any convex function ϕ and t ∈ R, s ∈ R, let ψ(s) = ϕ(st). Thus, using [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF], it is straightfroward

to show that p = prox ψ (s) ⇐⇒ p = t + prox ϕ (s -t) . (37)
Furthermore, it is easy to find that for ϕ(h

) = θ |h| prox ϕ (s) = T S θ (s). (38)
Applying (37) to ψ(s) = θ |s -t| with θ = γλ i , s = x[i] and t = y TH [i] yields [START_REF] Eckstein | A family of projective splitting methods for the sum of two maximal monotone operators[END_REF].

Note that now

rprox γΨ (x) = 2 y TH [i] + T S γλi (x[i] -y TH [i]) i∈I -x . (39)

Proximity operator of Φ

Clearly, Φ(x) = • TV • W (x) is a pre-composition of the TV-norm with the linear operator W . However, computing the proximity operator of Φ for an arbitrary W may be intractable. We then systematically assume that our frame is tight with constant c ∈ (0, ∞); i.e. W * W = c Id. An immediate consequence

is that W = c -1 W * .
The precise definition of the TV norm that we actually use depends on W x which belongs to L 2 (Ω).

Let X = L 2 (Ω) × L 2 (Ω) ⊂ L 2 (R 2 ) × L 2 (R 2
) and •, • X be the inner product in X . For any q ∈ [1, ∞], we denote by • q the L q norm on X . For any τ > 0, we define

B ∞ (τ ) as the closed L ∞ -ball of radius τ in X , B ∞ (τ ) def = z = (z 1 , z 2 ) ∈ X : |z(ξ)| ≤ τ, ∀ξ ∈ Ω . (40)
Let us also define the set

C(τ ) def = div(z) ∈ L 2 (Ω) : z ∈ C ∞ c (Ω × Ω) ∩ B ∞ (τ ) . (41)
More precisely the TV regularization term we use reads

u TV = sup Ω u(ξ) w(ξ) dξ : w ∈ C(1) . (42)
Let us remind that Φ(x) = W x TV . The expression of prox γΦ is given in statement (i) of the next lemma while the computation scheme to solve statement (ii) is established in Lemma 3.

Lemma 2 Let x ∈ ℓ 2 (I) and B ∞ (.) be as defined in [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF].

(i) We have

prox γΦ (x) = Id -W • Id -prox c -1 γ • TV • W (x) ; (43)
(ii) Furthermore,

prox c -1 γ • TV (u) = u -P C(c -1 γ) (u) , ∀u ∈ L 2 (Ω) , (44)
where C(.) is defined in (41) and P C(.) is the orthogonal projector on C(.).

Proof. Since W is bounded and linear, • TV is continuous and convex, and there is u ∈ L 2 (Ω) such that u TV < ∞, it is clear that all assumptions required in [26, Proposition 11] are satisfied. Applying the same proposition yields statement (i).

We focus next on (ii). Note that for any τ > 0, C(τ ) in ( 41) is a closed convex subset since B ∞ (τ ) is closed and convex, and the operator div is linear; thus the projection P C(τ ) is well defined.

Let us remind that the Legendre-Fenchel (known also as the convex-conjugate) transform of a function ϕ : H → R, ϕ ≡ ∞, where H is an Hilbert space, is defined by

ϕ ⋆ (h) = sup z∈dom(ϕ) h, z -ϕ(z) ,
One can see also [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]Lemma 2.10] for an alternate proof of [START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF].

Let ı S denote the indicator function of a nonempty set S, i. e.

ı S (z) = 0 if z ∈ S , +∞ otherwise .
It is easy to see that from [START_REF] Krissian | Oriented speckle reducing anisotropic diffusion[END_REF], [START_REF] Lions | Une méthode itérative de resolution d'une inéquation variationnelle[END_REF] and the definition of the conjugate that

c -1 γ • TV ⋆ (z) = ı C(c -1 γ) .
On the other hand, by Definition 1, it is easy to check that

P C(c -1 γ) = prox ı C(c -1 γ) .
Combining the last two equations yields prox c -1 γ . TV

⋆ = P C(c -1 γ) .
Identifying c -1 γ . TV with ϕ and c -1 γ . TV ⋆ with ϕ ⋆ , equation ( 45) leads to statement (ii). The proof is complete.

Note that our argument [START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF] for the computation of prox c -1 γ • TV (u) is not used in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF], which instead uses conjugates and bi-conjugates of the objective function.

Remark 3 In view of ( 44) and ( 41), one can see that the term between the middle parentheses in (43) admits a simpler form:

Id -prox c -1 γ • TV = P C(c -1 γ) .
Using [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] along with ( 43)-( 44) we easily find that

rprox γΦ (x) = Id -2W • P C(c -1 γ) • W (x) . ( 46
)
Calculation of the projection P C(.) in (44) in a discrete setting. In what follows, we work in the discrete setting. We consider that W ∈ R M ×N is the analysis matrix associated to a tight frame with

N < M = #I < ∞. Note that now W * = W T and thus W = c -1 W T .
Next we replace X by its discrete counterpart,

X = ℓ 2 (Ω) × ℓ 2 (Ω) where #Ω = N. (47)
We denote the discrete gradient by ∇ and consider Div : X → ℓ 2 (Ω) the discrete divergence defined by analogy with the continuous setting 2 as the adjoint of the gradient Div = -∇ * ; see e.g. [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF].

Unfortunately, the projection in [START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF] does not admit an explicit form. The next lemma provides an iterative scheme to compute the proximal points introduced in Lemma 2. In this discrete setting, the set C(.) introduced in (41) admits a simpler expression:

C γ c = Div(z) ∈ ℓ 2 (Ω) z ∈ B ∞ γ c , (48)
where B ∞ (.) is defined according to (40).

2 More precisely, let u ∈ ℓ 2 (Ω) be of size m × n, N = mn. We write

( ∇u)[i, j] = `u[i + 1, j] -u[i, j], u[i, j + 1] -u[i, j] ẃith boundary conditions u[m + 1, i] = u[m, i], ∀i and u[i, n + 1] = u[i, n], ∀i; then for z ∈ X , we have (Div(z))[i, j] = `z1 [i, j] -z 1 [i -1, j] ´+ `z2 [i, j] -z 2 [i, j -1] álong with z 1 [0, i] = z 1 [m, i] = z 2 [i, 0] = z 2 [i, n] = 0, ∀i.
Lemma 3 We adapt all assumptions of Lemma 2 to the discrete setting, as explained above. Given u ∈ ℓ 2 (Ω), consider the forward-backward iteration: z (0) ∈ B ∞ (1) and, for all k ∈ N,

z (k+1) = P B∞(1) z (k) + β k ∇ Div(z (k) ) -cu/γ where 0 < inf k β k ≤ sup k β k < 1/4, (49)
and, ∀(i, j) ∈ Ω,

P B∞(1) (z)[i, j] =    z[i, j] if |z[i, j]| ≤ 1; z[i, j] |z[i, j]| otherwise . ( 50
)
Then (i) (z (k) ) k∈N converges to a point ẑ ∈ B ∞ (1); (ii) u -c -1 γDiv(z (k) ) k∈N converges to u -c -1 γDiv(ẑ) = prox c -1 γ • TV (u) at the rate O(1/k).
Proof. Given u ∈ ℓ 2 (Ω), the projection ŵ = P C(c -1 γ) (u), where C(.) is given in [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF], is unique and satisfies ŵ = argmin

w∈C(c -1 γ) 1 2 u -w 2 = argmin w∈ℓ 2 (Ω) 1 2 c γ u -w 2 subject to w = Div(z) for z ∈ B ∞ (1)
ŵ = Div(ẑ) where ẑ = argmin z∈B∞(1)

1 2 c γ u -Div(z) 2 , ( 51
)
where . denotes the Euclidian norm on ℓ 2 (Ω). This problem can be solved using a projected gradient method (which is a special instance of the forward-backward splitting scheme) whose iteration is given by [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF]. This iteration converges to a minimizer of ( 51)-see [25, Corollary 6.5]-provided that the stepsize

β k > 0 satisfies sup k β k < 2/δ 2
, where δ is the spectral norm of the Div operator. It is easy to check that δ 2 ≤ 8-see e.g. [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF]. Hence statement (i).

Next we focus on statement (ii). Set

ω (k) = cγ -1 u -Div(z (k)
), ∀k ∈ N and ω = cγ -1 u -Div(ẑ).

Let J be the dual objective given in [START_REF] Nikolova | A variational approach and impulse noise[END_REF], namely

J = H • Div + G, (52)
with

H : w → 1 2 c γ u -w 2 and G = ı B∞(1)
.

Let DH denote the the usual derivative of H. For an optimal solution ẑ ∈ B ∞ (1), we define the Bregman-like distance as the functional

B(w) = G(w) -G(ẑ) + -∇(DH)(Div(ẑ)), w -ẑ X , ∀ w ∈ X , (53)
It is obvious that B(ẑ) = 0. We also have B(w) ≥ 0, ∀w ∈ X . This is checked by noting that the minimality condition corresponding to ( 52) is equivalent to -∇(DH)(Div(ẑ)) ∈ ∂G(ẑ). Applying the subgradient inequality to G proves the non-negativity claim.

As H • Div is differentiable, the Taylor distance is the remainder of the Taylor expansion of H • Div near ẑ

T(w) = H(Div(w)) -H(Div(ẑ)) --∇(DH)(Div(ẑ)), w -ẑ X , ∀ w ∈ X . ( 54
)
Reminding that Div = -∇ * and that H is convex shows that T(w) ≥ 0, ∀w ∈ X , along with T(ẑ) = 0.

It is easy to verify that

B(z (k) ) + T(z (k) ) = J(z (k) ) -J(ẑ) , ∀k ∈ N . ( 55
)
Using [START_REF] Polyak | Introduction to Optimization[END_REF] or [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]Theorem 4], the convergence rate over J satisfies

J(z (k) ) -J(ẑ) ≤ 2δ 2 R 2 k + 2 , ∀ k ∈ N , (56)
where 55), along with the facts that B(z (k) ) ≥ 0 and that δ 2 ≤ 8 leads to

R def = max z -ẑ : J(z) ≤ J(z (0) ) for z (0) ∈ B ∞ (1). R is obviously finite since B ∞ (1) is bounded. Using (
T(z (k) ) ≤ J(z (k) ) -J(ẑ) ≤ 16R 2 k + 2 , ∀ k ∈ N . ( 57
)
On the other hand, by the definition of H, which is a quadratic (strongly convex of modulus 1), we have

T(z (k) ) = H(Div(z (k) )) -H(Div(ẑ)) --∇(DH)(Div(ẑ)), z (k) -ẑ X = H(Div(z (k) )) -H(Div(ẑ)) -(DH)(Div(ẑ)), Div(z (k) ) -Div(ẑ) = 1 2 Div(z (k) ) -Div(ẑ) 2 = 1 2 ω (k) -ω 2 . ( 58
)
Piecing together ( 58) and ( 57), we obtain

ω (k) -ω 2 = 2T(z (k) ) ≤ 32R 2 k + 2 , ∀ k ∈ N , (59)
This completes the proof.

Note that computing prox • TV amounts to solving a discrete L 2 -TV-denoising (the criterion given in [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF]).

The forward-backward splitting-based iteration proposed in [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF] to compute the proximity operator of the TV-norm is different from the projection algorithm of [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF]. A similar iteration was suggested in [START_REF] Chambolle | Total variation minimization and a class of binary MRF models[END_REF] without a proof. The forward-backward splitting allows to derive a sharper upper-bound on the stepsize β k than the one proposed in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF]-actually twice as large. What is more, our proof is simpler than the one in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] since it uses well known properties of proximity operators, and we have a convergence rate on the iterates.

Our iteration to solve this problem is one possibility among others. While this paper was submitted, our attention was drawn to an independent work of [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF] who, using a different framework, derive an iteration similar to [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF] to solve the L 2 -TV-denoising. Another parallel work of [START_REF] Zhu | Duality-based algorithms for total-variation-regularized image restoration[END_REF] propose an application of gradient projection to solving the dual problem [START_REF] Nikolova | A variational approach and impulse noise[END_REF], and the authors of [START_REF] Weiss | total variation minimization under constraints in image processing[END_REF] applied the multi-step Nesterov scheme to [START_REF] Nikolova | A variational approach and impulse noise[END_REF]. See also [START_REF] Beck | Teboulle A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF] for yet another multi-step iteration to solve [START_REF] Nikolova | A variational approach and impulse noise[END_REF]. We are of course aware of max-flow/min-cut type algorithms, for instance the one in [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF]. We have compared our whole denoising procedure using our implementation of prox • TV and the max-flow based implementation that we adapted from the code available at [17]. We obtained similar results, although the max-flowbased algorithm was faster, mainly because it uses the ℓ 1 approximation of the discrete gradient, namely

( ∇u)[i, j] 1 = u[i + 1, j] -u[i, j] + u[i, j + 1] -[i, j] .
Let us remind that this approximation for the discrete gradient does not inherit the rotational invariance property of the L 2 norm of the usual gradient;

we observed that the quality of the resultant restorations is slightly reduced, compared to those involving the L 2 norm of the discrete gradient.

Comments on the Douglas-Rachford scheme for F y

The bottleneck of the minimization algorithm is in the computation of the proximity-operator of the TVnorm (Lemma 2). However, when inserted in the whole numerical scheme, this is not a real drawback as we explain it below. A crucial property of the Douglas-Rachford splitting scheme ( 33) is its robustness to numerical errors that may occur when computing the proximity operators prox Ψ and prox Φ , see [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF].

We have deliberately omitted this property in [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF] for the sake of simplicity. This robustness property has important consequences: e.g. it allows us to run the forward-backward sub-recursion (49) only a few iterations to compute an approximate of the TV-norm proximity operator in the inner iterations, and the Douglas-Rachford is still guaranteed to converge provided that these numerical errors are under control. More precisely, let a k ∈ ℓ 2 (I) be an error term that models the inexact computation of prox γΦ in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF], as the latter is obtained through [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF]. 4 Bias correction to recover the sought-after image

Recall from (5) that u o = log S o and set û = W x as the estimator of u o , where x is a minimizer of F y obtained from the Douglas-Rachford iteration. Unfortunately, the estimator û is prone to bias, i.e.

E [û] = u o -b û.
A problem that classically arises in statistical estimation is how to correct such a bias. More importantly is how this bias affects the estimate after applying the inverse transformation, here the exponential. Our goal is then to ensure that for the estimate Ŝ of the image satisfies E Ŝ = S o .

To this end we need the expectation and the variance of the log-noise n in ( 5)- [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF]. One can prove that

E [n] = ψ 0 (K) -log K , (60) Var
[n] = ψ 1 (K), (61)
where

ψ k (z) = d dz k+1 log Γ(z) (62)
is the polygamma function [START_REF] Abramowitz | Handbook of mathematical functions[END_REF].

Expanding exp û in the neighborhood of E [û], we have

exp û = exp (E [û]) 1 + (û -E [û]) + (û -E [û]) 2 2 + R 2 ,
and therefore

E [exp û] = exp (E [û])(1 + Var [û] /2 + E [R 2 ]) = S o exp (-b û)(1 + Var [û] /2 + E [R 2 ]) , (63)
where R 2 is the Lagrange remainder in the Taylor series. One can observe that the posterior distribution of û is nearly symmetric, in which case

E [R 2 ] ≈ 0. That is, b û ≈ log(1+Var [û] /2) to ensure unbiasedness.
Consequently, finite sample (nearly) unbiased estimates of u o and S o are respectively û+log(1+Var [û] /2), and exp (û) (1 + Var [û] /2). Var [û] can be reasonably estimated by ψ 1 (K), the variance of the noise n in (5) being given in [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF]. Thus, given the restored log-image û, our restored image reads:

Ŝ = exp (û) (1 + ψ 1 (K)/2) . (64)
The authors of [START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] propose a direct estimate of the bias b û using the obvious argument that the noise n in the log-transformed image has a non-zero mean ψ 0 (K)log K. A quick study shows that the functions (1 + ψ 1 (K)/2) and exp(log Kψ 0 (K)) are very close for K reasonably large.

We should emphasize that the bias correction approach we propose can be used in a more general setting.

Full algorithm to suppress multiplicative noise

Now, piecing together Lemma 1, Lemma 2 and Theorem 3, we arrive at the multiplicative noise removal algorithm:

Task: Denoise an image S contaminated with multiplicative noise according to (3).

Parameters:

The observed noisy image S, number of iterations N DR (Douglas-Rachford outer iterations) and N FB (Forward-Backward inner iterations), stepsizes µ k ∈ (0, 2), 0 < β s < 1/4 and γ > 0, tight-frame transform W and initial threshold T (e.g. T = 2 ψ 1 (K)), regularization parameters λ 0,1 associated to the sets I 0,1 .

Specific operators:

• Soft-thresholding T S γλi in [START_REF] Froment | Artifact free signal denoising with wavelets[END_REF].

• Projector P B∞(1) in [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF].

• The discrete gradient ∇ and divergence Div.

• ψ 1 (•) defined according to (62) (built-in Matlab function, otherwise see [START_REF] Press | Numerical recipes, the art of scientific computing[END_REF]).

Initialization:

• Compute v = log S and transform coefficients y = W v. Hard-threshold y at T to get y TH . Choose an initial x (0) .

Main iteration:

For k = 0 to N DR ,

(1) Inverse curvelet transform of x (k) according to u (k) = W x (k) .

(2) Initialize z (0) ; For s = 0 to N FB -1

z (s+1) = P B 1 ∞ (X ) z (s) + β s ∇ Div(z (s) ) -c γ u (k) .
(3) Set z (k) = z (NFB) .

(4) Compute w (k) = c -1 γ Div(z (k) ).

(5) Forward curvelet transform: α (k) = W w (k) . ( 6) From ( 46) compute r (k) = rprox γΦ (x (k) ) = x (k) -2α (k) . ( 7) By [START_REF] Hiriart-Urruty | Convex analysis and Minimization Algorithms[END_REF] compute

q (k) = rprox γΨ • rprox γΦ (x (k) ) = 2 y TH [i] + T S γλi r (k) [i] -y TH [i] i∈I -r (k) .
(8) Update x (k+1) using ( 33):

x (k+1) = 1 - µ k 2 x (k) + µ k 2 q (k) .

End main iteration

Output: Using (34) and [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF], get the denoised image Ŝ = exp W (x (NDR)α (NDR) ) (1 + ψ 1 (K)/2).

Remark 4 (Computation load)

The bulk of computation of our denoising algorithm is invested in applying W and its pseudo-inverse W . These operators are of course never constructed explicitly, rather they are implemented as fast implicit analysis and synthesis operators. Each application of W or W cost O(N log N ) for the second generation curvelet transform of an N -pixel image [START_REF] Candès | Fast discrete curvelet transforms[END_REF]. If we define N DR and N FB as the number of iterations in the Douglas-Rachford algorithm and the forward-backward sub-iteration, the computational complexity of the denoising algorithm is of order N DR N FB 2N log N operations.

Experiments

In all experiments carried out in this paper, our algorithm was run using second-generation curvelet tight frame along with the following set of parameters: ∀t, µ k ≡ 1, β s = 0.24, γ = 10 and N DR = 50. For comparison purposes, some very recent multiplicative noise removal algorithms from the literature are considered. We compare our method with the most recent successful algorithms, namely:

• BS Algorithm proposed in [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF]: this is a Stein-Block thresholding method in the curvelet domain, applied on the log transformed image. It is a sophisticated shrinkage-based denoiser that thresholds the coefficients by blocks rather than individually, and has been shown to be nearly minimax over a large class of images in presence of additive bounded noise (not necessarily Gaussian nor independent). This algorithm has one threshold parameter T, and in all our experiments, we set it to the theoretical value T = 4.50524 devised in [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF].

• AA Algorithm proposed in [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF]: it is sketched in [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF]; the algorithm was applied using 1000 iterations and stepsize dt = 0.1, as recommended by the authors;

• SO Algorithm proposed in [START_REF] Shi | A nonlinear inverse scale space method for a convex multiplicative noise model[END_REF]: the first algorithm proposed in [START_REF] Shi | A nonlinear inverse scale space method for a convex multiplicative noise model[END_REF], see [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF]. As recommended in [START_REF] Shi | A nonlinear inverse scale space method for a convex multiplicative noise model[END_REF] the stopping rule is to reach

k * such that k * = max{k ∈ N : Var [u k -u o ] ≥ Var [n] = ψ 1 (K)}
where u o is the underlying log-image and n the relevant noise; see [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] for the variance. We systematically used stepsize dt = 0.001 (except for Fig. 3 where dt = 0.0005).

The denoising algorithms were tested on three images: Cameraman (N = 256 2 pixels), an image of Fields (N = 512 2 pixels) and an aerial picture of the French city Nîmes (N = 512 2 pixels). All images were normalized so that their gray-scale is in the range [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]256]. For each image, a noisy observation is generated by multiplying the original image by a realization of noise according to the model in ( 3)-( 4

Qualitative results

The results are depicted in Fig. 3-6 for the Fields image, Fig. 7-9 for the Nîmes picture, and Fig. 10-12 for Cameraman. Our denoiser does a good job at restoring faint geometrical structures of the images even for low values of K, see for instance the results on Nîmes and Fields for K = 1 and K = 4. As expected, our hybrid method is less prone to staircasing artifacts and takes benefit from the curvelet transform for capturing efficiently the geometrical content of the images. Our algorithm performs among the best and even outperforms its competitors most of the time both visually and quantitatively as revealed by the PSNR and MAE values. Note also that a systematic behavior of AA and SO methods for low values of K is their tendency to lose some important details and the persistence of a low-frequency ghost as it can be seen on the residual images on Fig. 4. For SO method, the number of iterations necessary to satisfy the stopping rule rapidly increases when K decreases.

Quantitative results

The above visual results were confirmed by Monte Carlo simulations where, for each tested image and each value of K ∈ {1, 4, 10}, ten noisy realizations were generated. Then the compared algorithms were applied to the same noisy versions. The output PSNR and MAE were averaged over the ten denoised realizations.

For fair comparison, the parameters of SO and AA were tweaked manually to reach their best performance level. For SO, α = 0.25 was recommended by the authors. Their values are summarized in Table 1. Notice that the parameters of our method are very stable with respect to the image, except for Cameraman with small K. Curvelet based methods are indeed quite inadequate to this image since the thresholding step generates Gibbs-like oscillations and outliers that are strongly visible on the nearly uniform background (see the image restored by BS method). In order to remove these artifacts, we used different values for T , λ 0 and mostly λ 1 . Despite this drawback, our method yields quite good results on Cameraman.

The denoising performance results are tabulated in Table 2 where the best PSNR and MAE value is shown in boldface. By inspection of this table, the PSNR improvement brought by our approach can be quite high particularly for K = 1 (see e.g. the Fields image) and the visual resolution is quite respectable.

This is an important achievement since in practice K has a small value, usually 1, rarely above 4. This improvement becomes less salient as K increases which is intuitively expected. But even for K = 10, the PSNR of our algorithm can be higher by ∼ 0.45dB to more than 10dB compared to AA, SO and BS methods, depending on the image.

Fields 512 × 512 Cameraman 256 × 256 Fields 512 × 512

Algorithm Parameters K = 1 K = 4 K = 10 Ours T = 2 ψ 1 (K), λ 0 = 1.2, λ 1 = 10 T = 2 ψ 1 (K), λ 0 = 1.3, λ 1 = 10 T = 2 ψ 1 (K), λ 0 = 1.3,
Algorithm Parameters K = 1 K = 4 K = 10 Ours T = 2.6 ψ 1 (K), λ 0 = 1.8, λ 1 = 5.7 T = 2.5 ψ 1 (K), λ 0 = 1.8, λ 1 = 5.7 T = 2.1 ψ 1 (K), λ 0 = 1.3, λ 1 = 10 AA ρ = 120 ρ = 240 ρ = 390 SO ρ = 0.04, α = 0.25 ρ = 0.1, α = 0.25 ρ = 1, α = 0.25

Conclusions

This work proposes quite an original, efficient and fast method for multiplicative noise removal. The latter is a difficult problem that arises in various applications relevant to active imaging system, such as laser imaging, ultrasound imaging, SAR and many others. Multiplicative noise contamination involves inherent difficulties that severely restrict the main restoration algorithms.

The main ingredients of our method are: (1) consider the log-data to restore a log-image; (2) preprocess the log-data using and under-optimal hard-thresholding of its tight frame coefficients; (3) restore the logimage using a hybrid criterion composed of an ℓ 1 data-fitting for the coefficients and a TV regularization in the log-image domain; (4) restore the sought-after image using an exponential transform along with a pertinent bias correction. The resultant algorithm is fast, its consistency and convergence are proved theoretically.

The obtained numerical results are really encouraging since they outperform the most recent methods in this field.

Figure 1 :

 1 Figure 1: Noise distributions: K = 1 in (a), (b), (c) and K = 10 in (d), (e), (f). In (c) and (f) W is a curvelet frame. Hist in (c) and (f) stands for histogram.

(a) Noisy for K = 1 (Figure 2 :

 12 Figure 2: (a) Noisy Lena for K = 1. (b) Noisy Lena obtained via averaging, see (1), for K = 4. (c)-(h) Denoising of data v shown in (b) where the curvelet transform of v are hard-thresholded according to (15)-(17) for different choices of T where (see (61)). The displayed restorations correspond to v TH = exp (v TH ), as given in (17).

(I 1 )

 1 (a) Coefficients y[i] containing essentially noise-in which case the best we can do is to keep them null, i.e. x[i] = y TH [i] = 0; (b) Coefficients y[i] which correspond to faint edges and details in u o . Since y[i] is difficult to distinguish from the noise, the relevant x[i] should be restored using the edge-preserving prior conveyed by Φ. Let us emphasize that a careful restoration must find a nonzero x[i], since otherwise x[i] = 0 would generate Gibbs-like oscillations in û. The coefficients y[i] for i ∈ I 1 are of the following two types: (a) Large coefficients which carry the main features of the sought-after function û. They verify y[i] ≈ w i , u o and can be kept intact, i.e. x[i] = y TH [i] = y[i]. (b) Coefficients which are highly contaminated by noise, characterized by |y[i]| ≫ | w i , u o |. We call them outliers: if we had x

Theorem 2 [ 32 ]

 232 Let the assumptions of Theorem 1 hold. If x1 and x2 are two minimizers of F y (i.e. x1 ∈ G y and x2 ∈ G y ), then ∇ W x1 ∝ ∇ W x2 , a.e. on Ω, i.e. W x1 and W x2 have the same level lines.In words, the images û1 = W x1 and û2 = W x2 share the same level lines, i.e. they are obtained from each other by a local change of contrast. This is usually invisible to the naked eye.

  and that ϕ ⋆ is a closed convex function. If ϕ is convex, proper and lower semi-continuous, the original Moreau decomposition [46, Proposition 4.a] tells us that prox ϕ + prox ϕ ⋆ = Id .

  If the sequence of error terms (a k ) k∈N and step-sizes (µ k ) k∈N defined in Theorem 3 obey k∈N µ k a k < +∞, then the Douglas-Rachford algorithm (33) converges [25, Corollary 6.2]. In our case, noting that the convergence rate of Lemma 3 yields a k = O(1/k), one can easily derive a rule on the number of inner iterations at each outer iteration k such that the summability condition is verified.

2 Ŝ -S o 2 2

 22 ) with the choice µ = 1 and K ∈ {1, 4, 10}. For a N -pixel noise-free image S o and its denoised version by any algorithm Ŝ, the denoising performance is measured in terms of peak signal to noise ratio (PSNR) in decibels (dB) PSNR = 10 log 10 N max S omin S o dB , where max S omin S o gives the gray-scale range of the original image, and mean absolute-

10 T 10 T

 1010 ψ 1 (K), λ 0 = 1.5, λ 1 = = 2 ψ 1 (K), λ 0 = 1.5, λ 1 = = 2 ψ 1 (K), λ 0 = 1.3, λ 1 = 10 AA ρ = 60 ρ = 120 ρ = 130 SO ρ = 0.05, α = 0.25 ρ = 0.3, α = 0.25 ρ = 1.2, α = 0.25

Figure 3 : 1 Figure 4 :

 314 Figure 3: Fields (512 × 512) for K = 1 in (3). Restorations using different methods. Parameters: SO algorithm for ρ = 0.05, α = 0.25 and dt = 0.0005; AA algorithm for ρ = 125; our algorithm T = 2 ψ 1 (K), λ 0 = 1.2, λ 1 = 10. 22

Figure 5 :

 5 Figure 5: Fields (512 × 512) for K = 4 in (3). Restorations using different methods. Parameters: SO algorithm for ρ = 0.1, α = 0.25 and dt = 0.001; AA algorithm for ρ = 400; our algorithm T = 2 ψ 1 (K), λ 0 = 1.3, λ 1 = 10. 24

Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 6789 Figure 6: Fields (512 × 512) for K = 10 in (3). Restorations using different methods. Parameters: SO algorithm for ρ = 0.8, α = 0.25 and dt = 0.001; AA algorithm for ρ = 480; our algorithm T = 2 ψ 1 (K), λ 0 = 1.3, λ 1 = 10. 25

Figure 10 :

 10 Figure 10: Cameraman (256 × 256) for K = 1 in (3). Restorations using different methods. Parameters: SO algorithm for ρ = 0.04, α = 0.25 and dt = 0.001; AA algorithm for ρ = 125; our algorithm T = 2.6 ψ 1 (K), λ 0 = 1.8, λ 1 = 5.7. 29

Figure 11 :

 11 Figure 11: Cameraman (256 × 256) for K = 4 in (3). Restorations using different methods. Parameters: SO algorithm for ρ = 0.1, α = 0.25 and dt = 0.001; AA algorithm for ρ = 125; our algorithm T = 2.5 ψ 1 (K), λ 0 = 1.8, λ 1 = 5.7. 30

Figure 12 :

 12 Figure 12: Cameraman (256 × 256) for K = 10 in (3). Restorations using different methods. Parameters: SO algorithm for ρ = 1, α = 0.25 and dt = 0.001; AA algorithm for ρ = 125; our algorithm T = 2.1 ψ 1 (K), λ 0 = 1.3, λ 1 = 10. 31

Table 1 :

 1 Parameters used in the comparison study of Table2.

	Noisy: K = 1 in (3)	SO: PSNR=9.59, MAE=196
	AA: PSNR=15.74, MAE=76.66	BS: PSNR=22.52, MAE=35.22
	Ours: PSNR=22.89, MAE=33.67	Original (512×512)

Table 2 :

 2 Average PSNR and MAE (over ten noisy realizations) to denoise different images using the SO, AA BS and our algorithm as a function of K. The standard deviation of PSNR and MAE over the ten realizations are in parentheses. The best PSNR and MAE value is shown in boldface. The parameters used for each denoising algorithm are summarized in Table1.

	PSNR

Taking µ = 1 just amounts to re-scale S by µ.

ψ 1 (K), λ 0 = 1.3, λ 1 = 10.

Acknowledgments

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grants FREEDOM and NatImages (ANR07-JCJC-0048-01 and ANR08-EMER-009). The authors would like to thank J.-F. Aujol for providing us the Matlab code of his algorithm and giving us wise advices on the choice of the parameters. They also thank J. Shi and S. Osher for the fruitful discussions on their approach and the advices on the choice of the parameters, yet again.