Caroline Meynet

Cathy Maugis-Rabusseau

A sparse variable selection procedure in model-based clustering

Keywords: Variable selection, Gaussian mixture clustering, Non asymptotic penalized criterion, ℓ 1 -regularization method, Maximum likelihood estimator

Owing to the increase of high-dimensional datasets, the variable selection for clustering is an important challenge. In the context of Gaussian mixture clustering, we recast the variable selection problem into a general model selection problem. Our procedure first consists of using a ℓ 1 -regularization method to build a data-driven model subcollection. Second, the maximum loglikelihood estimators (MLEs) are obtained using the EM algorithm. Next a non asymptotic penalized criterion is proposed to select the number of mixture components and the relevant clustering variables simultaneously. A general model selection theorem for MLEs with a random model collection is established. It allows one to derive the penalty shape of the criterion, which depends on the complexity of the random model collection. In practice, the criterion is calibrated using the so-called slope heuristics. The resulting procedure is illustrated on two simulated examples. Finally, an extension to a more general modeling of irrelevant clustering variables is presented.

Une procédure parcimonieuse de sélection de variables pour la classification par mélanges gaussiens

Résumé : Au vu de l'augmentation du nombre de jeux de données de grande dimension, la sélection de variables pour la classification non supervisée est un enjeu important. Dans le cadre de la classification par mélanges gaussiens, nous reformulons le problème de sélection de variables en un problème général de sélection de modèle. Dans un premier temps, notre procédure consiste à construire une sous-collection de modèles grâce à une méthode de régularisation ℓ 1 . Puis, l'estimateur du maximum de vraisemblance est déterminé via un algorithme EM pour chaque modèle. Enfin un critère pénalisé non asymptotique est proposé pour sélectionner à la fois le nombre de composants du mélange et l'ensemble des variables informatives pour la classification. D'un point de vue théorique, un théorème général de sélection de modèles dans le cadre de l'estimation par maximum de vraisemblance avec une collection aléatoire de modèles est établi. Il permet en particulier de justifier la forme de la pénalité de notre critère, forme qui dépend de la complexité de la collection de modèles. En pratique, ce critère est calibré grâce à la méthode dite de l'heuristique de pente. Cette procédure est illustrée sur deux jeux de données simulées. Finalement, une extension, associée à une modélisation plus générale des variables non informatives pour la classification, est proposée.

Mots-clés :

Sélection de variables, Classification par mélanges gaussiens, Critère pénalisé non asymptotique, Méthode de régularisation ℓ 1 , Estimateur du maximum de vraisemblance.

Introduction

The goal of clustering methods is to discover clusters among n individuals described by p variables. Many clustering methods exist and roughly fall into two categories. The first one is distance-based clustering methods, including hierarchical clusterings and K-means type algorithms. The second category is model-based clustering methods: each cluster is represented by a parametric distribution, the entire dataset is modeled by a mixture of these distributions, and a criterion is used to optimize the fit between the data and the model. An advantage of model-based clustering is to provide a rigorous statistical framework to assess the number of clusters and the role of each variable in the clustering process. In this paper, we focus on clustering with Gaussian mixtures.

Variable selection for clustering is an important challenge, motivated by the increasing study of high-dimensional datasets. Since the structure of interest may often be contained into a subset of the available variables and many attributes may be useless or even harmful to detect a reasonable clustering structure, it is important to select the relevant clustering variables. In addition, removing the irrelevant variables enables to get simpler modeling and can largely enhance interpretability. Usually, two types of variable selection approaches are envisaged. On the one hand, the "filter" approaches select the variables before the cluster analysis (see for instance [START_REF] Dash | Feature selection for clustering -a filter solution[END_REF][START_REF] Jouve | A filter feature selection method for clustering[END_REF]. Their main weakness is the independence between the variable selection step and the clustering procedure. In contrast, the "wrapper" approaches combine variable selection and clustering. For distance-based clustering, one can cite the works of [START_REF] Fowlkes | Variable selection in clustering[END_REF], [START_REF] Inria Devaney | Efficient feature selection in conceptual clustering[END_REF] and [START_REF] Brusco | A variable selection heuristic for k-means clustering[END_REF] for instance. Wrapper methods are also developed in the model-based clustering framework. For instance, [START_REF] Maugis | Variable selection in model-based clustering: A general variable role modeling[END_REF], which is an extension of [START_REF] Raftery | Variable Selection for Model-Based Clustering[END_REF] and Maugis et al. (2009a), propose a general variable role modeling (relevant, redundant and irrelevant variables for clustering) and a backward stepwise algorithm is developed. In the same low dimensional context, Maugis and Michel (2011b,a) recast also the variable selection problem for clustering into a general model selection problem. For each model, the maximum likelihood estimator is considered and a data-driven penalized criterion is built for the model selection. Nevertheless, their procedure, called MM-MLE procedure in this paper, requires to consider the complete variable subset collection (or to preliminary order the variables). In the high dimensional context, some Bayesian methods have been developed as [START_REF] Tadesse | Bayesian variable selection in clustering highdimensional data[END_REF] and [START_REF] Kim | Variable selection in clustering via Dirichlet process mixture models[END_REF]. Rather than considering a Bayesian approach, [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] propose to take advantage of the sparsity property of ℓ 1 -penalization to perform automatic variable selection for high-dimensional data clustering. Their procedure, called PS-Lasso procedure in this paper, consists of using a Lasso method to select relevant clustering variables and estimate mixture parameters in the same exercise. Nevertheless, ℓ 1 -penalization induces shrinkage of the coefficients and thus biased estimators with high estimation risk. Moreover, they use a BIC-type criterion for the model selection which can be unsuitable for high-dimensional data.

In this paper, a global model selection procedure, called Lasso-MLE procedure, is proposed to simultaneously choose the number of clusters and the relevant clustering variables. This procedure is aiming to take advantage of MM-MLE and PS-Lasso procedures in low-and high-dimensional cases. Fol-lowing [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF], an ℓ 1 -penalized likelihood approach is considered to determine potential sets of relevant variables. This allows one to efficiently construct a data-driven model subcollection with reasonable complexity, even for high-dimensional situations. Contrary to Pan and Shen's approach, the evaluation of the maximum likelihood estimator (MLE) rather than the ℓ 1 -penalized maximum likelihood estimator for each model is considered to avoid estimation problems due to ℓ 1 -penalization shrinkage. Next, a non asymptotic penalized criterion is proposed to solve the model selection problem. Considering MLE and a random model collection require to extend the general model selection theorem of Massart (2007, Theorem 7.11) to our context. This extension allows one to justify the penalty shape of our criterion. In practice, the penalty depending on unknown constant(s) is calibrated using the so-called slope heuristics [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF][START_REF] Baudry | Slope heuristics: overview and implementation[END_REF]. To our knowledge, the proposed Lasso-MLE estimators have never been studied in model-based clustering. Yet, this idea has emerged in other frameworks. In regression, [START_REF] Connault | Calibration d'algorithmes de type Lasso et analyse statistique de données métallurgiques en aéronautique[END_REF] proposed such a procedure, which seems to have better performances than for the classical Lasso method. Such an estimator is also mentioned as LARS-OLS hybrid in Efron et al. (2004, p. 421). In the density estimation framework, [START_REF] Bertin | Adaptive dantzig density estimation[END_REF] consider such an idea to estimate densities decomposed on some dictionary. Note that in this paper, only finite Gaussian mixtures with common spherical covariance matrix are considered in order to simplify the reading. This assumption allows one to reduce the relevance of variables for the clustering on the mean component vectors. The extension of the Lasso-MLE procedure to general Gaussian mixtures is discussed in Section 7.

The paper is organized as follows. Clustering with Gaussian mixtures is recalled in Section 2.1 and the model collection is described in Section 2.2. Section 3 is devoted to the description of MM-MLE, PS-Lasso and Lasso-MLE procedures. A theoretical result for the Lasso-MLE estimator is stated in Section 4.1 and the penalty calibration for the data-driven penalized criterion is explained in Section 4.2. The three procedures are compared on two simulated examples in Section 5. Section 6 is devoted to the extension of the Lasso-MLE procedure in a curve clustering context. The paper is concluded with a brief discussion in Section 7. Technical aspects are given in Appendix.

2 Gaussian mixture models

Clustering with Gaussian mixtures

Consider a dataset Y = (Y 1 , . . . , Y n ) with Y i ∈ R p .
These data come from a probability distribution with unknown density s. This target s is estimated by a finite mixture model in a clustering purpose although s is not assumed to be a Gaussian mixture density itself. Model-based clustering consists of assuming that the data come from several subpopulations and that the overall population is a mixture of these subpopulations. In this paper, each mixture component is modeled by a Gaussian density. Thus the distribution s is modeled by a finite Gaussian mixture with K components:

y ∈ R p → K k=1 π k Φ(y | µ k , Σ k ).

Inria

The mixing proportions (π 1 , . . . , π K ) belong to

Π K = (a 1 , . . . , a K ) ∈ (0, 1) K ; K k=1 a k = 1 . The density Φ(• | µ k , Σ k )
is the p-dimensional Gaussian density with mean µ k and covariance matrix Σ k . After estimating the parameter vector, the data clustering can be obtained using the Maximum A Posteriori (MAP) principle:

i is assigned to Cluster k if πk Φ Y i | μk , Σk > πℓ Φ Y i | μℓ , Σℓ for all ℓ = k.

The model collection

Currently, statistics deals with problems where data are described by many variables. In principle, the more information one has about each individual, the better a clustering method is expected to perform. Nevertheless, some variables can be useless or even harmful to obtain a good data clustering. Thus, it is important to determine the relevant variables for the Gaussian mixture clustering process.

In this paper, spherical Gaussian mixtures are considered where covariances fulfill Σ k = σ 2 I p for all k = 1, . . . , K with σ 2 > 0. Then the clusters are characterized by the mean parameters and a variable j is irrelevant for the clustering if µ kj is independent of k. Without loss of generality, the data are assumed to have a null expectation: E[Y ij ] = 0. In practice, empirical centering of the data is performed to ensure this assumption. Since the expectation E[Y ij ] is estimated by K k=1 π k µ kj , a variable j is called irrelevant for the clustering if µ kj = 0 for all k = 1, . . . , K, otherwise it is called relevant for the clustering. In the sequel, J r denotes the subset of relevant variables and J c r := {1, . . . , p} \ J r is the irrelevant variable subset. Moreover, for all y ∈ R p , y [Jr ] denotes the restriction of y on J r .

Consequently, the following model collection indexed by the number of clusters K ∈ N * and a relevant variable subset J r , is considered to estimate the density s:

S (K,Jr ) =      y ∈ R p → s θ (y) = Φ(y [J c r ] | 0, σ 2 I) K k=1 π k Φ(y [Jr ] | µ k , σ 2 I) θ = (π 1 , . . . , π K , µ 1 , . . . , µ K , σ) ∈ Π K × R |Jr | K × R * +      .
(1) The dimension of a model S (K,Jr ) corresponds to the total number of free parameters estimated in the model:

D (K,Jr ) = K(1 + |J r |).

Three competitive procedures for variable selection in clustering

In this paper, a variable selection procedure is proposed to select relevant variables for improving data clustering. This procedure is a compromise between MM-MLE and PS-Lasso procedures proposed by Maugis and Michel (2011a,b) and [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] respectively. These three variable selection procedures are presented in this section and summarized in Table 1.

MM-MLE procedure

In their procedure, Maugis and Michel (2011a,b) consider a model collection {S (K,Jr ) } (K,Jr )∈M where S (K,Jr ) is defined by (1) and the model collection is indexed by M = N * × J , where J denotes the collection of the non empty subsets of {1, . . . , p}. For each (K, J r ) ∈ M, the maximum likelihood estimator

ŝ(K,Jr) = argmin t∈S (K,Jr ) γ n (t) with γ n (t) = - 1 n n i=1 ln[t(Y i )]
is computed using an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF]. In practice, Mixmod software [START_REF] Biernacki | Model-based cluster and discriminant analysis with the mixmod software[END_REF] or mclust [START_REF] Fraley | Enhanced software for model-based clustering, density estimation, and discriminant analysis: mclust[END_REF] can be used for instance. Next, they propose to solve the model selection problem with a data-driven penalized criterion. The theoretical construction of this criterion is the topic of [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF]. The resulting penalty is proportional to D/n, up to an unknown multiplicative constant κ. In practice, this penalty is calibrated using the slope heuristics [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF], recalled in Section 4.2. First, models are grouped according to their dimension D in order to obtain a model collection {S D } D∈D . For each dimension D ∈ D, let ŝD be the maximum likelihood estimator in S D and (K D , J D ) such that ŝD = ŝ(K D ,J D ) . Second, the slope κ is estimated to calibrate the penalty, based on the linear behavior of the function D/n → -γ n (ŝ D ) for large dimensions (see [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF], for more details). Third, the minimizer D of the penalized criterion

