S Amstutz
email: samuel.amstutz@univ-avignon.fr

AND A A Novotny
email: novotny@lncc.br

N Van Goethem
email: vangoeth@sissa.it

MINIMAL PARTITIONS AND IMAGE CLASSIFICATION USING A GRADIENT-FREE PERIMETER APPROXIMATION

Keywords: 2010 Mathematics Subject Classification. 49Q10, 49Q20, 49M25, 35J05, 35J25, 65K10 Image classification, deblurring, optimal partitions, perimeter approximation

In this paper a new mathematically-founded method for the optimal partitioning of domains, with applications to the classification of greyscale and color images, is proposed. Since optimal partition problems are in general ill-posed, some regularization strategy is required. Here we regularize by a non-standard approximation of the total interface length, which does not involve the gradient of approximate characteristic functions, in contrast to the classical Modica-Mortola approximation. Instead, it involves a system of uncoupled linear partial differential equations and nevertheless shows Γ-convergence properties in appropriate function spaces. This approach leads to an alternating algorithm that ensures a decrease of the objective function at each iteration, and which always provides a partition, even during the iterations. The efficiency of this algorithm is illustrated by various numerical examples. Among them we consider binary and multilabel minimal partition problems including supervised or automatic image classification, inpainting, texture pattern identification and deblurring.

Introduction

Image processing is a huge field of research ranging from code-based algorithms to advanced mathematical tools, and where mainly two classes of problems are addressed. The first one is image restoration whose aim is to remove all effects responsible for an image degradation: noise, blur, missing parts, etc. Another family of problems can be referred to as image segmentation, where the constituents of a given image (damaged or not) are identified: it can be different colors, intensities or texture regions. Image classification is a particular form of image segmentation, where it is emphasized that the image characteristics are sought within a prescribed number of components or phases, identified by labels.

The standard greyscale image processing problem can be stated as restore and/or segment f = Aū + ν where f : Ω → [0, 1] is the observed image, ū its idealized version (the undamaged image), A is a known or unknown mask operator (blur kernel or a projection operator away from the missing parts of ū), and ν is the noise. According to the application, one wishes to find a u which is either a continuous restoration of ū or a segmented version of ū.

Image restoration and image segmentation are hard problems per se, mainly due to the fact that the minimization problem min u∈H(Ω) J(u) := Au -f H(Ω) is in general ill-posed (w.r.t. f if A is known and also w.r.t. A if the mask is unknown) in the sense that small perturbations in the data may produce unbounded variations in the solution, but also because it is not necessarily a convex problem, especially when the set H(Ω) contains discrete levels. Of course, simultaneous segmentation and restoration of a blurred and noisy image is even a harder task, in particular in the case of blind deconvolution. Let us emphasize that an objective assessment of segmentation algorithms is hardly found, essentially because there is no unique ground-truth classification of an image against which the output of an algorithm may be compared.