D = argmin D∈D γ n (ŝ D ) + 2κ D n (2)
is determined and ( K, Ĵr ) := (K D , J D ) is selected. Finally, a data clustering is derived from ŝ( K, Ĵr ) by applying the MAP principle.

The MM-MLE procedure is only appropriate for low dimensional problem (p ≤ n). And, since the model collection is indexed by the complete collection of subsets J , the procedure is time consuming. When p is very large, the variables may be ordered but this preliminary step is limiting in practice. Thus, alternative variable selection procedures have to be considered for clustering high-dimensional data.

PS-Lasso procedure

Unlike Maugis and Michel (2011a), [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] do not consider a deterministic model collection. In light of the success of variable selection via ℓ 1 -penalization in regression, they rather construct a random model collection derived from a collection of potentially relevant variable subsets determined by a sparse procedure. The relevant variable selection and the parameter estimation are performed in the same process.

For all K ∈ N * , let S K = S (K,{1,...,p}) = {s θ ; θ ∈ Θ K } according to (1) with

Θ K := Π K ×(R p ) K × R * + .
To detect the relevant variables, Pan and Shen (2007) penalize the empirical contrast γ n (s θ ) by a ℓ 1 -penalty on the mean parameters proportional to |θ| 1 := p j=1 K k=1 |µ kj |. For all λ ∈ G K , which is a chosen grid of regularization parameters, the Lasso estimator is defined by

θL (K,λ) = argmin θ∈Θ K {γ n (s θ ) + λ|θ| 1 } . (3)

Inria

To compute θL (K,λ) , [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] construct an EM algorithm for ℓ 1penalized model-based clustering. Then the relevant variable subset selected by the Lasso estimator θL (K,λ) is

J (K,λ) = {j ∈ {1, . . . , p}; ∃ k ∈ {1, . . . , K} such that μkj = 0}
and the density s is estimated by the Lasso solution ŝL (K,J (K,λ) ) := s θL

(K,λ)
.

By varying λ ∈ G K and K ∈ N * , they get a model collection

{S (K,Jr ) } (K,Jr )∈M L where M L = {(K, J r ); K ∈ N * , J r ∈ J K }, J K = λ∈G K J (K,λ) and S (K,Jr ) is defined by (1).
Next, the BIC criterion is used to solve the model selection problem. First, models are grouped according to their dimension D in order to obtain a model collection {S D } D∈D . For each dimension D ∈ D, let ŝL D be the maximum likelihood estimator in S D and (K D , J D ) such that ŝL D = ŝL (K D ,J D ) . Second, the minimizer D of the BIC criterion

D = argmin D∈D γ n (ŝ L D ) + 1 2 D n ln n (4)
is determined and ( K, Ĵr ) = (K D , J D ) is selected. Finally, a data clustering is derived from the estimated Lasso parameter vector θL ( K, Ĵr ) by applying the MAP principle.

Our Lasso-MLE procedure

Motivation

From our point of view, taking advantage of the sparsity property of ℓ 1 -penalization to perform automatic variable selection in clustering for high-dimensional data is an interesting idea which is worth exploring. Nevertheless, although the model-based clustering is linked to the density estimation problem, the PS-Lasso procedure does not take the density estimation purpose into account. Let us explain the main weakness of the PS-Lasso procedure as regards estimation.

Let S be the set of all densities with respect to Lebesgue measure on R p . In a maximum likelihood approach, the loss function considered is the Kullback-Leibler information defined for all t ∈ S by KL(s, t) = R p ln s(y) t(y) s(y) dy if s dy is absolutely continuous with respect to t dy and +∞ otherwise. The density s is the unique minimizer of the Kullback-Leibler information on S.

Ideally, we want to estimate s by the so-called oracle ŝ(K ⋆ ,Jr ⋆ ) where

(K ⋆ , J r ⋆ ) = argmin (K,Jr )∈M KL(s, ŝ(K,Jr) ). (5)
In practice, (K ⋆ , J r ⋆ ) depending on the unknown density s, ŝ(K ⋆ ,Jr ⋆ ) is unattainable. But it is a benchmark to evaluate the quality of any estimator of s. In particular, it is the benchmark for MM-MLE procedure. Yet for high-dimensional data, the collection M is so rich that performing the selection of the exhaustive best subset over M is unfeasible. In this case, a natural idea is to consider a subset M ′ ⊂ M so that performing best subset selection over M ′ becomes practicable by aiming at argmin (K,Jr )∈M ′ KL(s, ŝ(K,Jr) ). [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] construct a subset M ′ = M L of M. But, rather than considering the family of MLE estimator {ŝ (K,Jr ) } (K,Jr )∈M L by aiming at argmin

(K,Jr )∈M L KL s, ŝ(K,Jr) , (6)
they consider the family of Lasso estimators ŝL (K,Jr ) (K,Jr )∈M L and aim at argmin

(K,Jr )∈M L KL s, ŝL (K,Jr ) . (7)
Thus the aim of the proposed procedure is to mimic the model defined by Equation ( 6), which is expected to be closer to the oracle (5) than the model defined by (7).

Description of our procedure

We propose the so-called Lasso-MLE procedure which is a compromise between the two previous procedures (see Sections 3.1 and 3.2). This procedure is decomposed into three main steps. The first step consists of constructing a subcollection of models. As [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF], an ℓ 1 -approach is considered to obtain {S (K,Jr ) } (K,Jr )∈M L . For a fix number of mixture components K and a regularization parameter λ, an EM algorithm is used to compute the Lasso solution. It differs from the usual EM algorithm (for computing the ML estimator) by the update of mean parameters. In the ℓ 1 procedures, the choice of the regularization parameter grid is often a difficulty. Contrary to [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF], which opt for a deterministic grid difficult to choose, we propose to construct a data-driven grid of regularization parameters by using the updating formulas of the mixture parameters in the EM algorithm computing the Lasso solutions. This construction is detailed in Appendix C. The second step consists of computing the MLE ŝ(K,Jr) using the standard EM algorithm for each model (K, J r ) ∈ M L . The third step is devoted to model selection. As in Maugis and Michel (2011a,b), a non asymptotic penalized criterion is proposed to solve the model selection problem. Its construction requires to extend the theoretical result to determine the penalty shape in the high-dimensional context and with a random model subcollection. Next, the penalty depending on unknown multiplicative constant(s) is calibrated using the slope heuristics. The construction of this penalized criterion is the topic of Section 4.

Model selection

An oracle inequality for the Lasso-MLE estimator

In the third step of the Lasso-MLE procedure, a model selection criterion is required to select the number of clusters K and the relevant variables subset

⇒ model collection {S (K,Jr ) } (K,Jr )∈M L
Table 1: Summary of the three procedures: PS-Lasso , MM-MLE and Lasso-MLE. K (resp. J r ) denotes the number of clusters (resp. a relevant variable subset). J r simultaneously. We follow the approach developed by [START_REF] Birgé | From model selection to adaptive estimation[END_REF] and [START_REF] Barron | Risk bounds for model selection via penalization[END_REF] which consists of defining a non asymptotic penalized criterion leading to an oracle inequality. In the context of density estimation, [START_REF] Barron | Risk bounds for model selection via penalization[END_REF] and Massart (2007, Theorem 7.11) propose a general model selection theorem for maximum likelihood estimation. Since this theorem is stated for a deterministic model collection, an extension is proposed for a random model subcollection (see Theorem 2 in Appendix A). Then, by applying Theorem 2 to the considered random model collection of the finite Gaussian mixtures, the following oracle inequality for the Lasso-MLE estimator is derived.

Theorem 1. Let {S (K,Jr ) } (K,Jr )∈M be the model collection defined by ( 1). Let M L be a random subcollection of index sets (selected by the Lasso) included in the whole collection M. Let A µ , a σ and A σ be absolute positive constants and consider the collection of bounded models {S B (K,Jr ) } (K,Jr )M L defined by

S B (K,Jr ) = s θ ∈ S (K,Jr ) ; θ ∈ Π K × [-A µ , A µ ] |Jr | K × [a σ , A σ ] . ( 8
)
Consider the maximum likelihood estimator

ŝ(K,Jr) = argmin s θ ∈S B (K,Jr ) γ n (s θ ).
Denote by

D (K,Jr ) = K(1 + |J r |) the dimension of the model S B (K,Jr ) . Define B(A µ , A σ , a σ , p) := 1 + ln A σ a σ 1 + A µ a σ + ln p. Let s (K,Jr ) ∈ S B (K,Jr ) such that KL(s, s (K,Jr ) ) ≤ 2 inf s θ ∈S B (K,Jr ) KL(s, s θ ) and let τ > 0 such that s (K,Jr ) ≥ e -τ s. ( 9
)
Let pen : M → R + . Suppose that there exists an absolute constant κ > 0 such that, for all (K, J r ) ∈ M,

pen (K, J r ) ≥ κ D (K,Jr ) n B 2 (A µ , A σ , a σ , p) + ln 1 1 ∧ B 2 (A µ , A σ , a σ , p) D (K,Jr ) n +(1 ∨ τ ) ln p D (K,Jr ) ∧ p . ( 10
)
Then, the estimator ŝ( K, Ĵr ) with

( K, Ĵr ) = argmin (K,Jr )∈M L γ n (ŝ (K,Jr ) ) + pen(K, J r ) satisfies E d 2 H s, ŝ( K, Ĵr ) ≤ C E inf (K,Jr )∈M L inf s θ ∈S B (K,Jr ) KL(s, s θ ) + pen(K, J r ) + 1 ∨ τ n (11)
for some absolute positive constant C.

The proof of Theorem 1 is given in Appendix B. Note that Condition ( 9) is required to control the second moment of loglikelihood ratios in order to apply Bernstein's inequality to bound the empirical process of ln(s/s (K,Jr ) ) (see Lemma 1 in Appendix A.). The larger the parameter τ , the larger the minimal penalty (10) and the less accurate Inequality (11). Since s (K,Jr ) is positive, there always exists some τ > 0 fulfilling Condition (17) in Theorem 2. It seems difficult to have an idea of the minimal convenient value of τ since it depends on the unknown true density s. Nonetheless, we may think that Condition (9) is satisfied for reasonable values of τ because s (K,Jr ) is expected to be close to s. Note that the constant 2 in the Kullback-Leibler constraint for s (K,Jr ) can be replaced by 1 + ε with ε > 0.

As in [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF], the mixture parameters are bounded (see the model collection ( 8)) in order to construct brackets over S (K,Jr ) and thus to upper bound the entropy number. The bracket construction is adapted from the one of [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] for our specific Gaussian mixtures (see Appendix B). We also obtain an Inequality (11) which is not exactly an oracle inequality since the Hellinger risk is upper bounded by the Kullback-Leibler bias. But contrary to MM-MLE procedure, our Lasso-MLE procedure runs on a small random subcollection of models and it remains feasible even for large p. Thus our estimator ŝ(K,Jr) is attainable in practice. Moreover, contrary to classical asymptotic criteria for which p is fixed and n tends to infinity, our result is non asymptotic and allows to study cases for which p increases with n. Since the ratio ln(p)/n appears in the right hand-side of Inequality (11) through the term pen(K, J r ), Theorem 1 ensures that our estimator achieves Inria good performance compared with the oracle as long as p < e n , which allows to consider many situations with p ≫ n.

As expected, the penalty ( 10) is proportional to the model dimension and thus penalizes models with high complexities. It also involves two additional logarithm terms. On the one hand, the first logarithm term is due to a lack of accuracy in the proof of Theorem 1. Specifically, only a global entropy bracketing control is obtained while a local version is only required for applying Theorem 2. It is sufficient since the local entropy is upper bounded by the global entropy but it is not optimal and yields extra logarithm terms. On the other hand, the second logarithm term quantifies the complexity of the model collection by taking into account the possible large number of models with identical dimension. In regression, [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] prove that, for complete variable selection, a penalty proportional to the dimension selects too complex models with high probability and that a logarithm term is necessary to select smaller models. This has been practically checked in many situations: for multiple change points detection in a regression framework [START_REF] Lebarbier | Detecting multiple change-points in the mean of Gaussian process by model selection[END_REF] or for histogram selection in a density estimation framework [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF] for instance. Nevertheless this logarithm term becomes unnecessary if the number of models with the same dimension is small enough. For instance, for finite Gaussian mixture models in a low-dimensional setting, Maugis and Michel (2011a) observe that a penalty proportional to the dimension,with no logarithm term, is sufficient to select a model close to the oracle. But in our high-dimensional context, the number of models having the same dimension is expected to grow. Nonetheless, thanks to the random preselection of relevant variables subsets, a complete variable selection is not performed here. However, it is difficult to know how rich is the random model collection. Thus, according to (10), we may retain that the penalty is pen(K, J r ) = κ 1 pen shape (K, J r ) with

pen shape (K, J r ) = D (K,Jr ) n 1 + κ 2 ln p D (K,Jr ) (12)
where κ 1 and κ 2 are two unknown constants. But if our random model subcollection is much poorer than the whole model collection and contains few models with the same dimension, the penalty (12) may be too pessimistic. In this case, a penalty proportional to the dimension

pen shape (K, J r ) = D (K,Jr ) n ( 13
)
might be sufficient to select a model with proper dimension.

Practical penalty calibration