Several mathematical models coexist in the literature to provide as output a restoration and/or a segmentation of a given image. Let us just mention the Mumford-Shah [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] and the TV-L 2 (the so-called Rudin-Osher-Fatemi model) and TV-L 1 functionals [START_REF] Chambolle | Total variation in imaging[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. It is today widely recognized [START_REF] Chambolle | An algorithm for total variation minimization and applications. Special issue on mathematics and image analysis[END_REF][START_REF] Chambolle | Theoretical foundations and numerical methods for sparse recovery[END_REF][START_REF] Vese | A study in the BV space of a denoising-deblurring variational problem[END_REF] that in order to obtain a solution which is smooth enough while preserving the edges, one should consider the problem min u∈H(Ω) J(u) + α|Du|(Ω) where |Du|(Ω) is the total variation of u in Ω, with Du, its distributional derivative, i.e., a measure with a diffuse part identified with ∇u outside the edges and a concentrated part on the edges, and with α a weight on the total variation of the image. For specific purposes, some authors have proposed to replace the total variation by anisotropic variants [START_REF] Dobrosotskaya | A wavelet-Laplace variational technique for image deconvolution and inpainting[END_REF][START_REF] Dobrosotskaya | Wavelet analogue of the Ginzburg-Landau energy and its Γ-convergence[END_REF][START_REF] Esedoḡlu | Decomposition of images by the anisotropic Rudin-Osher-Fatemi model[END_REF]. Segmentation problems, mainly in the bi-label case, can also be addressed by methods of moving interfaces, like level-sets [START_REF] Chan | Active contours without edges[END_REF][START_REF] Getreuer | Chan-Vese segmentation[END_REF][START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] or snakes [START_REF] Alvarez | A real time morphological snakes algorithm[END_REF]. A totally different (purely discrete) approach consists in viewing the image as a graph and segmentation as a minimal cut problem [START_REF] Chung | Spectral graph theory[END_REF][START_REF] Shi | Normalize cuts and image segmentation[END_REF][START_REF] Luxburg | A tutorial on spectral clustering[END_REF].

In this paper we propose a novel mathematically-founded method for image classification purposes. More precisely, our method assigns to each pixel of a given image a label, with a prescribed upper bound on the number of labels. Typically, each label corresponds to a grey or color level. The classification can be supervised, i.e., the intensity of grey or the color for each label is fixed, or unsupervised, meaning that these values are automatically chosen by the algorithm (as done, e.g., in [START_REF] Chambolle | A convex approach to minimal partitions[END_REF][START_REF] Getreuer | Chan-Vese segmentation[END_REF]). In Fig. 1 the resulting classification of a color image is shown with and without level updates, and where the effect of the parameter α is shown: a smaller value will provide a solution with a greater perimeter, that is with more details. Let us emphasize that the value of α should be specifically adapted for each kind of image to reconstruct. In particular a "large enough" value of α should be used for noisy images, in order to remove the spurious perimeter created by the noise. As an example, an unsupervised greylevel simultaneous classification and denoising of a picture with our method (and the TV-L 1 model) is shown in Fig. 2. Later, we will also show examples of black and white image deblurring with an assumed known blur kernel (as is the case in many applications, the blur being due to diffraction, motion, zoom, of a measure apparatus etc.). We will also apply our method to perform combined deblurring and denoising and to classify an image in terms of prescribed textures.

Our method consists in a multiphase (piecewise constant) joint classification and restoration of an image, as based on

• a gradient-free approximation of the total variation in the functional setting, • an optimal partitioning algorithm.

Here, we mean by functional setting the fact that we approximate the continuous total variation, instead of its discretized counterpart on a chosen mesh. It is also meant that the approximating functionals converge, as a small parameter ε tends to zero, in some variational sense (the Γ-convergence) to the total variation. By gradient-free we mean that, as opposed to Modica-Mortola/Ambrosio-Tortorelli-type regularizations [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Jung | Multiphase image segmentation via Modica-Mortola phase transition[END_REF][START_REF] Kohn | Local minimisers and singular perturbations[END_REF][START_REF] Modica | The gradient theory of phase transitions and the minimal interface criterion[END_REF][START_REF] Modica | Un esempio di Γ --convergenza[END_REF][START_REF] Oudet | Approximation of partitions of least perimeter by Γ-convergence: around Kelvin's conjecture[END_REF], the gradient of the image does not appear in the approximating functionals. Hence, from a numerical point of view, our approach is better suited to the reconstruction of discontinuous images. It has the other feature to force discrete prescribed values. This is of course a major advantage as to the applications to optimal partition problems (see [START_REF] Aujol | Optimal partitions, regularized solutions, and application to image classification[END_REF][START_REF] Chambolle | A convex approach to minimal partitions[END_REF]), which are in some sense primary to classification problems and whose basic concepts are recalled hereafter. Consider a bounded domain Ω of R 2 , a number N ∈ N, functions g 1 , ..., g N ∈ L 1 (Ω), and a parameter α > 0. A prototype problem of minimal partition reads: min Ω1,...,ΩN

N i=1

Ωi

g i (x)dx + α 2 Per(Ω i ) , (1.1)
where the minimum is searched among all partitions (Ω 1 , ..., Ω N ) of Ω by subsets of finite perimeter. Here, Per(Ω i ) is the relative perimeter of Ω i in Ω, i.e., under regularity conditions, Per(Ω i ) = |∂Ω i ∩Ω|.

Hence we have α 2

N i=1 Per(Ω i ) = α i<j |∂Ω i ∩ ∂Ω j |.
Optimal partition problems in imaging are known to be challenging. Other methods as found in the literature are based on convexification, whose major difficulty is to construct an auxiliary cost function which does not to depart too much from the convex envelope of the original one. In particular these methods solve the optimal partition problem exactly only for two levels (see [START_REF] Chambolle | A convex approach to minimal partitions[END_REF][START_REF] Zach | Fast global labelling for real-time stereo using multiple plabe sweeps[END_REF]). Other mathematical tools, such as topological asymptotic analysis [START_REF] Auroux | A topological asymptotic analysis for the regularized grey-level image classification problem[END_REF], can be applied to address image classification.

In essence, our method is not designed for the restoration of continuous images, but it applies whenever a segmented restoration is sought, two examples of which are shown in Figs. 1 and2. It is a priori also not well-suited to process images with texture parts that must be finely identified, since the optimal partitioning will rather identify homogeneous regions. Nevertheless it can be observed that the segmented image of Fig. 2 features some fine characteristics, as e.g., the shade on Lena's hat and some residual texture on its feather tuft. We are aware that more efficient methods exist to address these restoration problems, in particular if the image to restore shows sharp edges of connected subregions [START_REF] Aujol | Structure-texture image decomposition-modeling, algorithms,and parameter selection[END_REF][START_REF] Auroux | Image processing by topological asymptotic expansion[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

Nonetheless, our method is adequate for all problems whose underlying mathematical structure is related to optimal partition. In particular it is adapted to classify an image without sharp edges (a pioneering reference on this subject is [START_REF] Chan | Active contours without edges[END_REF]), a simple example of which is shown in Fig. 10(b). Let us remark that addressing denoising for this kind of images is in general cumbersome, since meaningful portions of the image are completely disconnected without being considered as noise. Our method shows particular good behavior in this case. Another example where our method is well suited is the following. Assume that a detail of an image is occluded, i.e., is missing for some reason, and that this detail shows a structure which cannot be recovered by, e.g., harmonic expansion [START_REF] Auroux | Image processing by topological asymptotic expansion[END_REF]. Then, an optimal partition formulation of the problem allows us to reconstruct the missing information, as shown in Fig. 5 for the triple and quadruple-point examples.

The paper is organized as follows. The leading ideas are introduced in Section 2 in the context of binary minimal partition. The method is then extended to the multilabel partitioning problem in Section 3. The image classification problem is addressed in Sections 4 and 5 for greyscale and color images, respectively. Section 6 deals with the classification of images based on the analysis of anisotropic textures. The deblurring problem is discussed in Section 7.

Binary minimal partition

2.1. Motivation: binary image classification. We shall introduce our method on a simple example of image processing which consists of binary image classification. Let us begin with some definitions and notation. Let Ω be an open rectangle of R 2 . We define the set

E = L ∞ (Ω, {0, 1})
of characteristic functions in Ω and the functional F : E → R ∪ {+∞} such that

F (u) = 1 2 |Du|(Ω) if u ∈ BV (Ω, {0, 1}), +∞ otherwise. (2.1)
We recall that the total variation of u ∈ L 1 (Ω) is defined as

|Du|(Ω) = sup{ u, div ξ : ξ ∈ C 1 c (Ω), |ξ(x)| ≤ 1 ∀x ∈ Ω}, (2.2)
and u is said of bounded variation, denoted u ∈ BV (Ω), when |Du|(Ω) < ∞. Throughout we use the notation

u, v := Ω u v dx
for every pair of functions u, v having suitable regularity. If u is the characteristic function of E ⊂ Ω, |Du|(Ω) is called the relative perimeter of E in Ω, and E is said of finite perimeter when u ∈ BV (Ω).

In this case it is shown that

|Du|(Ω) = H 1 (∂ M E ∩ Ω)
, where ∂ M E is the measure theoretical (or essential) boundary of E and H 1 is the one-dimensional Hausdorff measure on R 2 . We recall (see, e.g., [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF]) that ∂ M E is a subset of the topological boundary ∂E of E, defined by

∂ M E = x ∈ R 2 : lim sup ρ→0 L 2 (B(x, ρ) ∩ E) L 2 (B(x, ρ)) > 0 and lim sup ρ→0 L 2 (B(x, ρ) \ E) L 2 (B(x, ρ)) > 0 ,
with L 2 the Lebesgue measure on R 2 . Moreover, E admits a generalized unit inner normal n E at H 1 -almost every point x ∈ ∂ M E, see again [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF].

The following calculus illustrates our method on a simple case. Let f ∈ L ∞ (Ω, [0, 1]) be a given image and c 1 , c 2 ∈ [0, 1] be two classes representing different intensities of grey. We want to approximate f by the piecewise constant image u c 1 + (1 -u) c 2 , with u ∈ E. Given p ∈ [1, +∞[, the binary image classification problem is formulated as: find a minimizer for

I := min u∈E uc 1 + (1 -u)c 2 -f p L p (Ω) + αF (u) . (2.3)
Since u is a characteristic function, we have

Ω |uc 1 + (1 -u)c 2 -f | p dx = Ω |u(c 1 -f ) + (1 -u)(c 2 -f )| p dx = Ω u|c 1 -f | p dx + Ω (1 -u)|c 2 -f | p dx.
We can then redefine I, up to an additive constant, as

I = min u∈E {I(u) := u, g + αF (u)} , (2.4
)

with g = |c 1 -f | p -|c 2 -f | p . (2.5)
More generally, we shall address the numerical solution of (2.4) for an arbitrary function g ∈ L 1 (Ω). We refer to this problem as a binary minimal partition problem.

As already said, a number of specific methods provide satisfactory solutions to this problem, or some of its variants. We have quoted level-sets [START_REF] Chan | Active contours without edges[END_REF][START_REF] Getreuer | Chan-Vese segmentation[END_REF][START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], snakes [START_REF] Alvarez | A real time morphological snakes algorithm[END_REF] and graph cutting [START_REF] Chung | Spectral graph theory[END_REF][START_REF] Shi | Normalize cuts and image segmentation[END_REF][START_REF] Luxburg | A tutorial on spectral clustering[END_REF], and we must also mention the convex reformulation of (2.4) given in Prop. 2.6 of [START_REF] Chambolle | Total variation in imaging[END_REF], which can be solved for instance by Chambolle's algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications. Special issue on mathematics and image analysis[END_REF]. The main advantage of our approach, as we will see later, is its straightforward generalization to the multiphase case as well as to the presence of a mask operator.

Approximation of the perimeter term. Let us define

Ẽ = L ∞ (Ω, [0, 1]),
the convex hull of E, and the function F : Ẽ → R ∪ {+∞} such that

F (u) = F (u) if u ∈ E, +∞ otherwise.
It is shown in [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] that a suitable approximation of F is provided by the functional Fε defined as

Fε (u) = inf v∈H 1 (Ω) ε ∇v 2 L 2 (Ω) + 1 ε v 2 L 2 (Ω) + u, 1 -2v . (2.6)
The above minimization problem is easily solved, and we have the alternative expression

Fε (u) = 1 ε 1 -L ε u, u , (2.7)
where L ε u is the (weak) solution of the boundary value problem with unknown v ∈ H 1 (Ω):

-ε 2 ∆v + v = u in Ω, ∂ n v = 0 on ∂Ω. (2.8)
Note that, on choosing classically min(0, v) as test function in the variational formulation of (2.8), then considering the function 1 -v as unknown, we obtain that L ε u ∈ Ẽ for all u ∈ Ẽ. Substituting Fε for F leads to the approximate minimization problem:

I ε = min u∈ Ẽ I ε (u) := u, g + α Fε (u) .
(2.9) 2.3. Mathematical properties. We first recall some key properties of the functional Fε (Proposition 2.1 and Theorems 2.3, 2.4 and 2.5) proved in [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF]. Nevertheless, the proof of Theorem 2.3 given in [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] relies on a result from [START_REF] Solci | Variational models for phase separation[END_REF], and does not lend itself to the generalization we will need later. Therefore we provide an original and well-adapted proof of Theorem 2.3. However, we point out that this proof is valid only in dimension 2 for a rectangular domain Ω.

Let us start by defining the functional

F ε (u) = inf v∈H 1 (Ω) ε ∇v 2 L 2 (Ω) + 1 ε v -u 2 L 2 (Ω) . (2.10)
It is straightforward to check that, for every u ∈ L 2 (Ω),

Fε (u) = F ε (u) + 1 ε u, 1 -u . (2.11)
Proposition 2.1. The function Fε : Ẽ → R is continuous on Ẽ for the weak- * topology of L ∞ (Ω), and it is the relaxation (i.e. the lower semi-continuous envelope) of the function

u ∈ Ẽ → F ε (u) if u ∈ E, +∞ if u / ∈ E. (2.
12)

The notion of Γ-convergence (see, e.g., [START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF][START_REF] Braides | Γ-convergence for beginners[END_REF][START_REF] Maso | An introduction to Γ-convergence[END_REF]) is a powerful tool of calculus of variations in function spaces. Given a metrizable space (X, d) (in our case X = Ẽ endowed with the distance induced by the L 1 norm) one would like the maps

f → inf X f and f → argmin X f
to be sequentially continuous on the space of extended real-valued functions f : X → R ∪ {+∞}. Throughout we will simply index sequences by ε, to actually represent any positive sequence of indices

(ε k ) such that lim k→+∞ ε k = 0. Definition 2.2 (Γ-convergence). Let (f ε ) be a sequence of functions f ε : X → R ∪ {+∞} and f : X → R ∪ {+∞}.
We say that f ε Γ-converges to f if and only if, for all u ∈ X, the following two conditions hold:

(1) for all sequences (

u ε ) ∈ X such that d(u ε , u) → 0 it holds f (u) ≤ lim inf ε→0 f ε (u ε ), (2) there exists a sequence (ū ε ) ∈ X such that d(ū ε , u) → 0 and f (u) ≥ lim sup ε→0 f ε (ū ε ).
Theorem 2.3 (Γ-convergence of the approximating functionals). As ε → 0, the functionals Fε Γ-converge to F in Ẽ, strongly in L 1 (Ω).

Theorem 2.4 (Equicoercivity). Let (u ε ) be a sequence of functions of Ẽ such that sup ε Fε (u ε ) < +∞.

There exists u ∈ E such that u ε → u strongly in L 1 (Ω) for a subsequence.

Theorem 2.5. Let u ε be an approximate minimizer of (2.9), i.e.

I ε (u ε ) = u ε , g + α Fε (u ε ) ≤ I ε + λ ε ,
with u ε ∈ Ẽ and lim ε→0 λ ε = 0. Then we have

I ε (u ε ) → I.
Moreover, the sequence (u ε ) admits cluster points, and each of these cluster points is a minimizer of (2.4).

For multilabel problems, we will need special recovery sequences in the lim sup inequality of the Γ-convergence, which can be put together in order to form a partition of Ω. In fact, we will show a pointwise convergence result (Theorem 2.8). To do so, we will establish two preliminary results. The first one provides the Γ -lim inf inequality of the Γ-convergence, and has already been proved in [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] using [START_REF] Solci | Variational models for phase separation[END_REF]. For completeness we propose here a new, self-contained, and much shorter proof.

Proposition 2.6. Let u ∈ Ẽ and (u ε ) be a sequence of Ẽ such that u ε → u strongly in L 1 (Ω). Then we have lim inf ε→0 Fε (u ε ) ≥ F (u). Proof. First, from (2.11), we have that

Fε (u ε ) ≥ 1 ε u ε , 1 -u ε . Since 0 ≤ u ε ≤ 1, we also have u ε → u in L 2 (Ω)
, whereby lim inf ε→0 Fε (u ε ) = +∞ whenever u / ∈ E. We now assume that u ∈ E and that lim inf ε→0 Fε (u ε ) = ℓ < +∞. Up to extracting a (non-relabeled) subsequence, we further assume that ℓ = lim ε→0 Fε (u ε ). Like in [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] we introduce the potentials

W (s, t) = t 2 + s(1 -2t), W(t) = t 2 + min(0, 1 -2t), which clearly satisfy W (s, t) ≥ W(t) for all (s, t) ∈ [0, 1] × [0, 1] and W(1/2 -τ ) = W(1/2 + τ ) ≥ 0 for all τ ∈ R. Set v ε = L ε u ε . Then we have Fε (u ε ) = Ω ε|∇v ε | 2 + 1 ε W (u ε , v ε ) dx ≥ Ω ε|∇v ε | 2 + 1 ε W(v ε ) dx.
The elementary Young inequality entails

Fε (u ε ) ≥ 2 Ω |∇v ε | W(v ε )dx.
Let ψ : R → R be a function such that ψ ′ (t) = W(t) for all t ∈ R, and set

w ε = ψ • v ε . We obtain by the chain rule Fε (u ε ) ≥ 2 Ω |∇w ε |dx. (2.13) Next, from Fε (u ε ) ≥ F ε (u ε ) ≥ 1 ε u ε -v ε 2 L 2 (Ω) , we infer that lim ε→0 u ε -v ε L 2 (Ω) = 0.
Obviously, the same limit holds for the

L 1 norm. It follows that v ε -u L 1 (Ω) ≤ v ε -u ε L 1 (Ω) + u ε -u L 1 (Ω) → 0. Using that 0 ≤ W(t) ≤ 1/4 for all t ∈ [0, 1] we obtain that ψ is 1/2-Lipschitz on [0, 1], which implies that w ε → w := ψ • u in L 1 (Ω).
In view of (2.13), the lower semicontinuity of the total variation yields that w ∈ BV (Ω) with lim inf ε→0 Fε (u ε ) ≥ 2|Dw|(Ω). Yet, u ∈ E and w = ψ • u entail w ∈ BV (Ω, {α, β}) with α = ψ(0) and β = ψ(1). In addition, as u takes only 0 -1 values a.e. in Ω, we can write w = α + (β -α)u and u = (w -α)/(β -α), whereby u ∈ BV (Ω, {0, 1}). It also follows that Dw = (β -α)Du, hence |Dw|(Ω) = (β -α)|Du|(Ω). We finally compute

β -α = 1 0 ψ ′ (t)dt = 1 0 W(t)dt = 2 1/2 0 tdt = 1 4 , which leads to lim inf ε→0 Fε (u ε ) ≥ 1 2 |Du|(Ω) = F (u).
The proof is thus achieved.

The following proposition gives a uniform upper bound for the functionals F ε .

Proposition 2.7. For all u ∈ BV (Ω, [0, 1]) and all ε > 0 we have

F ε (u) ≤ 1 2 |Du|(Ω).
Proof. We proceed in two steps.

Step 1. Let us first assume that u ∈ Ẽ ∩ H 1 (Ω). From (2.7) and (2.11) we obtain that

F ε (u) = 1 ε u -v ε , u , (2.14)
with v ε = L ε u. The weak formulation of (2.8) gives the alternative expression

F ε (u) = ε Ω ∇v ε .∇udx. (2.15)
Let ũ be the extension of u to the whole R 2 obtained by performing successive reflections with respect to the main axes. By uniqueness we have v ε = Φ ε * ũ, where Φ ε is the fundamental solution of the operator -ε 2 ∆ + I, i.e., the solution in the sense of distributions of

-ε 2 ∆Φ ε + Φ ε = δ. (2.16) A change of variable in (2.16) indicates that Φ ε (x) = ε -2 Φ 1 (ε -1 x), whereby ∇v ε (x) = (∇Φ ε * ũ)(x) = 1 ε 3 R 2 ∇Φ 1 ( y ε )ũ(x -y)dy.
Then a change of variable in the above integral entails

ε∇v ε (x) = R 2 ∇Φ 1 (z)ũ(x -εz)dz.
Let ν be an arbitrary unit vector of R 2 . Using that 0 ≤ ũ ≤ 1 we derive that

ε∇v ε (x).ν = R 2 (∇Φ 1 (z).ν) ũ(x -εz)dz ≤ R 2
max(0, ∇Φ 1 (z).ν)dz =: M ν .

(2.17)

We shall now compute this latter bound. The fundamental solution of the normalized screened Poisson equation is known to be given by

Φ 1 (z) = 1 2π K 0 (|z|),
where K 0 is the modified Bessel function (see [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]) defined by

K 0 (r) = +∞ 0 exp(-r cosh t)dt.
Using polar coordinates with ν as first basis vector, we arrive at

M ν = 1 2π +∞ 0 2π 0 max(0, K ′ 0 (r) cos θ)rdθ dr.
Yet, we have

K ′ 0 (r) = - +∞ 0 exp(-r cosh t) cosh t dt ≤ 0, hence M ν = 1 2π +∞ 0 3π 2 π 2 cos θdθ K ′ 0 (r)rdr = - 1 π +∞ 0 K ′ 0 (r)rdr.
Using that K 0 has an exponential decay at +∞ and a logarithmic growth at 0, an integration by parts yields

M ν = 1 π +∞ 0 K 0 (r)dr = 1 π +∞ 0 +∞ 0 exp(-r cosh t)dt dr.
By Fubini's theorem, the order of integration can be changed, which leads to

M ν = 1 π +∞ 0 1 cosh t dt = 1 2 ,
this latter integral being easily computed by the change of variable s = e t . Coming back to (2.17), which holds true for any unit vector ν and a.e. x ∈ Ω, we infer

ε∇v ε L ∞ (Ω) ≤ 1 2 .
Finally, (2.15) yields

F ε (u) ≤ 1 2 ∇u L 1 (Ω) = 1 2 |Du|(Ω).
Step 2. We now consider an arbitrary function u ∈ BV (Ω, [0, 1]). By density of C ∞ (Ω) in BV (Ω) for the intermediate convergence (see, e.g., [START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF] Theorem 10.1.2), there exists a sequence of functions

u k ∈ C ∞ (Ω) such that u k → u in L 1 (Ω) and |Du k |(Ω) → |Du|(Ω). Set ūk = P [0,1] (u k ), with P [0,1] (t) = max(0, min(1, t)) the projection onto [0, 1]. As P [0,1] is 1-Lipschitz we have ūk -u L 1 (Ω) = P [0,1] (u k ) -P [0,1] (u) L 1 (Ω) ≤ u k -u L 1 (Ω)
, and since u k ∈ H 1 (Ω) we also have

|D ūk |(Ω) = Ω |∇ū k |dx = Ω |∇u k |χ {0≤u k ≤1} dx ≤ Ω |∇u k |dx = |Du k |(Ω).
As ūk ∈ Ẽ ∩ H 1 (Ω), we get from step 1

F ε (ū k ) ≤ 1 2 |D ūk |(Ω) ≤ 1 2 |Du k |(Ω). (2.18)
From the expression (2.14) and the continuity of

L ε : L 2 (Ω) → L 2 (Ω), it stems that F ε is continuous in L 2 (Ω), thus also in L 1 (Ω, [0, 1]
). Passing to the limit in (2.18) entails

F ε (u) ≤ 1 2 |Du|(Ω),
and the proof is complete.

We now arrive at the desired pointwise convergence result.

Theorem 2.8. For all u ∈ Ẽ it holds lim ε→0 Fε (u) = F (u). Proof. In view of Proposition 2.6 applied to the constant sequence u ε = u, it suffices to prove that lim sup ε→0 Fε (u) ≤ F (u).

(2.19)

If u / ∈ BV (Ω, {0, 1}), then (2.19) is obviously satisfied since F (u) = +∞. If now u ∈ BV (Ω, {0, 1}), then Fε (u) = F ε (u) for each ε > 0, thus (2.19) is a direct consequence of Proposition 2.7.
Note that, as announced in the beginning of this section, Theorem 2.3 is immediately retrieved by combining Proposition 2.6 and Theorem 2.8.

Finally, by the direct method of the calculus of variations, we straightforwardly prove the following existence result (see [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] for details).

Proposition 2.9. The infima of I and I ε in (2.4) and (2.9) are finite and attained in BV (Ω, {0, 1}) and Ẽ, respectively.

2.4. Algorithm. Our algorithm is based on a continuation method, namely we construct a decreasing sequence (ε m ) of positive numbers tending to zero, and, for each ε = ε m , we find an approximate minimizer of (2.9) using as initialization the solution obtained at iteration m -1. In the sequel, the subscript m will be dropped for simplicity.

2.4.1. Description of the algorithm in the function space setting. Plugging (2.6) into (2.9), we obtain that the subproblem at ε fixed consists is solving the following two-level minimization problem

I ε = min u∈ Ẽ inf v∈H 1 (Ω) u, g + α ε ∇v 2 L 2 (Ω) + 1 ε v 2 L 2 (Ω) + u, 1 -2v . (2.20)
The simple structure of this problem with respect to each variable u and v leads us to use an alternating minimization algorithm. We refer to [START_REF] Allaire | Shape optimization by the homogenization method[END_REF][START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] for efficient applications of such algorithms to problems of similar form in topology optimization. The superscript k is used to designate variables computed at iteration k. The iteration k, k ≥ 1, consists in the two steps described below.

(1) The minimization with respect to v is straightforward. It consists in solving the boundary value problem

-ε 2 ∆v k + v k = u k-1 in Ω, ∂ n v k = 0 on ∂Ω. (2.21)
(2) The minimization with respect to u is a linear programming problem in a convex set. Therefore a minimizer can always be found among the extreme points of Ẽ. More precisely here, we have to minimize at each point x ∈ Ω the linear function

φ(s) = s g(x) + α ε s(1 -2v k (x)) over s ∈ [0, 1]. Setting ζ k (x) = g(x) + α ε (1 -2v k (x)) ,
a solution is immediately found as

u k (x) = 0 if ζ k (x) ≥ 0, 1 otherwise.
In other terms, u k is the characteristic function of the level-set {ζ k < 0}, denoted by

u k = χ {ζ k <0} .
This algorithm ensures a decrease of the objective function at each iteration. Moreover, each cluster point (in the weak- * topology for u and the H 1 norm topology for v) is a stationary point. Of course, as the coupled problem in (u, v) is not convex, local minimizers are theoretically not excluded. An outstanding feature of this algorithm is that u k is always a characteristic function during the iterations.

Discrete version. For solving the boundary value problem (2.8) we use finite elements on a

Cartesian mesh with Q1 shape functions. The mesh nodes coincide with the image pixels, and without any loss of generality, the mesh size is fixed to 1. The discrete counterparts of the variables u and v are therefore vectors of R n where n is the number of pixels. Denoting by K and M the stiffness and mass matrices for the (negative) Laplacian, the discrete problem reads

I = min u∈R n min v∈R n M u • g + α εKv • v + 1 ε (M v • v + M u • (1 -2v)) .
In the above expression, the dot stands for the standard dot product of R n , and 1 = (1, ..., 1) T . In this framework the two steps of the algorithm consist of solving the linear system

(ε 2 K + M )v k = u k-1 , (2.22)
and setting

ζ k = M g + α ε (1 -2v k ) , (2.23
)

u k = χ {ζ k <0} .
(2.24) This algorithm enjoys the same monotonicity and convergence properties as its infinite-dimensional counterpart. The stopping criterion we choose corresponds to a relative variation of the vector u (in squared ℓ 2 norm) between two successive iterations smaller than some threshold, fixed to 10 -5 .

The linear system (2.22) is solved in an efficient way with the help of the fast Fourier transform (FFT), according to the following procedure. First, symmetries of the image u are performed in both axial directions, in such a way that we consider a domain Ω with double width and height. Then periodicity conditions are assumed on the boundary of Ω, which is a convenient way of implementing Neumann boundary conditions. In this framework, the matrix products Kv and M v represent Then we divide ε by two at each iteration of an outer loop, that is, we choose ε m = ε max /2 m . In order to approximate (2.21) properly, ε must not be taken significantly smaller than the grid size. Thus we stop the algorithm as soon as ε m ≤ ε min = 0.1. In fact, numerical tests show that almost no more evolution occurs when ε goes below this value.

2.4.4. Initialization. The initialization of u is performed by the expressions (2.23)-(2.24) with ε → +∞, that is, we set ζ 0 = M g, u 0 = χ {ζ 0 <0} .

Numerical examples.

In Figure 3, we present an example of binary image classification using the above procedure. The function g is defined by (2.5) with the two levels c 1 and c 2 corresponding to black and white. Since this particular example aims at denoising, p = 1 is chosen in (2.3)-(2.5).

In Figure 4, we present another example, also with black and white levels and L 1 norm. Intermediate iterations are displayed.

Multilabel minimal partition

In this section we extend the previous approach to the multilabel problem. We define the set

E N = (u 1 , ..., u N ) ∈ E N , N i=1 u i = 1 .
This is the set of N -tuples of characteristic functions (u 1 , ..., u N ) = (χ Ω1 , ..., χ ΩN ), where (Ω 1 , ..., Ω N ) form a partition of Ω. Given functions g 1 , ..., g N ∈ L 1 (Ω), the minimal partition problem (1.1) reads

I := min (u1,...,uN )∈EN I(u) := N i=1 u i , g i + αF (u i ) . (3.1)
In the same way as in Section 2, this problem is approximated by

I ε := min (u1,...,uN )∈ ẼN I ε (u) := N i=1 u i , g i + α Fε (u i ) (3.2) with ẼN = (u 1 , ..., u N ) ∈ ẼN , N i=1 u i = 1 .
The following result is an extension of its binary counterpart presented in Section 2, namely Theorem 2.5.

Theorem 3.1. Let u ε = (u ε i ) 1≤i≤N be an approximate minimizer of (3.2), i.e.,

I ε (u ε ) = N i=1 u ε i , g i + α Fε (u ε i ) ≤ I ε + λ ε ,
with u ε i ∈ ẼN and lim ε→0 λ ε = 0. Then we have

I ε (u ε ) → I.
Moreover, the sequence (u ε ) admits cluster points, and each of these cluster points is a minimizer of (3.1).

Proof. As usual for this kind of results, the proof relies on the Γ-convergence of the functionals u ∈ ẼN → N i=1 Fε (u i ) along with the compactness of the sequences of approximate minimizers (see, e.g., Theorem 12.1.1 of [START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF]). The lim inf inequality of the Γ-convergence immediately passes to the sum. As to the lim sup inequality, one has to construct a recovery sequence which belongs to ẼN , which is not automatically achieved by gathering independent recovery sequences for each variable u i . However, the pointwise convergence result of Theorem 2.8 shows that the constant sequence is a trivial recovery sequence.

From this latter property, there exists ūε ∈ ẼN such that ūε → (0, ..., 0, 1) in L 1 (Ω) N and N i=1 Fε (ū ε i ) → 0. Let (u ε ) be a sequence of approximate minimizers of (3.2). Using

N i=1 u ε i , g i + α Fε (u ε i ) ≤ N i=1 ūε i , g i + α Fε (ū ε i ) + λ ε
we derive that Fε (u ε i ) is uniformly bounded for every i. The equicoercivity property of Theorem 2.4 implies that, for a subsequence,

u ε i → u i ∈ E in L 1 (Ω) for each i. As u ε ∈ ẼN we infer that u = (u 1 , ..., u N ) ∈ E N .
3.1. Algorithm. For ε fixed we have to solve the approximate problem

I ε = min (u1,...,uN )∈ ẼN inf (v1,...,vN )∈H 1 (Ω) N N i=1 u i , g i + α ε ∇v i 2 L 2 (Ω) + 1 ε v i 2 L 2 (Ω) + u i , 1 -2v i .
We use again an alternating minimization algorithm with respect to the two N -tuples of variables (u 1 , ..., u N ) and (v 1 , ..., v N ). The superscript k is again used to designate these vectors at iteration k.

(1) The minimization with respect to (v 1 , ..., v N ) consists in solving the N boundary value problems

-ε 2 ∆v k i + v k i = u k-1 i in Ω, ∂ n v k i = 0 on ∂Ω. (3.3)
(2) The minimization with respect to (u 1 , ..., u N ) is a linear programming problem. Minimizers can be found by exploring the vertices of the polyhedron ẼN , that is, E N . The practical procedure is the following. Set

ζ k i = g i + α ε (1 -2v k i ).
At each point x ∈ Ω we find an index i(x) such that

ζ k i(x) (x) = min{ζ k 1 (x), ..., ζ k N (x)}. (3.4)
We then set

u k i (x) = 1 if i = i(x)
0 otherwise. The discrete counterpart is easily obtained in the same fashion as in the binary case. Again, an outstanding feature of this algorithm is that the functions u k 1 , ..., u k N are always characteristic functions of a partition of Ω.

Numerical validation.

The three examples of Fig. 5 are taken from and compared to [START_REF] Chambolle | A convex approach to minimal partitions[END_REF], where the authors introduce a special approximate convex envelope of (3.1), called local convex envelope. Let E 0 , E 1 , ..., E N be a given partition of Ω. We define g i , i = 1, ..., N , by

g i = 0≤j≤N j =i χ Ej = 1 -χ Ei .
This means that, in the set E i , i ≥ 1, the label i is favored, whereas in the set E 0 there is no preference, or, said otherwise, no information on which label to choose. From this point of view the problem is related to the wide field of image inpainting, see e.g. [START_REF] Dobrosotskaya | A wavelet-Laplace variational technique for image deconvolution and inpainting[END_REF] and the references therein. In the subsequent examples the domain Ω is the unit square containing 376 × 357 pixels, and we choose α = 20 so as to use the same data as in [START_REF] Chambolle | A convex approach to minimal partitions[END_REF].

3.2.1.

Example with 3 labels. The partition is as shown in Fig. 5(a) (left). The set E 0 is the black disc, while each E i , i = 1, 2, 3, is assigned to a specific color, namely red, green and blue, respectively. We retrieve the triple junction, which is known to be the theoretical solution [START_REF] Chambolle | A convex approach to minimal partitions[END_REF].

3.2.2.

Examples with 4 labels. We consider now the partitions of Fig. 5(b) and 5(c), with E 0 the black disc and 4 other subsets. In these cases, the solution is not unique, and the algorithm chooses a particular one. This choice stems from the selection of a particular minimizer in (3.4). This is in contrast with the results obtained in [START_REF] Chambolle | A convex approach to minimal partitions[END_REF] which, due to the convexification method employed, are mixtures of minimizers.

Multilabel image classification

4.1. Formulation as a minimal partition problem. We come back to the classification problem presented in Section 2.1, but this time with N grey levels c 1 , ..., c N ∈ [0, 1]. We are given an image f ∈ L ∞ (Ω, [0, 1]), and consider the piecewise constant image

w = N i=1 u i c i ,
where each u i is the characteristic function of a subset Ω i of Ω such that (Ω 1 , ..., Ω N ) forms a partition of Ω. We have for any L p norm on Ω:

w -f p L p (Ω) = N i=1 u i c i -f p L p (Ω) = N i=1 u i (c i -f ) p L p (Ω) = N i=1 Ω u i |c i -f | p dx.
The difference between the piecewise constant and original images is measured by:

w -f p L p (Ω) = N i=1 u i , g i , g i = |c i -f | p .
When the levels c i are fixed, we can directly apply the algorithm of Section 3.1.

4.2.

Update of levels. However, it is often desirable to determine the grey levels within the classes automatically. Thus, we include a third step in the alternating minimization algorithm, consisting of solving min

(c1,...,cN )∈[0,1] N N i=1 u i , g i = N i=1 Ω u i |c i -f | p dx.
This problem is separable in its variables c 1 , ..., c N , and each c i must satisfy Note that, since 0 ≤ f ≤ 1, the constraint 0 ≤ c ≤ 1 can be removed. We distinguish between the two cases of practical interest, namely p = 2 and p = 1.

c i ∈ argmin c∈[0,1] Ω u i |c -f | p dx.
• p = 2. This is a standard problem which results in computing the arithmetic mean

c i = Ω u i f dx Ω u i dx .
• p = 1. Minimizing a sum of L 1 distances classically amounts to computing a median, but here, due to the weights u i this is a little more involved. Details of the procedure in the discrete setting are given in Appendix A. Note that there may be several solutions. In this case we take the half-sum of the extreme points of the minimizing set. For the initialization, the levels are equidistributed in [0, 1], i.e., we choose

c 0 i = i -1 N -1 .

4.3.

Examples. We begin with two synthetic examples taken from [START_REF] Jung | Multiphase image segmentation via Modica-Mortola phase transition[END_REF], for comparison. In [START_REF] Jung | Multiphase image segmentation via Modica-Mortola phase transition[END_REF] the authors consider a multiphase "sine-sinc" extension of the classical Modica-Mortola perimeter approximation. Their functional approximates the same T V -L 2 model as ours, except that, in essence, it weighs the length term by the label jump. To enable comparison, we choose p = 2 and adjust our parameter α so as to obtain close results. The example depicted in Fig. 6 is a 3 label classification of a noisy image, and Fig. 7 shows a 5 label classification. In each case we have represented the "reconstructed" image N i=1 u i c i but, of course, our algorithm provides the N characteristic functions u 1 , ..., u N . However the fact that the u ′ i s are indeed characteristic functions is clearly visible in the images, in contrast to [START_REF] Jung | Multiphase image segmentation via Modica-Mortola phase transition[END_REF]. Examples with the L 1 norm are shown in Figs. 8-10. In Fig. 8, an example of denoising with 3 labels is shown. In Fig. 9, the classification of the swan picture is shown with 2 and 3 labels. This result is to be compared with Fig. 4. We have chosen a smaller value of α to take into account the decrease of the fidelity term due to the update of levels. In Fig. 10 two examples with disconnected geometric elements are featured. first regularize f by solving

-∆φ + φ = f in Ω, ∂ n φ = 0 on ∂Ω, (6.1)
then we set g i = -(∇φ.ξ i ) 2 . (6.2) In practice, (6.1) is solved by finite elements, in the same way as described in Section 2.4.2. Then (6.2) is computed at each node using a discrete gradient based on finite differences. We point out that the function g i can be interpreted as the topological sensitivity of the energy functional

1 2 ∇φ 2 L 2 (Ω) + 1 2 φ -f 2 L 2 (Ω)
with respect to the creation of an insulating crack normal to ξ i , see [START_REF] Amstutz | Crack detection by the topological gradient method[END_REF] and [START_REF] Auroux | Image processing by topological asymptotic expansion[END_REF] for various applications of this concept in image processing.

We then apply the multilabel minimal partition algorithm of Section 3 with the obtained functions (g 1 , ..., g N ).

In Figure 14(a) we show an example with 2 labels. We want to detect vertical and horizontal fluctuations, therefore we choose the vectors (ξ 1 , ξ 2 ) as

ξ 1 = (1, 0) T , ξ 2 = (0, 1) T .
In Figure 14(b) we show an example with 4 labels, where the vectors ξ i are given by

ξ 1 = (1, 0) T , ξ 2 = √ 2 2 (1, 1) T , ξ 3 = (0, 1) T , ξ 4 = √ 2

Image deblurring

In this section we apply a variant of our minimal partition model for binary image deconvolution. The blurring kernel is represented by a linear and continuous operator A : L 2 (Ω) → L 2 (Ω) such that A1 = 1. The given blurred (greyscale) image is f ∈ L ∞ (Ω, [0, 1]), and the reconstructed image is w = c 1 u + c 2 (1 -u), with u ∈ E and 0 ≤ c 1 < c 2 ≤ 1. We have

Aw = c 1 Au + c 2 A(1 -u) = (c 1 -c 2 )Au + c 2 ,
and the deblurring problem reads

I := min u∈E {I(u) = Φ(u) + αF (u)} (7.1) with Φ : L 2 (Ω) → L 2 (Ω) defined by Φ(u) = (c 1 -c 2 )Au + c 2 -f 2 L 2 (Ω) .
In our model the blurring kernel A is assumed to be known, which occurs in some practical applications, like when the blur is generated by an optical device. Algorithms for blind deconvolution can be found in [START_REF] Dobrosotskaya | A wavelet-Laplace variational technique for image deconvolution and inpainting[END_REF][START_REF] Esedoglu | Blind deconvolution of bar code signals[END_REF], using the Modica-Mortola functional or wavelet-based extensions.

In our case we approximate (7.1) by

I ε = min u∈ Ẽ I ε (u) = Φ(u) + α Fε (u) , (7.2)
and use the same continuation procedure with respect to ε as before. Using that the function Φ is at the same time continuous on Ẽ endowed with the L 1 norm topology and, by convexity, lower semi-continuous on Ẽ endowed with the weak- * topology of L ∞ (Ω), analogues to Theorem 2.5 and Proposition 2.9 can be straightforwardly derived, see [START_REF] Amstutz | Topology optimization methods with gradient-free perimeter approximation[END_REF] for full justifications. If A is compact, which happens whenever A is the convolution by a smoothing kernel, then the extension of the function Φ : E → R to Ẽ can be argued by the following Proposition.

Proposition 7.1. If A : L 2 (Ω) → L 2 (Ω) is compact, then Φ is continuous on Ẽ endowed with the weak- * topology of L ∞ (Ω), and it is the relaxation of the function

u ∈ Ẽ → Φ(u) if u ∈ E, +∞ if u / ∈ E.
Proof. The weak continuity is a direct consequence of the compactness of A. To obtain the relaxation result, it suffices (see e.g. Proposition 11.11 of [START_REF] Attouch | Variational analysis in Sobolev and BV spaces[END_REF]) to prove that for any u ∈ Ẽ there exists a sequence (u n ) ∈ E such that u n ⇀ u and Φ(u) ≥ lim sup n→+∞ Φ(u n ). The existence of a sequence of characteristic functions such that u n ⇀ u is obtained by a standard construction, see e.g. Proposition 7.2.14 of [START_REF] Henrot | Variation et optimisation de formes[END_REF], while the second assertion stems from the weak continuity of Φ.

For solving (7.2) at ε fixed, applying directly an alternating algorithm based on the formulation (2.6) would not be well-suited since the minimization with respect to u would be in itself a difficult problem. We could use instead a projected gradient algorithm, however, this would require a line search at every iteration. In fact, an efficient alternating algorithm can still be used thanks to the following lemma. Lemma 7.2. Let X be a Hilbert space and Φ : X → R be a differentiable function such that ∇Φ is λ-Lipschitz. We have for all u ∈ X:

Φ(u) = inf û∈X Φ(û) + ∇Φ(û), u -û + λ 2 u -û 2 .
Proof. On choosing û = u, we immediately get that inf û∈X Φ(û) + ∇Φ(û), u -û + λ 2 u -û 2 ≤ Φ(u).

The reverse inequality stems from Φ(u) -Φ(û) = In our experiments we take the operator A = A q 0 , where A 0 is the discrete convolution operator by the kernel We recall that the unknown u is symmetrized and periodized, hence the convolution is performed without boundary effect. In addition, the computation of the product Au = A q 0 u is efficiently performed through the FFT, without actually computing and storing the matrix of the operator A. We have that A = 1, therefore we choose λ = 2(c 1 -c 2 ) 2 . Assuming that c 1 ≤ f ≤ c 2 , the initialization is done by u = (f -c 1 )/(c 2 -c 1 ). The stopping criterion for the minimization of (7.3) is a threshold on the squared relative variation of u fixed to 10 -6 .

In Figure 15, f is of the form P [0,1] (A(u * c 1 + (1 -u * )c 2 + ν)), with u * a characteristic function and ν a random noise. The grey levels are black and white.

The example featured in Fig. 16 also shows that though the underlying minimal partition algorithm is intended for the classification of plane regions, the proposed method also permits to restore a blurred and noisy text. In order to retrieve the fine details of the text, it is necessary to choose the penalization coefficient α quite small, hence it is difficult to accommodate with a high level of noise. In Fig. 16(b), a better result is obtained when the decrease of ε is stopped before reaching the value ε min , although for this value the obtained image is not binary.

Concluding remarks

Whereas the computer vision community is very active in developing powerful algorithms, few of these approaches are theoretically justified. In this paper we propose a mathematically sound method to perform optimal partitions and apply it to image restoration and classification of greyscale and color images. Let us emphasize that the proposed continuation approach approximates the continuous total variation and not its discretized version. This is achieved through the solving of linear partial differential equations with constant coefficients, which is performed by finite elements and fast Fourier transforms. This task is from far the dominant part of the computational effort of our algorithms. Its

Figure 1 .

 1 Figure 1. Color image classification, from left to right. First: original image (611 × 471 pixels). Second: supervised solution for α = 5.10 -4 √ n (n is the number of pixels). Third and fourth: unsupervised solution for α = 5.10 -4 √ n and α = 10 -2 √ n, respectively. TV-model with H = L 1 , A = I and 5 labels.

Figure 2 .

 2 Figure 2. Greylevel image classification, from left to right. First: original picture (390 × 390 pixels). Second: damaged image with salt and pepper noise. Third: unsupervised solution for α = 4.10 -4 √ n. TV-model with H = L 1 , A = I and 4 levels of grey.

Figure 3 .

 3 Figure 3. Left: original black and white noisy image (320 × 320 pixels). From left to right: obtained solutions for α = 2.10 -4 √ n (incomplete denoising), α = 2.10 -3 √ n (optimal perimeter weight), α = 9.10 -3 √ n (excessive perimeter weight), respectively.

Figure 4 . 2 . 4 . 3 .

 4243 Figure 4. From left to right: original image (320×480 pixels), initialization (ε = ∞), and obtained solutions for ε = ε max /16 and ε = ε min , respectively (with α = 9.10 -3 √ n).

Figure 5 .

 5 Figure 5. 5(a) Triple point: given partition (left), obtained solution with our algorithm (middle), and solution from [17] (right); 5(b) and 5(c) Two quadruple points: given partition (left), obtained solution with our algorithm (middle), and solution from [17] (right).

Figure 6 .

 6 Figure 6. Greylevel image denoising with L 2 norm: original image (left, 100 × 100 pixels), solution with 3 labels and α = 10 -3 √ n (middle), solution taken from [28] (right).

Figure 7 .

 7 Figure 7. Greylevel image denoising with L 2 norm: original image (left, 240 × 240 pixels), solution with 5 labels and α = 10 -3 √ n (middle), solution taken from [28] (right).

Figure 8 .

 8 Figure 8. Greylevel image denoising with L 1 norm: original image (left), noisy image (middle), obtained solution with 3 labels and α = 5.10 -3 √ n (right).

Figure 9 .

 9 Figure 9. From left to right: original image (320×480 pixels) and obtained solutions with 2 and 3 labels, respectively (L 1 norm, α = 5.10 -3 √ n).

Figure 11 .

 11 Figure 11. Unsupervised color image classification with 2 labels and L 1 norm, α = 5.10 -3 √ n: original image (left, 327 × 500 pixels), obtained image (right).

Figure 12 .

 12 Figure 12. Unsupervised color image classification with L 1 norm, α = 5.10 -3 √ n: original image (left, 231×221 pixels), obtained image with 2 labels (middle), obtained image with 5 labels (right).

Figure 13 .

 13 Figure 13. Unsupervised color image classification with L 1 norm and 5 labels, from left to right: original image (320 × 584 pixels) and obtained images for α = 5.10 -3 √ n, α = 10 -2 √ n, α = 2.10 -2 √ n.

Figure 14 .

 14 Figure 14. (a) Anisotropy-based image classification with 2 labels: original image (left, 951 × 1000 pixels), obtained partition for α = 10 -5 √ n (right); (b) Anisotropybased image classification with 4 labels: original image (left, 512 × 512 pixels), obtained partition for α = 10 -6 √ n (right).

1 0 1 0 2 u -û 2 , 2 u -û 2 L 2 + α ε ∇v 2 L 2 (Ω) + 1 ε v 2 L 2 ( 1 -

 112222222221 ∇Φ(û + t(u -û)), u -û dt = ∇Φ(û), u -û + ∇Φ(û + t(u -û)) -∇Φ(û), u -û dt ≤ ∇Φ(û), u -û + λ and the proof is complete.From Lemma 7.2 and (2.6), (7.2) becomes:I ε = min u∈ Ẽ inf û∈L 2 (Ω) inf v∈H 1 (Ω) Φ(û) + ∇Φ(û), u -û + λ (Ω) (Ω) + u, 1 -2v .(7.3)We use an alternating minimization algorithm with respect to the three variables u, û, v. The minimization with respect to v amounts to solving (2.21). The minimization with respect to û is simply achieved by setting ûk = u k-1 . The minimization with respect to u is spatially uncoupled. A short calculation givesu k (x) = P [0,1] ûk (x) -1 λ ∇Φ(û k )(x) + α ε 2v k ) , with P [0,1] (ϕ) = max(0, min(1, ϕ)). Here we have ∇Φ(u) = 2(c 1 -c 2 )A * [(c 1 -c 2 )Au + c 2 -f ],with A * the adjoint operator of A, hence we can choose any λ ≥ 2(c 1 -c 2 ) 2 A * A .

Figure 15 .

 15 Figure 15. Deblurring and denoising: original image (left, 257 × 257 pixels), damaged image with blur and noise effects (middle), reconstructed image for α = 2.10 -4 √ n (right).

Figure 16 .

 16 Figure 16. Deblurring of the text "w ∈ Ξ". 16(a) Original image (left, 245 × 600 pixels), blurred image (middle), reconstructed image for α = 2.10 -5 √ n (right). 16(b) blurred image with noise (left), its restoration with α = 2.10 -5 √ n for ε = 1.99 (middle) and for ε = ε min (right).

Acknowledgements. The work of the third author has been supported in part by Fundação para a Ciência e a Tecnologia (Ciência 2007) and in part by the ERC Advanced Grant "Quasistatic and Dynamic Evolution Problems in Plasticity and Fracture" (grant agreement no. 290888).

n: original image (left, 600 × 600 pixels), obtained images with 2 labels (middle) and with 3 labels (right).

Multilabel classification of color images

The original image f is represented by the three channels (f 3 representing the intensity of red, green and blue, respectively. Each phase Ω i is associated to a color (c i1 , c i2 , c i3 ) ∈ [0, 1] 3 in the same RGB system. The reconstructed image w = (w 1 , w 2 , w 3 ) is given by

where u i is the characteristic function of Ω i . We have for each channel

The difference between the segmented and original images is measured by

We then apply the same algorithm as in Section 4. Note that the geometrical variable u = (u 1 , ..., u N ) as well as the auxiliary variable v = (v 1 , ..., v N ) remain N -dimensional vectors, and the update of levels is separable in the channels.

The first example (see Figure 11) is a two-label problem, initialized with a pure black phase and a pure white phase.

In the second example (see Figure 12), we first choose 2 labels, then 5 labels. In this latter case, the phases are initialized by the pure colors black, red, green, blue and white. We observe that only 4 labels remain at convergence, which is an effect of the perimeter penalization.

In the third example (see Figure 13), we choose again 5 labels. Results obtained with different values of α are depicted.

Anisotropy-based image classification

We come back to a greyscale image f ∈ L ∞ (Ω, [0, 1]). We are given N vectors ξ 1 , ..., ξ N ∈ S 2 , where S 2 is the unit sphere of R 2 . In order to detect fluctuations oriented along the direction ξ i , we implementation could certainly be further improved, but code optimization was not the aim of this work.

We consider the present paper as a first step of an ongoing work where an image partition algorithm and a gradient-free approximation of the perimeter is proposed to address piecewise constant image restoration and/or classification. In particular, we have applied our algorithm to denoising, deblurring and supervised texture identification problems, as well as for inpainting. Future extensions of the method could be multilabel deblurring and unsupervised texture identification.

Appendix A. Weighted median Let (x 1 , ..., x n ) ∈ R n , (α 1 , ..., α n ) ∈ R n be given. We assume that the x i 's are numbered in increasing order. We want to minimize

This function is clearly convex and affine on each interval [x i , x i+1 ]. Therefore, the minimizing set is an interval of the form [x l , x r ], 1 ≤ l ≤ r ≤ n. The subdifferential of V at the point x j is

Then x j is a minimizer of V if and only if 0 ∈ ∂V (x j ), i.e.

-α j ≤ i<j α i -i>j α i ≤ α j .

(A.1)

We obtain x l and x r by checking for each x j if it satisfies (A.1).


Yaozhong Hu
email: hu@math.ku.edu

David Nualart
email: nualart@math.ku.edu

Samy Tindel
email: samy.tindel@univ-lorraine.fr

Fangjun Xu
email: fangjunxu@gmail.com

DENSITY CONVERGENCE IN THE BREUER-MAJOR THEOREM FOR GAUSSIAN STATIONARY SEQUENCES

Consider a Gaussian stationary sequence with unit variance X = {X k ; k ∈ N ∪ {0}}. Assume that the central limit theorem holds for a weighted sum of the form

, where f designates a finite sum of Hermite polynomials. Then we prove that the uniform convergence of the density of V n towards the standard Gaussian density also holds true, under a mild additional assumption involving the causal representation of X.

Introduction

Let X = {X k ; k ∈ N∪{0}} be a centered Gaussian stationary sequence with unit variance. For all v ∈ Z, we set ρ(v) = E[X 0 X |v| ]. Therefore ρ(0) = 1 and |ρ(v)| ≤ 1 for all v. Let γ be the standard Gaussian probability measure and f ∈ L 2 (γ) be a fixed deterministic function such that E[f (X 1 )] = 0. We expand f in the orthonormal basis of Hermite polynomials {H k ; k ≥ 0}, which are more specifically defined in Section 2.2. In particular, if f has Hermite rank d ≥ 1, it admits the following Hermite expansion:

f (x) = ∞ j=d a j H j (x), with a d = 0. Define V n = 1 √ n n-1 k=0 f (X k ).
Then the celebrated Breuer-Major Theorem (see [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF] or Theorem 7.2.4 in [START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF]) can be written as follows:

Theorem 1.1. Suppose that v∈Z |ρ(v)| d < ∞ and set σ 2 = ∞ j=d j!a 2 j v∈Z ρ(v) j , which is assumed to be in (0, ∞). Then the convergence:

V n Law --→ N (0, σ 2 ) (1)
holds true as n tends to infinity.

We shall be in fact interested in a particular case of Theorem 1.1 for finite linear combinations of Hermite polynomials, which is stated here for convenience: Corollary 1.2. Consider 2 ≤ d ≤ q < ∞ and a family of real numbers {a j ; j = d, . . . , q}. Let H j be the jth order Hermite polynomial, and assume that σ 2 ∈ (0, ∞), where σ 2 ≡ q j=d j!a 2 j v∈Z ρ(v) j . Set

V d,q n = 1 √ n n-1 k=0 q j=d a j H j (X k ). (2)
Then V d,q n Law --→ N (0, σ 2 ) as n tends to infinity. In particular, we have:

lim n→0 E V d,q n 4 = 3σ 4 . (3)
Remark 1.3. The relation between Gaussian convergence in law for sequences in a fixed Wiener chaos and behavior of the fourth moment has been extensively studied since the seminal paper [START_REF] Nualart | Central limit theorems for sequences of multiple stochastic integrals[END_REF]. We will need only a small part of the information available on the topic, such as relation [START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF].

Due to its importance, Breuer-Major theorem has been extended and refined in several directions. Important generalizations can be found in Arcones [START_REF] Arcones | Limit theorems for nonlinear functionals of a stationary Gaussian sequence of vectors[END_REF] (multidimensional case), Chambers and Slud [START_REF] Chambers | Central limit theorems for non-linear functionals of stationary Gaussian processes[END_REF] and Giraitis and Surgailis [START_REF] Giriatis | CLT and other limit theorems for functionals of Gaussian processes[END_REF]. A proof of Theorem 1.1 using a combination of Stein's method with Malliavin calculus was given by Nourdin, Peccati and Podolskij in [START_REF] Nourdin | Quantitative Breuer-Major theorems[END_REF], where one can find explicit bounds in the total variation and Wasserstein distances. We refer the reader to the monograph by Nourdin and Peccati [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF] for a more detailed account on this topic.

We shall mainly be concerned here by convergences of densities, and here again the relationship between fourth moment behavior and various type of convergences of random variables in a fixed Wiener chaos have been thoroughly studied in the recent past. The interested reader is referred to [START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF] for further details, but we will use here the following recent criterion (see [START_REF] Nourdin | Fisher Information and the Fourth Moment Theorem[END_REF]Corollary 1.2] and [START_REF] Hu | Convergence of densities of some functionals of Gaussian processes[END_REF]Corollary 4.6]).

Theorem 1.4. Let {F n ; n ∈ N} be a sequence of random variables belonging to a fixed chaos H q with q ≥ 2. Suppose E[F 2 n ] = 1 and lim n→∞ E[F 4 n ] = 3. Let p Fn be the density of the random variable F n and let φ(x) = (2π) -1/2 exp(-|x| 2 /2) be the density of the standard Gaussian distribution on R.

(i) Suppose that for some ǫ > 0,

sup n E DF n -4-ε < ∞.
Then, there exists a constant c such that for all n ≥ 1,

sup x∈R |p Fn (x) -φ(x)| ≤ c E[F 4 n ] -3.
(ii) Suppose that for all p ≥ 1,

sup n E DF n -p < ∞.
Then, for any m ≥ 0, there exists a constant c m such that for all n ≥ 1,

sup x∈R |p (m) Fn (x) -φ (m) (x)| ≤ c m E[F 4 n ] -3.
The goal of the current paper is to apply the criterion given by Theorem 1.4 in order to get convergence of density in the landmark of Breuer-Major theorem. In order to do this we need a uniform estimate on the negative moments of the Malliavin derivative of the sequence, and this is the contents of our main result.

Theorem 1.5. Let X be a Gaussian stationary sequence whose spectral density f ρ satisfies

f ρ ∈ L 1/2 ([-π, π]) and log(f ρ ) ∈ L 1 ([-π, π]) (see the hypothesis 2.
1 and examples in the next section). Let V d,q n be the random variable defined by (2), and assume the hypothesis of Corollary 1.2 to be satisfied. Then for any p ≥ 1, there exists n 0 such that

sup n≥n 0 E DV d,q n -p < ∞. (4)
In the case of a fixed Wiener chaos we can obtain the following consequence.

Corollary 1.6. Under the conditions of Theorem 1.5, if q = d, and we define

F n = V d,d n /σ n , where σ 2 n = E[(V d,d n ) 2 ],
then, for all m ≥ 0 there exists an n 0 (depending on m) such that sup

n≥n 0 sup x∈R |p (m) Fn (x) -φ (m) (x)| ≤ c m E[F 4 n ] -3.
In the case q = d, Theorem 1.4 cannot be applied. In the reference [START_REF] Hu | Convergence of densities of some functionals of Gaussian processes[END_REF] one can find results on the uniform convergence of density for general random variables similar to those stated in Theorem 1.4, but they require a uniform lower bound for the negative moments of the random variable | DF n , -DL -1 F n H |, and our approach does not seem to work in this case because it is not clear how to express DF n , -DL -1 F n H as a sum of squares. Nevertheless, condition (4) allows us to derive the uniform convergence of the densities and their derivatives from a general result proved below (see Proposition 2.6) although in this case we have no information about the rate of convergence. Corollary 1.7. Under the conditions of Theorem 1.5, if we define F n = V d,q n /σ n , where

σ 2 n = E[(V d,q n ) 2 ],
then, for all m ≥ 0 we have lim

n→∞ sup x∈R |p (m) Fn (x) -φ (m) x)| = 0.
Notice that a particular case of Theorem 1.5 has been established in [START_REF] Nourdin | Fisher Information and the Fourth Moment Theorem[END_REF], for q = 2 and X k = B k+1 -B k for a fractional Brownian motion B with Hurst parameter H ∈ (0, 1). The proof of the existence of negative moments for DF n there is based on the Volterra representation of B, which leads to long computations. In comparison our current Theorem 1.5 is more general, since it is valid for a wide class of Gaussian stationary sequences. Its proof is also significantly simplified. These are achieved by the introduction of two new ingredients in the proof, namely:

• A general formula to compute conditional expectations for random variables of the form H q (X k ). • Related to the previous item, we heavily resort to the causal representation of X k , which is particularly convenient in order to compute conditional expectations.

Here is how our paper is structured: we give some preliminary results concerning Gaussian stationary sequences and related Malliavin calculus in Section 2. We then prove our main Theorem 1.5 in Section 3.

Preliminaries

This section is devoted to some preliminaries on causal or moving average representations for Gaussian stationary sequences, as well as Malliavin calculus tools which will be used in the sequel.

2.1. Moving average representation. The classical results on time series presented here are borrowed from [START_REF] Beran | Statistics for long-memory processes[END_REF][START_REF] Brockwell | Time series: theory and methods[END_REF], to which we refer for further details. Start from our Gaussian stationary sequence {X k ; k ∈ N ∪ {0}} with covariance function ρ. We will work under the following assumptions: Hypothesis 2.1. We suppose that ρ admits a spectral density called f ρ , defined by

f ρ (λ) = 1 2π k∈Z ρ(k) e ıkλ , λ ∈ [-π, π]. (5)
We further assume that

f ρ ∈ L 1/2 ([-π, π]) and log(f ρ ) ∈ L 1 ([-π, π]). (6)
Notice that the condition log(f

ρ ) ∈ L 1 ([-π, π]
) is referred to as purely nondeterministic property in the literature.

The interest of dealing with purely nondeterministic sequences is that they admit a socalled causal representation which is particularly convenient for conditional expectation computations. Let us state a result in this direction, which is taken from [4, Chapter 5]. Proposition 2.2. Let X be a Gaussian stationary sequence satisfying Hypothesis 2.1. Then for each k ∈ N ∪ {0} the random variable X k can be decomposed as

X k = j≥0 ψ j w k-j , (7)
where (w k ) k∈Z is a discrete Gaussian white noise and the coefficients ψ j are deterministic. With a slight abuse of notation, extend the sequence ψ to (ψ j ) j∈Z by setting ψ -j = 0 for j ≥ 0. Then one can choose ψ such that it enjoys the following properties:

(i) The sequence ψ admits a spectral density f ψ such that

f ψ = f 1/2 ρ 2π . (ii) In particular, ψ 0 = 1 2π π -π f 1/2 ρ (λ) dλ and ψ 0 > 0. (iii) For all k 1 , k 2 ∈ N we have ρ(k 1 -k 2 ) = k 1 ∧k 2 l=-∞ ψ k 1 -l ψ k 2 -l . Proof.
As mentioned above, those results are classic and borrowed from [START_REF] Brockwell | Time series: theory and methods[END_REF], see in particular formula (5.7.9) therein. Let us justify briefly the fact that ψ 0 > 0. Indeed, the purely nondeterministic assumption imposes log(f ρ ) ∈ L 1 ([-π, π]), and thus f ρ > 0 almost everywhere. This easily yields

ψ 0 = π -π f 1/2 ρ (λ) dλ > 0. Item (iii) is obtained by writing E[X k 1 X k 2 ]
with the expression of X k given by [START_REF] Giriatis | CLT and other limit theorems for functionals of Gaussian processes[END_REF].

Let us now turn to examples for which our standing assumptions of Hypothesis 2.1 are met. The following proposition provides two typical and classical cases for which a spectral density exists and satisfies some integrability properties. Proposition 2.3. Let ρ be the covariance function of X. We have the following statements. (i) If ρ ∈ ℓ 1 , then the spectral density f ρ exists and is a nonnegative L 2 function defined on [-π, π]. Condition (5) is thus fulfilled. (ii) If lim k→∞ |k| α ρ(k) = c ρ for some α ∈ (0, 1) and some positive constant c ρ , then the spectral density exists, is strictly positive almost everywhere and satisfies lim λ→0 |λ| 1-α f ρ (λ) = c f . In particular, condition (5) is satisfied.

We now give two specific and important examples which satisfy Hypothesis 2.1.

Example 2.4. The so-called autoregressive fractionally integrated moving-average (ARFIMA) processes are introduced in [START_REF] Granger | An introduction to long-memory time series models and fractional differencing[END_REF] and [START_REF] Hosking | Fractional differencing[END_REF]. Denote by B the one lag backward operator (BX k = X k-1 ). Let φ(z) and θ(z) be two polynomials which have no common zeros and such that the zeros of φ lie outside the closed unit disk {z , |z| ≤ 1}. Suppose that X k is given by

φ(B)X k = (Id -B) -d θ(B)w k , (8)
where -1 < d < 1/2, and where the operator (Id -B) -d is defined by:

(Id -B) -d = ∞ j=1 η j B j with η j = Γ(d + j) Γ(j + 1)Γ(d)
.

Also notice that in (8) the sequence (w k ) k∈Z is a discrete Gaussian white noise. It is wellknown (see [START_REF] Palma | Long-memory time series[END_REF], Theorem 3.4 and equation (3.19)) that under the above conditions, {X k , k ∈ N} admits a spectral density whose exact expression is:

f (λ) = 1 2π 2 sin λ 2 -2d |θ(e -iλ )| 2 |φ(e -iλ )| 2 .
It is thus readily checked that the conditions ( 5) and ( 6) are satisfied, and hence X k has a causal representation.

Example 2.5. Our second example is the fractional Gaussian noise. Let {B t , t ≥ 0} be a fractional Brownian motion of Hurst parameter H ∈ (0, 1). Then

{X k = B k+1 -B k , k ∈ N ∪ {0}} is a stationary Gaussian process with correlation ρ(k) = 1 2 (k + 1) 2H -2k 2H -(k -1) 2H .
Its spectral density (see e.g. [START_REF] Beran | Statistics for long-memory processes[END_REF], equation (2.17)) is:

f (λ) = 1 2π ∞ k=-∞ ρ(|k|)e iλ = 2c f (1 -cos(λ)) ∞ j=-∞ |2πj + λ| -2H-1 , λ ∈ [-π, π] ,
where

c f = (2π) -1 sin(πH)Γ(2H + 1). If H ≤ 1/2, it is clear that ∞ k=-∞ |ρ(|k|)| < ∞. This implies sup λ∈[-π,π] |f (λ)| < ∞ . Thus f ∈ L 1/2 . If 1/2 < H < 1, then 0 ≤ f (λ) ≤ 2c f (1 -cos(λ))|λ| -2H-1 + 2c f j =0 |2πj + λ| -2H-1 , λ ∈ [-π, π] .
The first term is in

L 1 since H < 1. When j = 0, π -π |2πj + λ| -2H-1 dλ ≤ Cj -2H for some positive constant C. Thus π -π j =0 |2πj + λ| -2H-1 dλ < ∞,
owing to the fact that H > 1/2. Therefore, we have f ∈ L 1 and hence f ∈ L 1/2 . Summarizing we have f ∈ L 1/2 for all H ∈ (0, 1). This also implies log

+ f (λ) ∈ L 1 . To see log -f (λ) ∈ L 1 , we notice that f (λ) ≥ 2c f (1 -cos(λ))|λ| -2H-1 . So log -f (λ) ≤ C + log (1 -cos(λ))|λ| -2H-1 which is in L 1 .
In conclusion, the sequence X satisfies Hypothesis 2.1.

Malliavin calculus.

We start by briefly recalling some basic notation and results connected to Gaussian analysis and Malliavin calculus. The reader is referred to [START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF][START_REF] Nualart | The Malliavin calculus and related topics[END_REF] for details or missing proofs.

2.2.1.

Wiener space and generalities. Let H be a real separable Hilbert space with inner product •, • H . The norm of H will be denoted by • = • H . Recall that we call isonormal Gaussian process over H any centered Gaussian family W = {W (h) : h ∈ H}, defined on a probability space (Ω, F , P) and such that E[W (h)W (g)] = h, g H for every h, g ∈ H.

In our application the underlying Gaussian family will be a discrete Gaussian white noise (w k ) k∈Z . The space H is given here by H = ℓ 2 (Z) (the space of square integrable sequences indexed by Z) equipped with its natural inner product. Set {ε j ; j ∈ Z} for the canonical basis of ℓ 2 (Z), that is ε j k = δ j (k). We thus identify w j with W (ε j ). Assume from now on that our underlying σ-algebra F is generated by W .

For any integer q ∈ N ∪ {0}, we denote by H q the qth Wiener chaos of W . We recall that H 0 is simply R whereas, for any q ≥ 1, H q is the closed linear subspace of L 2 (Ω) generated by the family of random variables {H q (W (h)), h ∈ H, h H = 1}, with H q the q-th Hermite polynomial given by

H q (x) = (-1) q e x 2 2 d q dx q e -x 2 2 . ( 9
)
Let S be the set of all cylindrical random variables of the form

F = g(W (h 1 ), . . . , W (h n )),
where n ≥ 1, h i ∈ H, and g is infinitely differentiable such that all its partial derivatives have polynomial growth. The Malliavin derivative of F is the element of L 2 (Ω; H) defined by

DF = n i=1 ∂g ∂x i (W (h 1 ), . . . , W (h n )) h i .
By iteration, for every m ≥ 2, we define the mth derivative D m F . This is an element of L 2 (Ω; H ⊙m ), where H ⊙m designates the symmetric mth tensor product of H. For m ≥ 1 and p ≥ 1, D m,p denote the closure of S with respect to the norm • m,p defined by

F p m,p = E[|F | p ] + m j=1 E D j F p H ⊗j . Set D ∞ = ∩ m,p D m,
ϕ : R n → R is in C 1 b
(that is, belongs to the set of continuously differentiable functions with a bounded derivative) and if {F i } i=1,...,n is a vector of elements of D 1,2 , then ϕ(F 1 , . . . , F n ) ∈ D 1,2 and

Dϕ(F 1 , . . . , F n ) = n i=1 ∂ϕ ∂x i (F 1 , . . . , F n ) DF i . (10
E[ DF n -p ] < ∞. (11)
Then, for all m ≥ 0 the derivative p (m) n of the density of F n , converges uniformly and in in L p (R) for all p ≥ 1 to the corresponding derivative of the density of F ∞ .

Proof. First notice that by condition [START_REF] Nourdin | Fisher Information and the Fourth Moment Theorem[END_REF], the random variable F n has an infinitely differentiable density p n , whose derivatives can be expressed by

p (m) n (x) = E 1 {F >x} G (m) n , (12)
where the random variables G n ] = 0, we get:

sup n p (m) n (x) ≤ sup n P(|F n | > |x|) E |G (m) n | 2 1 2 ≤ c m,2 sup n P(|F n | > |x|) ≤ c m,2 sup n E[|F n | q ] 1 2 |x| -q 2 ≤ c|x| -q 2 ,
for all q ≥ 1 and for some constant c depending on q and m. By the results of [START_REF] Nourdin | Convergence in total variation on Wiener chaos[END_REF], the densities p n converge in L 1 (R) to the density p ∞ of F ∞ . The boundedness in L p (R) and the uniform bound of p n imply that this convergence holds in L p (R) for any p ≥ 1. From a compactness argument we can deduce the convergence of all the derivatives in L p (R) for any p ≥ 1. Finally, the uniform convergence is also easy to establish.

2.4.

A key lemma. Our future computations will heavily rely on an efficient way to compute conditional expectations. Towards this aim, we state here some general results. Let us start with a decomposition for Hermite polynomials: Lemma 2.7. For any q ≥ 1, let H q be the polynomial defined by relation [START_REF] Hosking | Fractional differencing[END_REF]. Consider y, z ∈ R and two real parameters a, b with a 2 + b 2 = 1. Then the following relation holds true:

H q (ay + bz) = q ℓ=0 q ℓ a q-ℓ b ℓ H q-ℓ (y) H ℓ (z). ( 13
)
Proof. By the definition of the Hermite polynomials, we have

e aty-(at) 2 2 = ∞ i=0 (at) i H i (y), and e tbz-(bt) 2 2 = ∞ j=0 (bt) j H j (z). (14)
In the same way, we also obtain:

e t(ay+bz)-t 2 /2 = ∞ q=0
t q H q (ay + bz).

Since a 2 + b 2 = 1, we obviously have e aty-(at) 2 2 e tbz-(bt) 2 2 = e t(ay+bz)-t 2 /2 . Thus, multiplying the right hand sides of ( 14) we recover the right hand side of (15), namely:

∞ q=0 t q H q (ay + bz) = ∞ i=0 (at) i H i (y) ∞ j=0 (bt) j H j (z),
which easily yields the desired identity [START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF].

With this preliminary result in hand, we are ready to state our result on conditional expectations: Proposition 2.8. Let Y and Z be two centered Gaussian random variables such that Y is measurable with respect to a σ-algebra G ⊂ F and

Z is independent of G. Assume that E[Y 2 ] = E[Z 2 ] = 1.
Then for any q ≥ 1, and real parameters a, b such that a 2 + b 2 = 1, we have:

E[H q (aY + bZ)|G] = a q H q (Y ). (16)
Proof. Apply identity [START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF] in order to decompose H q (aY + bZ). Then identity [START_REF] Nourdin | Convergence in total variation on Wiener chaos[END_REF] follows easily from the fact that Y is G-measurable, Z is independent from G and Hermite polynomials have 0 mean under a centered Gaussian measure except for H 0 ≡ 1.

2.5. Carbery-Wright inequality. In the proof of Theorem 1.5, we will make use of the following inequality due to Carbery and Wright [5, Theorem 8], which is recalled here for convenience:

Proposition 2.9. Let X = (X 1 , . . . , X n ) be a Gaussian random vector in R n and Q : R n → R a polynomial of degree at most m. Then there is a universal constant c > 0 such that:

E[|Q(X 1 , . . . , X n )|] 1 m P(|Q(X 1 , . . . , X n )| ≤ x) ≤ c m x 1 m , for all x > 0. (17)
3. Proof of Theorem 1.5

In this section, we will prove our main result, which amounts to show the inequality (4). This will be done into several steps.

Step 1: Computation of the Malliavin norm. Recall that V d,q n is defined by relation [START_REF] Beran | Statistics for long-memory processes[END_REF], and set for the moment f = q j=d a j H j . Invoking relation [START_REF] Hu | Convergence of densities of some functionals of Gaussian processes[END_REF], plus the fact that Dw j = ε j with the notation of Section 2.2.1, we get:

DV d,q n = 1 √ n n-1 k=0 f ′ (X k ) j≥0 ψ j ε k-j = 1 √ n n-1 l=-∞ n-1 k=l + ψ k-l f ′ (X k ) ε l , (18)
where l + = max{l, 0}. Invoking Proposition 2.3, item (iii), it is thus readily checked that:

DV d,q n , DV d,q n H = 1 n n-1 k 1 ,k 2 =0 f ′ (X k 1 )ρ(k 1 -k 2 )f ′ (X k 2 ),
where we recall that ρ is the covariance function of the Gaussian stationary sequence {X k ; k ≥ 0}. This is consistent with the expression found in [START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF]Chapter 5]. However, in order to write the above expression as sum of some squares, we will start directly from expression [START_REF] Nualart | Central limit theorems for sequences of multiple stochastic integrals[END_REF]. Since {ε l ; l ∈ Z} is an orthonormal basis of ℓ 2 (Z) we obtain:

DV d,q n , DV d,q n H = 1 n n-1 ℓ=-∞ n-1 k=ℓ + ψ k-ℓ f ′ (X k ) 2 .
Rearranging terms (namely, change kℓ to k and then nℓ -1 to m), we end up with:

DV d,q n , DV d,q n H ≥ 1 n n-1 ℓ=0 n-ℓ-1 k=0 f ′ (X ℓ+k ) ψ k 2 = 1 n n-1 m=0 m k=0 f ′ (X n-1-(m-k) ) ψ k 2 ≡ A n .
As a last preliminary step we resort to the fact that X = {X k ; k ∈ N ∪ {0}} is a Gaussian stationary sequence, which allows to assert that A n is identical in law to B n with

B n := 1 n n-1 m=0 m k=0 f ′ (X m-k ) ψ k 2 = 1 n n-1 m=0 m k=0 f ′ (X k ) ψ m-k 2 .
We will now bound the negative moments of B n .

Step 2: Block decomposition. We now wish to apply the Carbery-Wright inequality [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] in order to get bounds for negative moments of B n . However, relation [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] only applies to moments of small order, and this is why we proceed to a decomposition of B n into smaller blocks.

Fix thus an integer N ≥ 1 and let M = [n/N] be the integer part of n/N. Then n ≥ NM and as a consequence,

B n = 1 n n-1 m=0 m k=0 f ′ (X k ) ψ m-k 2 ≥ 1 n N -1 i=0 (i+1)M -1 m=iM m k=0 f ′ (X k ) ψ m-k 2
.

For i = 0, . . . , N -1, define

B i n = 1 n (i+1)M -1 m=iM m k=0 f ′ (X k ) ψ m-k 2 so that B n ≥ N -1 i=0 B i n .
Then it is readily checked that:

(B n ) -p 2 ≤ N -1 i=0 (B i n ) -p 2N . ( 19
)
Recall once again the representation of the sequence X in [START_REF] Giriatis | CLT and other limit theorems for functionals of Gaussian processes[END_REF], and denote by F k the filtration generated by {w ℓ : ℓ < k}. Then starting from [START_REF] Palma | Long-memory time series[END_REF] we obtain:

E (B n ) -p 2 ≤ E N -1 i=0 (B i n ) -p 2N = E E (B N -1 n ) -p 2N |F (N -1)M N -2 i=0 (B i n ) -p 2N . (20)
Step 3: Application of Carbery-Wright. Let us go back to the particular situation of f = q j=d a j H j , which means in particular that f ′ = q j=d j a j H j-1 . We are now in a position to apply a conditional version of inequality [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] to the block (B N -1 n ) -p 2N in (20). First, we notice

E (B N -1 n ) -p 2N |F (N -1)M ≤ 1 + p 2N 1 0 P B N -1 n ≤ x| F (N -1)M x -p 2N -1 dx . ( 21
) Since B N -1
n is a polynomial of order m = 2(q -1), Carbery-Wright's inequality (17) yields:

P B N -1 n ≤ x|F (N -1)M ≤ c x 1 2(q-1) E B N -1 n |F (N -1)M 1 2(q-1) . ( 22
)
Step 4: Estimates for the conditional expectation. We now estimate the conditional expectation

E[B N -1 n |F (N -1)M ].
We have:

E B N -1 n |F (N -1)M = 1 n N M -1 m=(N -1)M E   m k=0 f ′ (X k ) ψ m-k 2 F (N -1)M   ≥ 1 n N M -1 m=(N -1)M A m , (23)
where we have set

A m = Var m k=(N -1)M f ′ (X k ) ψ m-k F (N -1)M .
Furthermore, notice that

f ′ (X k ) = f ′ k ℓ=-∞ ψ k-i w i = f ′ (Y k + Z k ),
where

Y k = (N -1)M -1 i=-∞ ψ k-i w i is F (N -1)M -measurable and Z k = k i=(N -1)M ψ k-i w i is inde- pendent of F (N -1)M . Recalling that f ′ = q
j=d j a j H j-1 , we can thus resort to Lemmas 2.7 and 2.8. This gives:

H q-1 (X k ) -E[H q-1 (X k )|F (N -1)M ] = q j=d j-1 ℓ=1 j 2 a j j -1 ℓ σ j-1-ℓ Y k H j-1-ℓ ( Y k )σ ℓ Z k H ℓ ( Z k ), where σ Y k = [Var(Y k )] 1/2 , σ Z k = [Var(Z k )] 1/2 , Y k = Y k /σ Y k and Z k = Z k /σ Z k . Therefore, A m = E     m k=(N -1)M q j=d j-1 ℓ=1 a j,ℓ,k H j-1-ℓ ( Y k ) H ℓ ( Z k ) ψ m-k   2 F (N -1)M   = E     q-1 ℓ=1 m k=(N -1)M q j=(ℓ+1)∨d a j,ℓ,k H j-1-ℓ ( Y k ) H ℓ ( Z k ) ψ m-k   2 F (N -1)M   ,
where we have set a j,ℓ,k = j 2 a j j-1

ℓ σ j-1-ℓ Y k σ ℓ Z k .
Recall that the random variables Y k are F (N -1)M -measurable while the random variables Z k are independent of F (N -1)M . By decorrelation properties of Hermite polynomials we thus get:

A m = q-1 ℓ=1 E     m k=(N -1)M q j=(ℓ+1)∨d a j,ℓ,k H j-1-ℓ ( Y k ) H ℓ ( Z k ) ψ m-k   2 F (N -1)M  
and we trivially lower bound this quantity by taking the term corresponding to ℓ = q -1. In this situation the sum q j=(ℓ+1)∨d is reduced to the term corresponding to j = q, and since a q,q-1,k = q 2 a q σ q-1 Z k we obtain:

A m ≥ E     m k=(N -1)M q 2 a q σ q-1 Z k H q-1 ( Z k ) ψ m-k   2 F (N -1)M   = q 4 a 2 q E     m k=(N -1)M σ q-1 Z k H q-1 ( Z k ) ψ m-k   2   .
We now invoke the identity E[H p ( Zk 1 )H p ( Zk 2 )] = 1 p! (E[ Zk 1 Zk 2 ]) p in order to obtain

A m ≥ q 5 a 2 q q! m k 1 ,k 2 =(N -1)M σ q-1 Z k 1 σ q-1 Z k 2 E Z k 1 Z k 2 q-1 ψ m-k 1 ψ m-k 2 .
Furthermore, it is readily checked that:

E Z k 1 Z k 2 = 1 σ Z k 1 σ Z k 2 k 1 ∧k 2 i=(N -1)M ψ k 1 -i ψ k 2 -i ,
and thus

A m ≥ q 5 a 2 q q! m k 1 ,k 2 =(N -1)M   k 1 ∧k 2 i=(N -1)M ψ k 1 -i ψ k 2 -i   q-1 ψ m-k 1 ψ m-k 2 = q 5 a 2 q
q! m i 1 ,...,i q-1 =(N -1)M m k 1 ,k 2 =max(i 1 ,...,i q-1 )

ψ m-k 1 ψ m-k 2 q-1 j=1
ψ k 1 -i j ψ k 2 -i j = q 5 a 2 q q! m i 1 ,...,i q-1 =(N -1)M   m k=max(i 1 ,...,i q-1 )

ψ m-k q-1 j=1 ψ k-i j   2 .
Here again, this sum of squares is trivially lower bounded by taking the term corresponding to i 1 = • • • = i q-1 = m, which yields:

A m ≥ c a,q,ψ with c a,q,ψ ≡ q 5 a 2 q q! ψ 2q 0 > 0.

(24)

Step 5: Conclusion. In the remainder of the proof the constants c a,q,ψ,N and so can change from line to line without further mention. Plugging relation (24) into (23) and recalling that N is a given integer whose exact value will be fixed below, we get:

E B N -1
n |F (N -1)M ≥ M c a,q,ψ n ≥ c a,q,ψ,N > 0, as long as N stays bounded. We then insert back this inequality into (21) and ( 22) in order to get:

P B N -1 n ≤ x|F (N -1)M ≤ 1 + p c a,q,ψ,N 2N 1 0
x 1 2(q-1) -p 2N -1 dx = c a,q,ψ,N,p < ∞, where we have chosen N such that p 2N < 1 2(q-1) . Iterating this bound into (20), we have thus obtained:

E (B n ) -

,

  for any m ≥ 1. From this formula it follows that the derivatives p (m) n are uniformly bounded and also uniformly bounded in L p (R) for all p ≥ 1. Indeed, by[START_REF] Nourdin | Convergence in total variation on Wiener chaos[END_REF] Lemma 2.4] we have sup n E[|F n | r ] < ∞ for all r ≥ 1. This uniform bound on the moments, together with the equivalence of the • m,p norms in any finite sum of Wiener chaoses and condition[START_REF] Nourdin | Fisher Information and the Fourth Moment Theorem[END_REF] imply that sup n G (m) n L p (Ω) = c m,p < ∞ for all m ≥ 0. Then, we can write from[START_REF] Nourdin | Stein's method on Wiener chaos[END_REF] supn sup x p (m) n (x) ≤ sup n E |G (m) n | = c m,1 < ∞,and, using the fact that E[G (m)