According to ( 12) and ( 13), the penalty is known up to multiplicative constant(s) and in practice, we have to choose between these two penalty shapes. To fill in the gap between theory on penalization and practical calibration penalty, [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] proposed the so-called slope heuristics. In the context of Gaussian homoscedastic least squares regression with fixed design, they show that there exists a minimal penalty, namely such that the dimension and the risk of the models selected with smaller penalties become large. Moreover, they propose the rule of thumb of the slope heuristics: a penalty equal to twice the minimal penalty allows us to select a model close to the oracle model. This rule of thumb is theoretical proved in few specific frameworks [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF][START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF][START_REF] Arlot | Data-driven calibration of linear estimators with minimal penalties[END_REF][START_REF] Lerasle | Optimal model selection for stationary data under various mixing conditions[END_REF][START_REF] Lerasle | Optimal model selection in density estimation[END_REF] and largely valid in practice (see for instance [START_REF] Lebarbier | Detecting multiple change-points in the mean of Gaussian process by model selection[END_REF]Maugis and Michel, 2011a;[START_REF] Caillerie | Model selection for simplicial approximation[END_REF][START_REF] Verzelen | Data-driven neighborhood selection of a Gaussian field[END_REF]. For determining the minimal penalty, the data-driven slope estimation (DDSE) procedure proposed by [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] is used. This method is based on the existence of a linear behavior between the penalty shape and the contrast value for the most complex models (the empirical bias gets stable for the most complex models and the minimal penalty, corresponding to an empirical "estimation error" term, is given by the behavior of -γ n (ŝ (K,Jr ) )). In our context, the DDSE procedure first provides a graphical way to preliminary choose the penalty shape between ( 12) and ( 13), and second allows us to calibrate the penalty by estimating the linear slope. The reader is referred to [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] for more details about the DDSE procedure. Note that if the penalty shape ( 12) is required, the DDSE procedure is adapted with a double regression for calibrating the two constants κ 1 and κ 2 (instead of one constant usually).

Applications

In this section, the Lasso-MLE procedure is studied on two simulated examples and is compared with PS-Lasso and MM-MLE procedures. The Adjusted Rand Index (ARI) is used to measure the similarity between data clusterings and the variables declared relevant for the clustering by the three procedures are compared. We also compare the data-driven model selection criterion with two widely-used criteria: the Akaike Information Criterion (AIC) of [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] and defined by γ n (ŝ D )+ D n ; the Bayesian Information Criterion (BIC) of [START_REF] Schwarz | Estimating the dimension of a model[END_REF] defined by γ n (ŝ D ) + ln n 2 D n . Moreover, since the density s is known for simulated datasets, the model selected by one criterion can be compared with the associated oracle model to judge the quality of this criterion. The oracle model is defined by

D oracle = argmin D∈D KL (s, ŝD ) = argmin D∈D - x∈R p ln [ŝ D (x)] s(x) dx .
In practice, the oracle model is obtained by approximating the integral by a Monte Carlo procedure.

First simulated example

This first simulated dataset is in the spirit of an example in Maugis and Michel (2011a). The dataset consists of n = 200 observations described by p variables. The data are simulated according to a mixture of four spherical Gaussian distributions Φ(• | µ k , I) where µ 1 = -µ 3 = (3, 2, 1, 0.7, 0.3, 0.2, 0.1, 0.07, 0.05, 0.025, 0 p-10 ), µ 2 = 0 p , µ 4 = (3, -2, 1, -0.7, 0.3, -0.2, 0.1, -0.07, -0.05, -0.025, 0 p-10 ).

Inria

The vector 0 l denotes the null vector of length l. The mixing proportions are (π 1 , π 2 , π 3 , π 4 ) = (0.3, 0.2, 0.2, 0.3). The relevant clustering variables are the first ten variables (J r ⋆ = {1, . . . , 10}). Note that the four subpopulations are progressively gathered together into a unique Gaussian distribution, thus the discriminant power of the relevant variables decreases with respect to the variable index (see Figure 1). In this simulation study 20 datasets are simulated for each value of p ∈ {30, 200, 1000}, and mixtures with K ∈ {2, . . . , 6} components are considered.

G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 -4 -2 0 2 4 6 VAR1 VAR2 VAR11 VAR10 VAR9 VAR8 VAR7 VAR6 VAR5 VAR4 VAR3
Figure 1: Boxplots of the first eleven variables (VAR1,...,VAR11) on the four mixture components (G1,G2,G3,G4).

The Lasso-MLE procedure is compared with PS-Lasso and MM-MLE procedures in the low-dimensional case (p = 30) and only with PS-Lasso procedure in the high-dimensional cases (p = 200, 1000). Note that for the MM-MLE procedure, the variables are preliminary ordered (the collection of relevant variable subsets is {(1, 2, . . . , d), 1 ≤ d ≤ p}) as in Maugis and Michel (2011a) because the model collection is too rich to perform complete variable selection. Figure 2 checks that the function D/n ∈ [0, 1] → -γ n (ŝ D ) has a linear behavior for most complex models, which justifies the use of penalty shape (13). Note that the values of the estimated slopes by DDSE for Lasso-MLE and MM-MLE procedures are very similar here. The results are summarized in Table 2.

Low-dimensional case (p = 30 ≪ n)

With the Lasso-MLE procedure, the data-driven criterion (DDSE), BIC and the oracle globally select the true number of cluster K ⋆ = 4. With the variable selection, the DDSE criterion allows us to select a relevant variable subset closer to the oracle than BIC. Moreover, note that the model collections considered by both procedures are different. With the MM-MLE procedure, the ordered variable subsets is considered, containing the true relevant variable subset J r ⋆ . This explains in particular why the selected models never contains false relevant variables. On the contrary, the Lasso-MLE procedure considers data-driven relevant variable subsets for each simulation. In particular, J r ⋆ may not belong to this data-driven collection. Moreover, the relevant variable subsets collection obtained by varying the regularization parameter in the Lasso procedure may remove the least relevant variables (j ∈ {7, . . . , 10}, see Figure 1) before some noisy variables. With both procedures, the oracle model does not coincide with the true model. Furthermore, AIC often selects too many components and more relevant variables than the other criteria.

With the PS-Lasso procedure, the oracle model sometimes overestimates the number of mixture components. It contains most (yet not all) true relevant variables but also many false relevant variables. This tendency to select too many variables has already been widely noted in the regression framework [START_REF] Zhao | On model selection consistency of lasso[END_REF][START_REF] Zou | The adaptive lasso and its oracle properties[END_REF][START_REF] Yuan | On the non-negative garrotte estimator[END_REF][START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF][START_REF] Connault | Calibration d'algorithmes de type Lasso et analyse statistique de données métallurgiques en aéronautique[END_REF]. For this procedure, the errors of variable selection are reduced by using BIC.

High-dimensional case (p = 200, 1000)

For the two high-dimensional scenarios where p = 200 and p = 1000, the Lasso-MLE procedure and PS-Lasso procedure are compared. First, results in Table 2 show that the variable selection problem for the PS-Lasso procedure (oracle and BIC) gets worse when p grows: the selected variable subset contains more and more false relevant variables at the expense of true relevant variables. Moreover, the selection of the component number is deteriorated. With the Lasso-MLE procedure, the models selected by AIC and, to a lesser extend by BIC, contain more (true relevant and false relevant) variables when p grows. On the opposite, the oracle model and the model selected by the data-driven penalized criterion remain stable. Moreover, the average multiplicative factor 2κ in the data-driven penalty pen(D) = 2κ D n is equal to 1.47(±0.10), 2.27(±0.19) and 3.68(±0.57) when p = 30, 200 and 1000 respectively. This factor globally increases with respect to p and thus the associated penalty becomes stronger as p increases, whereas the fixed BIC penalty (0.5 ln(n) = 2.65) tends to underpenalize as p grows.

Globally the data clustering slightly deteriorates as p grows. The biggest deterioration is for PS-Lasso procedure. The ARI for each model in PS-Lasso and Lasso-MLE procedures is represented in Figure 3. For PS-Lasso procedure, the models achieving the best clustering have moderate or high dimension while BIC selects a less complex model. On the opposite, for Lasso-MLE procedure, small models achieve good clustering. In particular, the DDSE, which selects such a small model, leads to a satisfactory data clustering. Figure 3

Second simulated example

This second simulated example is proposed in [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF]. The dataset consists of n = 200 observations described by p = 1000 variables. The data are simulated according to a mixture of two Gaussian distributions π 1 Φ(•|0 p , I) + (1π 1 )Φ(•|µ 2 , I) where µ 2 = (1.5, . . . , 1.5, 0 950 ) and π 1 = 0.85. The relevant variables are the first fifty variables (J r ⋆ = {1, . . . , 50}). We perform 20 simulations of the dataset. For each simulation, models with K ∈ {1, 2, 3} clusters are considered. In this high-dimensional context, the Lasso-MLE procedure is compared with the PS-Lasso procedure. The results are summarized in Table 3.

Table 3 shows that the Lasso oracle model, and to a lesser extend the model selected by BIC, contain many false relevant variables and may overestimate the number of mixture components. This confirms that the PS-Lasso procedure is not suited to recover the true model and the true relevant variables. Moreover, BIC data clustering is disappointing. In contrast, the Lasso-MLE oracle model always coincides with the true model and leads to a very good data clustering. The DDSE achieves better performance than BIC and AIC. Table 3: Averaged number of true relevant (TR) and false relevant (FR) variables (± standard deviation); number of times a clustering with K = 1, 2 and 3 components is selected; Averaged ARI (± standard deviation) over the 20 simulations.

Extension for curve clustering

For some clustering problems, one may be interested in providing a sparse parameter estimation for each cluster besides clustering the data. This is typically the case for curve clustering problems involving a smooth representative curve per cluster (see for instance Misiti et al., 2007a;[START_REF] Auder | Projection-based curve clustering[END_REF]. This sparse curve reconstruction allows to improve the curve clustering interpretability and may be used for prediction. For instance, this problem is encountered in the study of electricity consumption where representative curves are desirable to forecasting. For such problems, curves are preliminary decomposed into some appropriate basis such as a wavelet basis. Then, the data are the coefficients associated to the curve decomposition into the basis. If the Lasso-MLE procedure is used on the empirically centered coefficient matrix to obtain a curve partition, representative curves per cluster are not directly available. One solution to estimate a representative profile for a given cluster would be to take the mean of the ob-served curves assigned to this cluster. Nonetheless, this process produces noised profiles and it is not able to provide smooth profiles. To overcome this problem, an alternative to the Lasso-MLE procedure which avoids empirical centering of the data and which is suited to clustering involving smooth representative curves is proposed. This alternative is based on a more general variable role modeling, implying a richer model collection than (1). The theoretical model selection criterion is adapted and we highlight that a logarithm factor in the penalty is now observed, due to the model collection richness.

Description of the alternative procedure

In this section, data are not assumed to fulfill the assumption E[Y ij ] = 0. A variable j is now irrelevant for the clustering if µ 1j = . . . = µ Kj := ν j , ν j being not necessary equal to zero. With the previous modeling, given a number K of clusters and relevant variables subset J r , K|J r | + |J c r | ≥ p mean parameters have to be estimated, even when p ≫ n.

Therefore, selecting the relevant clustering variables J r is not sufficient to get sparse models and an additional dimensional reduction step is needed before the estimation step. Here, the irrelevant variables J c r are assumed to be decomposed into the zero irrelevant variables J 0 (j ∈ J 0 if ν j = 0) and non zero irrelevant variables J c 0 (j ∈ J c 0 if ν j = 0). Moreover, we assume that the irrelevant variables are predominantly zero irrelevant (|J 0 | ≫ |J c 0 |). In particular, this assumption is fulfilled for curve clustering using sparse representations of signals in some appropriate basis such as a wavelet basis. This new variable role modeling leads to the new model collection {S (K,Jr ,J0) } (K,Jr ,J0)∈M LL where S (K,Jr ,J0) is defined by

         y ∈ R p → s θ (y); s θ (y) = Φ(y [J0] |0, σ 2 I) Φ(y [J c 0 ] |ν, σ 2 I) K k=1 π k Φ(y [Jr ] |µ k , σ 2 I) θ = (π 1 , . . . , π K , ν, µ 1 , . . . , µ K , σ) ∈ Π K × R |J c 0 | × R |Jr | K × R * +          (14)
The dimension of the model

S (K,Jr ,J0) is D (K,Jr ,J0) = K(1 + |J r |) + |J c 0 |.
In practice, the random model collection S (K,Jr ,J0) is determined using the following alternative procedure. First, a model collection {S (K,Jr ) } (K,Jr )∈M L is constructed using the Lasso-MLE procedure on the centered data. Second, for each J c r , the empirical contrast restricted to the Y i[J c r ] 's is penalized by an ℓ 1 -penalty proportional to ν 1 with various regularization parameter values, in order to detect the zero irrelevant variables. Next, Theorem 1 can be easily adapted for the model collection {S (K,Jr ,J0) } (K,Jr ,J0)∈M LL in order to construct a data-driven non asymptotic penalized criterion to solve the model selection problem. As in Section 4, the same sufficient penalty shape ( 10) is established with model dimensions D (K,Jr ,J0) . But this alternative procedure depending on two embedding Lasso algorithms, the model collection is richer than the model collection {S (K,Jr ) } (K,Jr )∈M L obtained with the Lasso-MLE procedure. As a consequence, the logarithm term, which takes into account the richness of the model collection and the number of models having the same dimension, is now observed in the penalty shape (see Figure 5). As explained in Section 4.2, the penalty is calibrated using the DDSE method.

Inria

Functional data clustering example

The functional dataset consists of n = 200 noisy curves simulated as follows: the wavelet coefficients {Y ij , i = 1, . . . , n, j = 1, . . . , p} with p = 1086 are simulated according to the Gaussian mixture 0.85 Φ(• | µ 1 , I) + 0.15 Φ(• | µ 2 , I) with mean vectors µ 1 = (0 25 , 1.5 25 , 0 p-50 ) and µ 2 = (1.5 50 , 0 p-50 ) (a l denotes the vector of length l whose coordinates equal a). Each wavelet coefficients vector Y i corresponds to the decomposition of a curve g i in the symmlet-4 basis at level 10 ( [START_REF] Misiti | Wavelets and their Applications[END_REF]. In particular, µ 1 and µ 2 are the decomposition of f 1 and f 2 which are the discretization of two functions on a fine time grid {t 1 , . . . , t T } containing T = 1024 points, represented at Figure 4 (top left). Thus, the first 25 variables are relevant for the clustering while the variables 26 to 50 are non zero irrelevant. In order to obtain a clustering of this 200 curves, our alternative method is used on the wavelet coefficients. Next, an inverse wavelet transform of the estimated mean vectors is performed to obtain a denoised representative curve for each cluster. Models with K ∈ {1, 2, 3} mixture components are considered. The alternative procedure is compared with the Lasso-MLE procedure. Both procedures select a model with two classes. The two estimated representative curves obtained by the Lasso-MLE procedure with model collection {S (K,Jr ) } (K,Jr )∈M L described in Section 3.3.2 and by the alternative procedure with model collection {S (K,Jr ,J0) } (K,Jr ,J0)∈M LL described in Section 6.1 are displayed in Figure 4. As expected, the curve estimations obtained by the

Inria

Lasso-MLE procedure is very noisy because no thresholding of the smallest mean coefficients is performed. Note that the PS-Lasso procedure suffers from the same drawback and moreover, it gives shrunk curve profile estimations since the non zero mean wavelet coefficients are shrunk due to ℓ 1 -penalization. On the contrary, the alternative procedure leads to smoother estimated representative curves. This is all the more true for the estimations obtained by the penalty shape with the logarithm term rather than without the logarithm term. This highlights that the model collection {S (K,Jr ,J0) } (K,Jr ,J0)∈M LL is rich enough to require a penalty with a logarithm factor and that a penalty proportional to the dimension underpenalizes, resulting in a less smooth curve estimation with some extra peaks. The slope graphs obtained for one simulation of this dataset are represented in Figure 5.

Discussion

In this paper, the Lasso-MLE procedure is proposed to cluster data, detecting the relevant clustering variables. This procedure is especially suited to study high-dimensional datasets. It is based on a compromise between a ℓ 1regularization procedure and a MLE procedure. The variable selection and the data clustering problems are recast into a general model selection problem. A ℓ 1 -approach is used to perform automatic variable selection and thus to deduce a reasonable random model collection. A data-driven penalized criterion is then built to solve the model selection problem and a data clustering is deduced to the associated ML estimator using the MAP principle.

Only spherical Gaussian mixtures are considered in this paper to focus on the mean vectors of mixtures. Nevertheless the Lasso-MLE procedure may be extended for Gaussian mixtures with more general variance matrices. First a more general Lasso procedure can be adapted to the ℓ 1 -procedure of [START_REF] Zhou | Penalized model-based clustering with unconstrained covariance matrices[END_REF]. Second, Theorem 1 may be extended for a random collection of general Gaussian mixture models by using the control of the bracketing entropy for the general Gaussian mixture families established in [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF]. The main difficulty lies in the definition of a relevant clustering variable, which depends on the behavior of the mean vectors but also the definition of variance matrices for the studied mixtures. This could have consequences on the penalty shape in particular. This extension will be the topic of a future work.

This use of a ℓ 1 -approach to only produce a reasonable collection of relevant variable subsets is an interesting idea which may be beneficial for other variable selection procedures. For instance, [START_REF] Maugis | Variable selection in model-based clustering: A general variable role modeling[END_REF] propose a procedure based on a general variable role modeling but the associated algorithm Selvar-ClustIndep is too slow when the number of variables increases. A ℓ 1 variable selection procedure may be used to initialize the forward version of this algorithm to improve its stability and/or give a way of variable subsets to follow in order to dramatically speed up the algorithm.

A A model selection theorem for MLEs in random models

The Hellinger distance between two nonnegative integrable functions t and u is the norm √ t -√ u , denoted d H (t, u). Consider S the set of all densities on R p . An ε-bracketing for a subset S of S with respect to d H is a set of integrable function pairs (l 1 , u 1 ), . . . , (l N , u N ) such that for each t ∈ S, there exists j ∈ {1, . . . , N } such that l j ≤ t ≤ u j and d H (l j , u j ) ≤ ε. The bracketing number N [.] (ε, S, d H ) is the smallest number of ε-brackets necessary to cover S and the bracketing entropy is defined by

H [.] (ε, S, d H ) = ln[N [.] (ε, S, d H )].
Let {S m } m∈M be some at most countable model collection such that S m ⊂ S for all m ∈ M. We shall say that {S m } m∈M fulfills Property (P) if, for all m ∈ M, H [.] (ε, S m , d H ) is integrable at 0 and if there exists a function

Ψ m on R + such that Ψ m is nondecreasing, ξ → Ψ m (ξ)/ξ is nonincreasing on ]0, +∞[, and for ξ ∈ R + and u ∈ S m , denoting S m (u, ξ) = {t ∈ S m ; d H (t, u) ≤ ξ}, ξ 0 H [.] (ε, S m (u, ξ), d H ) dε ≤ Ψ m (ξ). (15)
Theorem 2. Let s ∈ S be an unknown density to be estimated from a n-sample (Y 1 , . . . , Y n ). Consider {S m } m∈M some at most countable deterministic model collection fulfilling Property (P). Let {x m } m∈M be some family of nonnegative numbers such that

m∈M e -xm = Σ < ∞. ( 16
)
For every m ∈ M, consider Ψ m defined by (P) and ξ m such that Ψ m (ξ m ) = √ nξ 2 m . Let s m ∈ S m such that KL(s, s m ) ≤ 2 inf t∈Sm KL (s, t) and let τ > 0 such that s m ≥ e -τ s.

(

) 17
Introduce {S m } m∈ M some random subcollection of {S m } m∈M . Let ρ ≥ 0 and consider the collection of ρ-MLEs {ŝ m } m∈ M :

γ n (ŝ m ) ≤ inf t∈Sm γ n (t) + ρ.
Let pen : M → R + . Suppose that there exists an absolute constant κ > 0 such that, for all m ∈ M,

pen (m) ≥ κ ξ 2 m + (1 ∨ τ ) x m n . ( 18
)
Let ρ ′ ≥ 0. Then, any penalized likelihood estimator ŝ m with m ∈ M such that

γ n (ŝ m) + pen ( m) ≤ inf m∈ M {γ n (ŝ m ) + pen (m)} + ρ ′ (19) satisfies E d 2 H (s, ŝ m) ≤ C E inf m∈ M inf t∈Sm KL(s, t) + pen(m) + (1 ∨ τ ) Σ 2 n + ρ + ρ ′
(20) for some absolute positive constant C.

Inria

To prove Theorem 2, we provide an inequality for the moments of order 2 of loglikelihood ratios. This inequality is based on the following Claim 3. Claim 3. Let τ > 0. For all x > 0, consider f (x) = x(ln x) 2 ,, h(x) = x ln xx + 1 and φ(x) = e xx -1. Then, for all 0 < x ≤ e τ ,

f (x) ≤ τ 2 φ(-τ ) h(x). ( 21
)
Proof. First note that f (1) = h(1) = 0, thus we just need to prove (21) for x = 1. Define

ψ : R → R, y → φ(y)/y 2 if y = 0, 1/2 if y = 0 and ϕ : R → R, y → φ(y)/y if y = 0, 0 if y = 0.
Let us first check that ψ is nondecreasing on R. Since e y = 1 + y + y 2 /2 + o y→0 (y 2 ), the functions ψ and ϕ are continuous on R and, for y = 0, ψ(y) = (ϕ(y)ϕ(0))/(y -0) is the difference quotient at 0 for ϕ. Thus, we just need to prove that ϕ is a convex function to derive that ψ is nondecreasing. By differentiating twice ϕ, we get that ϕ ′′ (y) = 2e y g(y)/y 3 with g(y) = 1y + y 2 /2e -y . The function g is nondecreasing because g ′ (y) = -1 + y + e -y ≥ 0. But g(0) = 0. So, g(y) ≤ 0 for all y ≤ 0 and g(y) ≥ 0 for all y ≥ 0. It implies that ϕ ′′ (y) ≥ 0 for all y ∈ R and ϕ is convex. Now, let 0 < x ≤ e τ , x = 1. Put y =ln x. Then, y ≥ -τ and, since ψ is nondecreasing, ψ(y) ≥ ψ(-τ ). Moreover, x = 1, so y = 0 and ψ(y) = φ(y)/y 2 . Thus, φ(y)

y 2 ≥ φ(-τ ) τ 2 . ( 22
)
Taking the definition of φ and y =ln x into account, (22) leads to

ln x -1 + 1 x ≥ φ(-τ ) τ 2 (ln x) 2 .
We get (21) by multiplying the last inequality by x > 0.

Lemma 1. Let P and Q be two probability measures with P ≪ Q. Assume that there exists τ > 0 such that ln ( dP/dQ ∞ ) ≤ τ . Then,

ln dP dQ 2 dP ≤ τ 2 e -τ + τ -1 KL(P, Q). ( 23
)
Proof. Since ln (dP/dQ) ≤ τ , we can apply Claim 3 to x = dP/dQ:

f dP dQ ≤ τ 2 φ(-τ ) h dP dQ .
Integrating with respect to Q and taking the definition of f , φ and h into account, we get

dP dQ ln dP dQ 2 dQ ≤ τ 2 e -τ + τ -1 ln dP dQ dP -dP + dQ thus ln dP dQ 2 dP ≤ τ 2 e -τ + τ -1 KL(P, Q).
Proof of Theorem 2. For the sake of simplicity, we assume that ρ = ρ ′ = 0.

For any measurable function g, denote by ν n the recentred process defined by

ν n (g) = 1 n n i=1 {g(Y i ) -E[g(Y i )]} . ( 24
)
For all m ∈ M, consider s m such that KL(s, s m ) ≤ 2 inf t∈Sm KL(s, t). Define the functions

g m = - 1 2 ln s m s , ĝm = - 1 2 ln ŝm s , fm = -ln s + ŝm 2s . ( 25
) Fix m ∈ M. Introduce M(m) = {m ′ ∈ M, γ n (ŝ m ′ ) + pen (m ′ ) ≤ γ n (ŝ m ) + pen (m)}
and let m ′ ∈ M(m). By definition of ŝm ,

γ n (ŝ m ′ ) + pen (m ′ ) ≤ γ n (ŝ m ) + pen (m) ≤ γ n (s m ) + pen (m) ,
and according to the definition of g m and ĝm ′ ,

2 n n i=1 ĝm ′ (Y i ) + pen (m ′ ) ≤ 2 n n i=1 g m (Y i ) + pen (m) . (26)
Then, by concavity of the logarithm, we have fm ′ ≤ ĝm ′ , and (26) gives

2 n n i=1 fm ′ (Y i ) ≤ 2 n n i=1 g m (Y i ) + pen (m) -pen (m ′ ) . (27)
By taking ( 24) and ( 25) into account, Inequality (27) implies

2 KL s, s + ŝm ′ 2 ≤ KL (s, s m ) + pen (m) -pen (m ′ ) + 2 ν n (g m ) -ν n fm ′ .
(28) Our purpose is now to control both ν n (g m ) and -ν n ( fm ′ ).

To bound -ν n ( fm ′ ), we refer to the proof of Theorem 7.11 in [START_REF] Massart | Concentration inequalities and model selection[END_REF]. It is proved that there exists κ ′′ > 0 such that for all u > 0, for all m ′ ∈ M(m), for all y > ξ m ′ , the following inequality holds except on a set with probability less than 2e -u :

-ν n ( fm ′ ) y 2 + √ s - √ ŝm ′ 2 ≤ κ ′′ ξ m ′ + u/n y + u ny 2 . ( 29
)
Let us now focus on controlling ν n (g m ). From ( 24) and ( 25), we have

ν n (g m ) = n i=1 X i -E[X i ], X i := 1 2n ln s(Y i ) s m (Y i ) . ( 30
)

Inria

To get an upper bound of ν n (g m ), we apply Bernstein's Inequality [START_REF] Massart | Concentration inequalities and model selection[END_REF]. This inequality requires to control the moments of order k for all k ≥ 2 of X i defined by (30). Such a control is provided by Lemma 1 on condition that ln( s/s m ∞ ) ≤ τ . Assume that ( 17) is fulfilled. Then, ln( s/s m ∞ ) ≤ τ and we deduce from Lemma 1 that

R p ln s(y) s m (y) 2 s(y) dy ≤ τ 2 e -τ + τ -1 KL(s, s m ).
On the one hand, τ 2 /(e -τ + τ -1) ∼ τ →∞ τ , so there exists A > 0 such that τ 2 /(e -τ + τ -1) ≤ 2τ for all τ ≥ A. On the other hand, τ → τ 2 /(e -τ + τ -1) is continuous on ]0, A] and τ 2 /(e -τ + τ -1) ∼ τ →0 2, so there exists B > 0 such that τ 2 /(e -τ +τ -1) ≤ B for all τ ∈]0, A]. Thus, for all τ > 0, τ 2 /(e -τ +τ -1) ≤ δ(1 ∨ τ ) with δ = 2 ∨ B, and

R p ln s(y) s m (y) 2 s(y) dy ≤ δ(1 ∨ τ ) KL(s, s m ). (31)
From ( 30), ( 31) and the assumption ln( s/s m ∞ ) ≤ τ , we derive that

n i=1 E X 2 i ≤ n (2n) 2 R p ln s(y) s m (y) 2 s(y) dy ≤ δ(1 ∨ τ ) KL(s, s m ) 4n (32)
and that for all integers ≥ 3,

n i=1 E (X i ) k + ≤ n (2n) k R p ln s(y) s m (y) k + s(y) dy ≤ n (2n) k R p ln s(y) s m (y) k 1 {s(y)≥sm(y)} s(y) dy ≤ n (2n) k R p ln s(y) s m (y) k-2 ln s(y) s m (y) 2 1 {s(y)≥sm(y)} s(y) dy ≤ n (2n) k τ k-2 R p ln s(y) s m (y) 2 1 {s(y)≥sm(y)} s(y) dy ≤ n (2n) k τ k-2 δ(1 ∨ τ ) KL(s, s m ) ≤ 1 2 τ 2n k-2 δ(1 ∨ τ ) KL(s, s m ) 2n . (33)
From ( 32) and ( 33), we can apply Bernstein's Inequality with

v := δ(1 ∨ τ ) KL(s, s m ) 2n , c := τ 2n . ( 34
)
It gives that, for every positive u, except on a set with probability less than e