p 2 ≤

 2 c N a,q,ψ,N,p , which is a finite quantity. Finally recall from Step1 that E[(B n ) -p 2 ] = E[ DV d,q n -pH ], which finishes the proof.

  p . One can then extend the definition of D m to D m,p . When m = 1, one simply write D instead of D 1 . As a consequence of the hypercontractivity property of the Ornstein-Uhlenbeck semigroup (see, e.g.,[START_REF] Nourdin | Normal Approximations using Malliavin Calculus: from Stein's Method to the Universality[END_REF] Theorem 2.7.2]), all the • m,p -norms are equivalent in any finite sum of Wiener chaoses.Finally, let us recall that the Malliavin derivative D satisfies the following chain rule: if

  ) 2.3. Convergence of densities. Suppose that F is a random variable in D ∞ such that E[ DF -p ] < ∞ for all p ≥ 1.Then, we know that F has an infinitely differentiable density and there are explicit formulas for the density and its derivatives (see[START_REF] Nualart | The Malliavin calculus and related topics[END_REF] Proposition 2.1.5]). Using this result, we can establish the following criteria for convergence of densities for random variables in a finite sum of Wiener chaos.

Proposition 2.6. Let {F n ; n ∈ N} be a sequence of random variables belonging to a finite sum of Wiener chaos ⊕ M k=1 H k , which converges in law to a nonzero random variable F ∞ . Suppose that, for all p ≥ 1 sup n

Y. Hu is partially supported by a grant from the Simons Foundation #209206. D. Nualart is supported by the NSF grant DMS1208625. F. Xu is supported by the 111 Project(B14019). This project has been carried out while S. Tindel was on sabbatical at the University of Kansas. He wishes to express his gratitude to this institution for its warm hospitality.