-u , ν n (g m ) ≤ √ 2vu + cu. (35)
Let z > 0 to be chosen later. Using that z 2 + KL(s, s m ) ≥ 2z KL(s, s m ) and z 2 + KL(s, s m ) ≥ z 2 , we get from (35) that except on a set with probability less than e -u ,

ν n (g m ) z 2 + KL(s, s m ) ≤ √ 2vu + cu z 2 + KL(s, s m ) ≤ √ vu z 2 KL(s, s m ) + cu z 2 . ( 36
)
Let us gather ( 29) and ( 36): There exists κ ′′ > 0 such that, for every positive u, for all m ′ ∈ M(m), for all z > 0 and for all y ≥ ξ m ′ , except on a set with probability less than 3e -u ,

-ν n fm ′ y 2 + √ s - √ ŝm ′ 2 ≤ κ ′′ ξ m ′ y + u/n y + u ny 2 (37) and ν n (g m ) z 2 + KL(s, s m ) ≤ √ vu z 2 KL(s, s m ) + cu z 2 . ( 38
)
Now, let x > 0. Let x m and x m ′ be defined by ( 16). We apply ( 37) and ( 38) to u = x + x m + x m ′ and we choose adequately y and z by defining for some constants γ and β to be specified later,

y m,m ′ := γ -1 ξ 2 m ′ + x + x m + x m ′ n . (39)
and

z m,m ′ := β -1 v 2 KL(s, s m ) + c (x + x m + x m ′ ). (40)
Using that a 2 + b 2 ≥ a 2 , we get that except on a set with probability less than 3e

-(x+xm+x m ′ ) , -ν n ( fm ′ ) ≤ κ ′′ (2γ + γ 2 ) y 2 m,m ′ + √ s -ŝm ′ 2 (41) and ν n (g m ) ≤ (β + β 2 )(z 2 m,m ′ + KL(s, s m )). (42)
We can now come back to Inequality (28). Injecting ( 41) and ( 42) into (28) yields Massart, 2007, Lemma 7.23) and choosing γ such that 2κ ′′ (2γ + γ 2 ) = 2 ln 2 -1 := α, we get

2 KL s, s + ŝm ′ 2 ≤ KL (s, s m ) + pen (m) -pen (m ′ ) +2(β + β 2 )[z 2 m,m ′ + KL(s, s m )] +2κ ′′ (2γ + γ 2 ) y 2 m,m ′ + √ s -ŝm ′ 2 . Putting κ(β) := (1 + 2(β + β 2 )), using the inequality KL(s, (s + ŝm ′ )/2) ≥ (2 ln 2 -1) √ s - √ ŝm ′ 2 (
α √ s -ŝm ′ 2 ≤ κ(β) KL (s, s m ) + pen (m) -pen (m ′ ) + 2(β + β 2 )z 2 m,m ′ + αy 2 m,m ′ .

Inria

From ( 40), ( 39) and ( 34), we deduce that

α √ s -ŝm ′ 2 ≤ κ(β) KL(s, s m ) + pen(m) -pen(m ′ ) + (β + β 2 )β -2 δ(1 ∨ τ ) 2 + τ x + x m + x m ′ n + αγ -2 ξ 2 m ′ + x + x m + x m ′ n . Since τ ≤ 1 ∨ τ , if we choose β such that (β + β 2 )(δ/2 + 1) = αγ -2 , we get α √ s -ŝm ′ 2 ≤ κ(β) KL(s, s m ) + pen(m) -pen(m ′ ) + αγ -2 ξ 2 m ′ + αγ -2 β -2 (1 ∨ τ ) + 1 x + x m + x m ′ n . Put κ = αγ -2 (β -2 + 1). Then, since 1 ≤ 1 ∨ τ , α √ s -ŝm ′ 2 ≤ κ(β) KL(s, s m ) + pen(m) -pen(m ′ ) + αγ -2 ξ 2 m ′ + κ (1 ∨ τ ) x + x m + x m ′ n ≤ κ(β) KL(s, s m ) + pen(m) + κ (1 ∨ τ ) x m n + αγ -2 ξ 2 m ′ + κ(1 ∨ τ ) x m ′ n -pen(m ′ ) + κ(1 ∨ τ ) x n ≤ κ(β) KL(s, s m ) + pen(m) + κ (1 ∨ τ ) x m n + κ ξ 2 m ′ + (1 ∨ τ ) x m ′ n -pen(m ′ ) + κ(1 ∨ τ ) x n .
Now, assume that Condition (18) on the penalty function is fulfilled for this value of κ. Then, for all x > 0, for every m ∈ M and m ′ ∈ M(m), except on a set with probability less than 3e

-(x+xm+x m ′ ) , α √ s -ŝm ′ 2 ≤ κ(β) KL(s, s m ) + 2 pen(m) + κ(1 ∨ τ ) x n . (43)
It only remains to sum up the tail bounds (43) over all the possible values of m ∈ M and m ′ ∈ M(m) by taking the union of the different sets of probability less than 3e -(x+xm+x m ′ ) . For all x > 0, except on a set with probability less than

3 m∈M,m ′ ∈M(m) e -(x+xm+x m ′ ) ≤ 3e -x (m,m ′ )∈M×M e -(xm+x m ′ ) = 3e -x m∈M e -xm 2 = 3Σ 2 e -x ,
we have simultaneously for all m ∈ M and m ′ ∈ M(m),

α √ s -ŝm ′ 2 ≤ κ(β) KL(s, s m ) + 2 pen(m) + κ(1 ∨ τ ) x n . (44)
Inequality ( 44) is in particular satisfied for all m ∈ M and m ′ ∈ M(m) and, since m defined by ( 19) belongs to M(m) for all m ∈ M, we deduce from ( 44) that for all x > 0, except on a set with probability less than 3Σ 2 e -x ,

α √ s -ŝ m 2 ≤ inf m∈ M {κ(β) KL(s, s m ) + 2 pen(m)} + κ(1 ∨ τ ) x n ≤ inf m∈ M 2κ(β) inf t∈Sm KL(s, t) + 2 pen(m) + κ(1 ∨ τ ) x n . (45)
By integrating (45) over x > 0, we finally get that there exists an absolute constant C > 0 such that

E √ s -ŝ m 2 ≤ C E inf m∈ M inf t∈Sm KL(s, t) + pen(m) + (1 ∨ τ ) Σ 2 n .

B Sketch of the proof of Theorem 1

To deduce Theorem 1 from Theorem 2, the control of the entropy bracketing of S B (K,Jr ) for the Hellinger distance is required. For that, the proof of Maugis and Michel (2011b) is adapted for our specific mixtures. Next, a function ψ (K, J ) fulfilling Property (P) is deduced. We just give a sketch of the proof of Theorem 1 and we refer to [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] for more details.

B.1 Control of the entropy bracketing H

[.] (ε, S B (K,Jr) , d H )
To apply Theorem 2, the first step is to control the bracketing entropy of the Gaussian mixture families S B (K,Jr ) . Note that Theorem 2 only requires to control the local bracketing entropy H [.] (ε, S B (K,Jr ) (u, ξ), d H ). Nevertheless, it is difficult to characterize the subset S B (K,Jr ) (u, ξ) in function of the parameters of its mixtures. Thus, we rather control the global entropy bracketing H [.] (ε, S B (K,Jr ) , d H ), which is sufficient since the local bracketing entropy is upper bounded by the global bracketing entropy.

Inria

Proof. The key idea is that the control of the bracketing entropy of S B (K,Jr ) can be recast into the control of the bracketing entropies of the associated mixture component density families. Specifically, from (8), each mean vector µ k of a p-dimensional Gaussian mixture density in S B (K,Jr ) can be decomposed into a |J c r |-dimensional null mean vector and a |J r |-dimensional free mean vector:

S B (K,Jr ) =      K k=1 π k Φ • | µ k , σ 2 I ; ∀ k : µ k[J r ] ∈ [-A µ , A µ ] |Jr | , µ k[J c r ] = 0, ∀ k : π k > 0, K k=1 π k = 1, σ ∈ [a σ , A σ ]      . ( 47
)
Consider the (K -1)-dimensional simplex Π K defined by Π K := (π 1 , . . . , π K ) ∈ (0, 1) K ;

K k=1 π k = 1
and the family of K-tuples of p-dimensional Gaussian densities .

F (K,Jr ) =      Φ • | µ 1 , σ 2 I , . . . , Φ • | µ K , σ 2 I ; ∀ k : µ k[J r ] ∈ [-A µ , A µ ] |Jr | , µ k[J c r ] = 0, σ ∈ [a σ , A σ ]      .
From Lemma 2, all the matter is to calculate an upper bound of the bracketing entropy of F (K,Jr ) .

Let f = (f 1 , . . . , f K ) := Φ • | µ 1 , σ 2 I , . . . , Φ • | µ K , σ 2 I ∈ F (K,Jr ) . We want to find an ε-bracket for f . We shall consider shrunk and dilated Gaussian densities.

Step 1. Construction of a net for the variance Let δ ∈ (0, 1] to be chosen later. Let ⌈x⌉ denotes the smallest integer greater than or equal to x. We construct a regular net for the variance σ 2 ∈ [a 2 σ , A 2 σ ]. For l ∈ {2, . . . , r}, we define σ 2 l = (1 + δ) 1-l 2 A 2 σ where

r =     4 ln Aσ aσ √ 1 + δ ln(1 + δ)     (48)
is chosen so that σ 2 r < a 2 σ < σ 2 r-1 ≤ . . . ≤ σ 2 2 = A 2 σ .

cardinal of B ε (F (K,Jr ) ). From ( 48) and ( 50 .

From ( 48), (1+δ) (r-2)/4 ≤ (1+δ) 1/4 A σ /a σ ≤ 2 1/4 A σ /a σ and r -1 ≤ 4(A σ /a σ + 1/2)/δ, so

B ε F (K,Jr ) ≤ 4 2 5/4 A µ √ c ′ a σ K|Jr | A σ a σ + 1 2 δ -(1+K|Jr |) ≤ 4 2 5/4 A µ √ c ′ a σ K|Jr | A σ a σ + 1 2 √ cp ε 1+K|Jr | . (51)
Finally, Proposition 1 is derived from Lemma 2 and (51).

B.2 Determination of a function Ψ (K,Jr)

This section is devoted to the determination of a function Ψ (K,Jr ) defined by Property (P). From Proposition 1, for all ξ > 0, But from ( 46) and the fact that D (K,Jr ) = K(1 + |J r |), we have ln (C(A µ , A σ , a σ , K, J r , p)) we get that the function Ψ (K,Jr ) defined on R ⋆ + by Ψ (K,Jr ) (ξ) = D (K,Jr ) ξ B(A µ , A σ , a σ , p) + ln 1 ξ ∧ 1 satisfies (15). Besides, Ψ (K,Jr ) is nondecreasing and ξ → Ψ (K,Jr ) (ξ)/ξ is nonincreasing, so Ψ (K,Jr ) is convenient.

≤

B.3 Lower bound of the penalty function

Finally, according to the lower bound (18) of the penalty function, we need to find an upper bound of ξ * satisfying Ψ (K,Jr ) (ξ * ) = √ n ξ 2 * and to calculate the weights x (K,Jr ) to take into account the richness of the family S B (K,Jr ) . This can be done along the proofs of [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] by replacing the dimension of the models considered by [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] by the dimension D (K,Jr ) of our models S B (K,Jr ) . This leads to the two following lemmas:

Figure 2 :

 2 Figure 2: For one simulation with p = 1000, graphical validation of the penalty shape (pen(D) = 0.91D/n).

Figure 3 :

 3 Figure3: For one simulation with p = 200, ARI values for each model with K = 4 clusters obtained with the PS-Lasso and Lasso-MLE procedures are represented according to D/n. D/n is equal to 0.1 (resp. 0.59) for the model selected by PS-Lasso procedure with BIC (resp. by the Lasso oracle) and equal 0.14 for the Lasso-MLE procedure with DDSE and the oracle.

Inria

Figure 4 :

 4 Figure 4: Top left : true curves for the two clusters. Top right : curve estimations obtained by the Lasso-MLE procedure with a penalty without logarithm term. Bottom : curve estimations obtained by the alternative procedure with a penalty without (left) and with (right) a logarithm term.

  ŝ) -3.07(D/n) ln(K(1 + p)/D)

Figure 5 :

 5 Figure 5: For one simulation, slope graphs obtained by considering the model collection {S (K,Jr ,J0) } (K,Jr ,J0)∈M LL with a penalty without (top) and with ( bottom) a logarithm term. The resulting calibrated penalties are pen(D) = 2×1.75D/n and pen(D) = 2(4.42D/n+3.07(D/n) ln(K(1+p)/D)) respectively.

Proposition 1 .

 1 Put D (K,Jr ) = K(1 + |J r |). For all ε ∈ (0, 1], N [.] ε, S B (K,Jr ) , d H ≤ C(A µ , A σ , a σ , K, J r , p) µ , A σ , a σ , K, Jr, p) := 4 (2πe) c = sh(1) + 49/128 and c ′ = 5(1 -2 -1/4 )/8. Hence, H [.] ε, S B (K,Jr ) , d H ≤ ln [C(A µ , A σ , a σ , K, J r , p)] + D (K,Jr ) ln 1 ε .

Following

  the arguments developed by[START_REF] Maugis | Sélection de variables pour la classification non supervisée par mélanges gaussiens[END_REF] (proof of Proposition 7.A.2), it is easy to show that the study of the bracketing entropy of S B (K,Jr ) can be recast into the study of the bracketing entropy of Π K and F (K,Jr ) :Lemma 2. For all ε ∈ (0, 1],N [.] ε, S B (K,Jr ) , d H ≤ N [.] ε 3 , Π K , d H N [.] ε 3 , F (K,Jr ) , d H where N [.] (ε, Π K , d H ) ≤ K(2πe) K/2 1 ε K-1

ξ 0 H

 0 [.] (ε, S B (K,Jr ) , d H ) dε ≤ ξ ln (C(A µ , A σ , a σ , K, J r , p))In order to control the last term of the right-hand side of the last inequality, we apply the following technical result taken from[START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF]:Lemma 3. (Maugis and Michel, 2011b) For all ξ ∈ ]0, 1], Jr ) , d H ) dε ≤ ξ ln (C(A µ , A σ , a σ , K, J r , p)) + D (K,Jr ) (ξ ∧ 1) √ π + ln 1 ξ ∧ 1 ≤ D (K,Jr ) ξ √ π + ln (C(A µ , A σ , a σ , K, J r , p)) D (K,Jr ) + ln 1 ξ ∧ 1 .

Lemma 4 .

 4 Consider ξ * such that Ψ (K,Jr ) (ξ * ) = √ n ξ 2 * . Then, ξ 2 * ≤ D (K,Jr ) n 2B 2 (A µ , A σ , a σ , p) + ln 1 1 ∧ B 2 (A µ , A σ , a σ , p) D (K,Jr ) n .(52)Inria

  also confirms that variable selection is useful to get a better data clustering: introducing the relevant variables into the models improves the clustering but adding after the irrelevant variables deteriorates the clustering.

									Lasso
	0.88								Lasso + BIC
									Lasso oracle
	0.86								Lasso-MLE
									Lasso-MLE + DDSE or + oracle
	0.84
	0.82
	0.8
	0.78
	0.76
	0.74
	0.72
	0.7	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2

RR n°????

Acknowledgements

The authors are grateful to Gilles Celeux (INRIA Saclay-Île-de-France) and Pascal Massart (Université Paris-Sud) for helpful discussions and valuable comments.

Step 2. Construction of a net for the mean vectors Let l be the unique integer in {2, . . . , r} such that σ 2 l+1 < σ 2 ≤ σ 2 l . For all k ∈ {1, . . . , K}, let ν k ∈ R p to be specified later. Consider the functions defined on R p by

l I . Put l = (l 1 , . . . , l K ) and u = (u 1 , . . . , u K ). We now determine δ and (ν 1 , . . . , ν K ) so that l and u form an ε-bracket for f . On the one hand, by using the calculation of the Hellinger distance between two multivariate Gaussian densities (Maugis and Michel, 2011b, Corollary 3) and by upper bounding some usual functions, we get that, for all k ∈ {1, . . . , K}, d 2 H (l k , u k ) ≤ cp 2 δ 2 where c = sh(1) + 49/128. Thus, we take δ = ε/( √ cp) so that d H (l k , u k ) ≤ ε. On the other hand, by using the ratio of two multivariate Gaussian densities (Maugis and Michel, 2011b, Corollary 2), the definition of σ l and σ l+1 , the inequality ln(1 + δ) ≥ δ/2 for all δ ∈ (0, 1] and the concavity of δ → 1 -(1 + δ) -1/4 , we get that a sufficient condition for

where

For all k ∈ {1, . . . , K}, for all j ∈ J r , choose

Then, ν k (l) fulfills ( 49) and we get a net for the mean vectors.

Step 3. Upper bound of the number of ε-brackets for F (K,Jr )

From Step 1 and Step 2, the family

is an ε-bracket covering for F (K,Jr ) . Therefore, an upper bound of the number of ε-brackets necessary to cover F (K,Jr ) is deduced from an upper bound of the Inria Lemma 5. Consider the weight family x (K,Jr ) (K,Jr )∈M defined by

Then, we have (K,Jr )∈M e -x (K,Jr ) ≤ 1.

From ( 18), ( 52) and ( 53), we can apply Theorem 2 as soon as there exists κ > 0 such that pen (K, J r ) satisfies for all (K, J r ) ∈ M:

Applying Theorem 2 leads to Theorem 1.

C Grid of regularization parameters [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] propose to use a deterministic regular grid, similar for each number K of mixture components. In order to find the maximum value of regularization parameters, they conduct a few Lasso algorithms by increasing the value of λ until obtaining a model whose mean parameters all equal zero but this method can reveal quite time-consuming. Moreover, the choice of the grid step is difficult. To construct a grid of regularization parameters, another approach is here considered. The key idea is to construct a data-driven grid G K (depending on K) of regularization parameters by using the updating formulas of the mixture parameters estimation in the EM algorithm computing the Lasso solutions. More precisely, first the EM algorithm is used with λ = 0 to determine the ML estimator (π 0 1 , . . . , π0 K , μ0 1 , . . . , μ0 K , σ0 ). Then for all λ > 0 the parameter estimation at the r-th iteration of the EM algorithm is considered: for all k ∈ {1, . . . , K}, for all j ∈ {1, . . . , p}, π

where ν

ik and sign(.) is the sign function. From ( 54) and ( 55), we remark that r+1) .

Moreover, we hope that for all k ∈ {1, . . . , K} and for all j ∈ {1, . . . , p}, the values π λ kj in order to ensure a mixture with very sparse mean component vectors (the factor 2 is quite arbitrary; in practice the algorithm stops as soon as the null solution is reached, even is λ = 2 is not reached yet). In practice, this method is time-efficient since it only requires to run one EM algorithm for λ = 0.

Inria

RESEARCH CENTRE SACLAY -ÎLE-DE-FRANCE


Fotios Sidiroglou

Wilfried Blanc

Matt Kilburn

Thinh Nguyen

Ann Roberts

Bernard Dusardier

Gérard Monnom

Greg W Baxter

Improving the Radial Dopant Distribution in Silica Optical Fibres

Keywords: component, Fibre Frabrication, MCVD, Erbium Doped Optical Fibres

A novel approach for controlling and flattening the radial distribution of the erbium ions within an aluminium enriched silica core glass as part of the modified chemical vapour deposition and solution doping techniques is presented.

INTRODUCTION

Rare-earth (RE) doped optical fibres constitute a key component in the production of optical amplifiers, fibre lasers and sensor devices [START_REF] Digonnet | Rare-earth-doped fiber lasers and amplifiers[END_REF]. Although a number of fabrication techniques have been developed over the past few years (an extensive summary of which can be found in Ref [START_REF] Mendez | Specialty optical fibres handbook[END_REF]), the majority of fibres are still manufactured with the coupling of the modified chemical vapor deposition (MCVD) [START_REF] Macchesney | Preparation of low loss optical fibers using simultaneous vapor phase deposition and fusion[END_REF] and solution doping techniques [START_REF] Townsend | Solution doping technique for fabrication of rare earth doped optical fibers[END_REF]. This approach is easy to implement and offers great flexibility in doping REs with varying concentrations and in different hosts. However, control of dopant concentration and subsequently dopant distribution within the fabricated preform, which depends on a number of fabrication parameters (i.e., surface adsorption, solution concentration, soot density, etc), is still a major concern in this process. As a result a number of fabrication related issues (clustering, RE concentration depletion dip, etc) are still present, while the technology remains practically unchanged since its introduction.

During the solution doping process and as applied to MCVD, an unsintered layer of the core material is first obtained. The resulting porous tube is then soaked with an alcoholic/aqueous solution containing REs and other dopants, so that the dopants will impregnate into the pores of the core layer. The solution is removed from the tube and the impregnated layer is first dried before the tube is collapsed to form the fibre perform. Additional doped layers can be built up by repeating the aforementioned process prior to the collapse phase. However, after the solution doping step, the silica tube is always dried and collapsed to produce the preform.

In an alternative approach, we report on the fabrication of an erbium (Er 3+ ) doped optical fibre where additional silica layers were deposited following the sintering phase and prior to the preform collapse stage. Subsequent elemental characterisation has revealed a homogeneous dopant distribution profile that exhibits no signs of any dopant concentration depletion dip around the central core region. Information about the fibre fabrication and characterisation processes is presented in this paper.

II. FIBRE PREFORM FABRICATION

This experimental fibre was manufactured in an attempt to produce a "standards" fibre sample aimed to test the imaging properties of a near field scanning optical microscopic (NSOM) system. The latter had been modified to probe the Er 3+ ion distribution in optical fibres [START_REF] Sidiroglou | Micro-characterisation of RE doped optical fibres[END_REF]. The need for the existence of specific features that could test the ability of this imaging technique led to the fibre design shown in Fig. 1. It can be observed that the objective was to produce a fibre that consisted of two concentric rings containing Er 3+ , interleaved with silica layers. Using the MCVD and solution doping technique, the porous layers of the phosphorous (P) containing outer ring were deposited and later impregnated with an alcoholic solution of ErCl 3 .6H 2 O salts (0.01 mol/l in methanol) for 2 hours. The inner ring porous layers were impregnated with an alcoholic solution of ErCl 3 .6H 2 O (0.1 mol/l) and AlCl 3 .6H 2 O salts (1 mol/l) for a period of 2 hours as well. Four layers of pure silica were deposited to separate the outer and inner (core) ring layers, while two layers of pure silica were deposited after sintering the inner core region. The tube was then collapsed into a 10 mm-diameter solid preform, out of which the fibre was later drawn. [START_REF] Sidiroglou | Micro-characterisation of RE doped optical fibres[END_REF] displaying the relative ion distribution for all elements present in the produced fibre were also obtained from the application of a Nano-secondary ion mass spectroscopy (Nano-SIMS) scheme [START_REF] Sidiroglou | Simultaneous multidopant investigation of rare-earthdoped optical fibers by an ion microprobe[END_REF].

IV. RESULTS AND DISCUSSION

As evident from the set of data obtained as part of the characterisation process (Figs 2, 4, 5), this fibre exhibits a uniform elemental distribution in the central core region. In contrast to the conceptualised fibre design, the produced fibre appears to have a homogeneous RI profile with no signs of the interleaved silica regions that were deposited in order to form the two distinct Er 3+ containing rings of Fig 1. We believe that this attribute is due to the highly mobile nature of Al in silica, which continued to diffuse within the host glass matrix even after the sintering phase and the subsequent incorporation of the extra silica layers. In a recently published report, Tang et al studied the Al incorporation in sintered glass layers and collapsed preforms during a multiplesolution doping technique in optical fibre fabrication [START_REF] Tang | Multiple solution-doping in optical fibre fabrication I -Aluminium doping[END_REF][START_REF] Tang | Multiple solution-doping in optical fibre fabrication: II -Rare earth and aluminium co-doping[END_REF]. By measuring and comparing the peak Al content in the fabricated preform cores after sintering and then after the collapse stage, it was found that the peak Al-content in the preform core was on average about 15% less than for the sintered glass layer before preform collapse. They proposed that Al was redistributed during the collapse phase, since calculations showed that the overall Al-content in the preform core was only slightly less than in the overall volume of the sintered layer, indicating that Al was redistributed outwards from the centre of the core. Although the precise mechanisms for why is this happens are uncertain, a similar redistribution of the Al must have taken place during the fabrication of this fibre too. It is quite clear from the distribution profile that Al is redistributed not only outwards from the core centre (i.e., towards the P containing ring), but also inwards, indicating that Al mobility encouraged during the high temperatures of the extra silica layers deposition and the preform collapse. This is in agreement with previous reported data showing that alumina particles up to 80 µm in size can be dissolved completely in SiO 2 within 100 s at 1855 °C [START_REF] Risbud | Dependence of Phase Composition on Nuclei Available in SiO2-A12O3 Mixtures[END_REF].

More importantly, it can be observed that the RE ions have also diffused in a similar manner to that of Al. This observation is in agreement with past studies on RE doped silica fibres, where it has been empirically established that Al co-doping can be used for controlling the radial profiles of REs1. As a result, it becomes evident that the proposed fibre fabrication method promotes the redistribution of the Er3+ ions along the central core region. The latter is known to be integral for the successful operation of lasing and amplifying fibre devices.

Another clear benefit from the application of this fibre fabrication approach is related to the production of large mode area fibres. Manufacture of such fibres requires deposition of many supplementary doped layers fabricated with the repeat of the solution doping and sintering stages. This recurring process frequently leads to the formation of undesirable artefacts such as bubbles. With the technique proposed here, it is only necessary to create one additional dopant layer. Furthermore, high elemental concentration could be achieved with the implementation of the multiple solution doping technique [START_REF] Tang | Multiple solution-doping in optical fibre fabrication I -Aluminium doping[END_REF][START_REF] Tang | Multiple solution-doping in optical fibre fabrication: II -Rare earth and aluminium co-doping[END_REF]. Careful pre-selection of the number of additional silica layers to be deposited after the solution soaking and sintering stages can allow further control on the final core to cladding ratio.

Figure 1 .

 1 Figure 1. Fibre design

  PREFORM AND FIBRE CHARACTERISATIONA commercial refractive index (RI) profiler was employed to measure the RI difference (∆n) in the produced fibre (Fig2). Knowing the correspondence between index increase and oxide concentration in silica glass, the resulting ∆n was used to infer the final incorporated Al 2 O 3 and P 2 O 5 concentrations within the fibre. Note that the index increase due to Er 3+ was neglected due to its relatively low concentration and minor impact in raising the final RI profile around the Al doped region. The Al 2 O 3 and P 2 O 5 concentrations were estimated to be 8.0 mol% and 0.3 mol% respectively. Information about the Er 3+ ion concentration was deduced from the 1.5 µm absorption spectra (Fig 3) and in conjunction with the absorption cross section coefficients reported in Ref 2. Considering an attenuation coefficient of 230 dB/m and an absorption cross-section at 1.53 µm of 5×10 -25 m 2 , the Er 3+ ion concentration in the fibre was estimated about 2500 ppm mol Er 2 O 3 .

Figure 2 .

 2 Figure 2. Refractive index difference of Er 3+ doped fibre Additional information regarding the relative Er 3+ ion distribution in this fibre was obtained from the investigation of this fibre with the aid of a fluorescence intensity based confocal optical microscopy technique [6]. Two dimensional cross section images (Fig 4) as well as transverse linescans (Fig[START_REF] Sidiroglou | Micro-characterisation of RE doped optical fibres[END_REF] displaying the relative ion distribution for all elements present in the produced fibre were also obtained from the application of a Nano-secondary ion mass spectroscopy (Nano-SIMS) scheme[START_REF] Sidiroglou | Simultaneous multidopant investigation of rare-earthdoped optical fibers by an ion microprobe[END_REF].

Figure 3 .

 3 Figure 3. Spectral attenuation of Er 3+ doped fibre

Figure 4 .

 4 Figure 4. Cross section images of the relative elemental distribution obtained with the application of the NanoSIMS technique7. (a) Al, (b) P and (c) Er.

Figure 5 .

 5 Figure 5. Transverse profile through the centre of the fibre showing the relative Al (continuous line), P (dashed line) and Er (dotted line) ion distributions. Also shown is the RI difference (dashed-dotted line)

ACKNOWLEDGMENT

M. Ude and S. Trzesien from LPMC are greatly acknowledged for the preform and fibre fabrication.