Charles Bouveyron

Camille Brunet

Discriminative variable selection for clustering with the sparse Fisher-EM algorithm

Introduction

With the exponential growth of measurement capacities, the observed data are nowadays frequently high-dimensional and clustering such data remains a challenging problem. In particular, when considering the mixture model context, the corresponding clustering methods show a disappointing behavior in high-dimensional spaces. They suffer from the well-known curse of dimensionality [START_REF] Bellman | Dynamic Programming[END_REF] which is mainly due to the fact that model-based clustering methods are dramatically overparametrized in high-dimensional spaces. Moreover, even though we dispose of many variables to describe the studied phenomenon, most of the time, only a small subset of these original variables are in fact relevant.

Several recent works have been interested to simultaneously cluster data and reduce their dimensionality by selecting relevant variables for the clustering task. A common assumption to these works is that the true underlying clusters are assumed to differ only with respect to some of the original features. The clustering task aims therefore to group the data on a subset of relevant features. This presents two practical advantages: clustering results should be improved by the removing of non informative features and the interpretation of the obtained clusters should be eased by the meaning of retained variables. In the literature, variable selection for clustering is handled in two different ways.

On the one hand, some authors such as [START_REF] Law | Simultaneous Feature Selection and Clustering Using Mixture Models[END_REF][START_REF] Liu | Bayesian clustering with variable and transformation selection[END_REF][START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF][START_REF] Raftery | Variable selection for model-based clustering[END_REF] tackle the problem of variable selection for model-based clustering within a Bayesian framework. In particular, the determination of the role of each variable is recast as a model selection problem. A first framework was proposed by Raftery and Dean [START_REF] Raftery | Variable selection for model-based clustering[END_REF] in which two kinds of subsets of variables are defined: a subset of relevant variables and a subset of irrelevant variables which are independent from the clustering but which can be explained from the relevant variables through a linear regression. An extension of the previous work has then been proposed by Maugis et al. [START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF] who consider two kinds of irrelevant variables: the ones which can be explained by a linear regression from a subset of the clustering variables and finally a set of irrelevant variables which are totally independent of all the relevant variables. The models in competition are afterward compared with the integrated log-likelihood via a BIC approximation.

Even though these approaches present good results in most practical situations, their computational times are nevertheless very high and can lead to an intractable procedure in the case of high-dimensional data.

On the other hand, penalized clustering criteria have also been proposed to deal with the problem of variable selection in clustering. In the Gaussian mixture model context, several works, such as [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF][START_REF] Wang | Variable selection for model-based high dimensional clustering and its application to microarray data[END_REF][START_REF] Xie | Penalized model-based clustering with clusterspecific diagonal covariance matrices and grouped variables[END_REF][START_REF] Zhang | A flexible and efficient algorithm for regularized fisher discriminant analysis[END_REF] in particular, introduced a penalty term in the log-likelihood function in order to yield sparsity in the features. The penalty function can take different forms according to the constraints imposed on the structure of the covariance matrices. The introduction of a penalty term in the log-likelihood function was also used in the mixture of factor analyzers approaches, such as in [START_REF] Galimberti | Penalized factor mixture analysis for variable selection in clustered data[END_REF][START_REF] Xie | Penalized mixtures of factor analyzers with application to clustering high-dimensional microarray data[END_REF]. More recently, Witten and Tibshirani [START_REF] Witten | A framework for feature selection in clustering[END_REF] proposed a general non-probabilistic framework for variable selection in clustering, based on a general penalized criterion, which governs both variable selection and clustering. It appears nevertheless that the results of such procedures are usually not sparse enough and select a large number of the original variables, especially in the case of highdimensional data.

Other approaches focus on simultaneously clustering the data and reducing their dimensionality by feature extraction rather than feature selection. We can cite in particular, the subspace clustering methods [START_REF] Bouveyron | High Dimensional Discriminant Analysis[END_REF][START_REF] Ghahramani | The EM algorithm for factor analyzers[END_REF][START_REF] Mcnicholas | Parsimonious Gaussian mixture models[END_REF][START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF][START_REF] Montanari | Heteroscedastic Factor Mixture Analysis[END_REF][START_REF] Yoshida | A mixed factor model for dimension reduction and extraction of a group structure in gene expression data[END_REF] which are based on probabilistic frameworks and model each group in a specific and low-dimensional subspace. Even though these methods are very efficient in practice, they present nevertheless several limitations regarding the understanding and the interpretation of the clusters. Indeed, in most of subspace clustering approaches, each group is modeled in its specific subspace which makes difficult a global visualization of the clustered data. Even though some approaches [START_REF] Baek | Mixtures of Factor Analyzers with Common Factor Loadings: Applications to the Clustering and Visualisation of High-Dimensional Data[END_REF][START_REF] Montanari | Heteroscedastic Factor Mixture Analysis[END_REF] model the data in a common and low-dimensional subspace, they choose the projection matrix such as the variance of the projected data is maximum and this can not be sufficient to catch discriminative information about the group structure.

To overcome these limitations, Bouveyron and Brunet [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] recently proposed a new statistical framework which aims to simultaneously cluster the data and produce a low-dimensional and discriminative representation of the clustered data. The resulting clustering method, named the Fisher-EM algorithm, clusters the data into a common latent subspace of low dimensionality which best discriminates the groups according to the current fuzzy partition of the data. It is based on an EM procedure from which an additional step, named F-step, is introduced to estimate the projection matrix whose columns span the discriminative latent space. This projection matrix is estimated at each iteration by maximizing a constrained Fisher's criterion conditionally to the current soft partition of the data. As reported in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF], the Fisher-EM algorithm turned out to outperform most of the existing clustering methods while providing a useful visualization of the clustered data. However, the discriminative latent space is defined by "latent variables" which are linear combinations of the original variables. As a consequence, the interpretation of the resulting clusters according to the original variables is usually difficult. An intuitive way to avoid such a limitation would be to keep only large loadings variables, by thresholding for instance. Even though this approach is commonly used in practice, it has been particularly criticized by Cadima [START_REF] Cadima | Loadings and correlations in the interpretation of the principal components[END_REF] since it induces some misleading information. Furthermore, it often happens when dealing with high-dimensional data that a large number of noisy or non-informative variables are present in the set of the original variables. Since the latent variables are defined by a linear combination of the original ones, the noisy variables may remain in the loadings of the projection matrix and this may produce a deterioration of the clustering results.

To overcome these shortcomings, three different approaches are proposed in this work for introducing sparsity in the Fisher-EM algorithm and thus select the discriminative variables among the set of original variables. The remainder of this document is organized as follows. Section 2 reviews the discriminative latent mixture model of [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] and the Fisher-EM algorithm which was proposed for its inference. Section 3 develops three different procedures based on ℓ 1 penalties for introducing sparsity into the Fisher-EM algorithm. The first approach looks for the best sparse approximate of the solution of the F-step of the Fisher-EM algorithm. The second one recasts the optimization problem involved of the F-step as a lasso regression-type problem. The last approach is based on a penalized singular value decomposition (SVD) of the matrix involved in the constrained Fisher criterion of the F-step. Numerical experiments are then presented in Section 4 to highlight the practical behavior of the three sparse versions of the Fisher-EM algorithm and to compare them to existing approaches. In section 5, a sparse version of the Fisher-EM algorithm is applied to the segmentation of hyperspectral images. Section 6 finally provides some concluding remarks and ideas for further works.

The DLM model and the Fisher-EM algorithm

In this section, we briefly review the discriminative latent mixture (DLM) model [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] and its inference algorithm, named the Fisher-EM algorithm, which models and clusters the data into a common latent subspace. Conversely to similar approaches, such as [START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF][START_REF] Mcnicholas | Parsimonious Gaussian mixture models[END_REF][START_REF] Montanari | Dimensionally reduced mixtures of regression models[END_REF][START_REF] Montanari | Heteroscedastic Factor Mixture Analysis[END_REF][START_REF] Yoshida | A mixed factor model for dimension reduction and extraction of a group structure in gene expression data[END_REF], this latent subspace is assumed to be discriminative and its intrinsic dimension is strictly bounded by the number of groups.

The DLM model

Let {y 1 , . . . , y n } ∈ R p denote a dataset of n observations that one wants to cluster into K homogeneous groups, i.e. adjoin to each observation y i a value z i ∈ {1, . . . , K} where z i = k indicates that the observation y i belongs to the kth group.

On the one hand, let us assume that {y 1 , . . . , y n } are independent observed realizations of a random vector Y ∈ R p and that {z 1 , . . . , z n } are also independent realizations of a random variable Z ∈ {1, . . . , K}. On the other hand, let E ⊂ R p denote a latent space assumed to be the most discriminative subspace of dimension d ≤ K -1 such that 0 ∈ E and K < p. Moreover, let {x 1 , . . . , x n } ∈ E denote the actual data, described in the latent space E of dimension d, which are in addition

X Y Z π µ k Σ k W = [U, V ] ǫ Ψ Figure 1: Graphical summary of the DLM [Σ k β] model.
presumed to be independent realizations of an unobserved random vector X ∈ E. Finally, the observed variable Y ∈ R p and the latent variable X ∈ E are assumed to be linked through a linear transformation:

Y = UX + ε, ( 1
)
where U is a p × d orthonormal matrix common to the K groups and satisfying U t U = I d . The p-dimensional random vector ε stands for the noise term which models the non discriminative information and which is assumed to be distributed according to a centered Gaussian density function with a covariance matrix Ψ (ε ∼ N (0, Ψ)). Besides, within the latent space, X is assumed, conditionally to Z = k, to be Gaussian :

X |Z=k ∼ N (µ k , Σ k ) (2)
where µ k ∈ R d and Σ k ∈ R d×d are respectively the mean vector and the covariance matrix of the kth group. Given these distribution assumptions and according to equation [START_REF] Baek | Mixtures of factor analyzers with common factor loadings: applications to the clustering and visualisation of high-dimensional data[END_REF],

Y |X,Z=k ∼ N (UX, Ψ), (3)
and its marginal distribution is therefore a mixture of Gaussians:

f (y) = K k=1 π k φ(y; m k , S k ), (4)
∆ k =                   Σ k 0 0 β 0 . . . . . . 0 β                            d ≤ K -1                (p -d)
These last conditions imply that the discriminative and the non-discriminative subspaces are orthogonal, which suggests in practice that all the relevant clustering information remains in the latent subspace. This model is referred to by DLM

[Σ k β]
in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] and a graphical summary is given in Figure 1.

A family of parsimonious models

Several other models can be obtained from the DLM

DLM [Σ k β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K 2 (K -1)/2 + K DLM [Σ k β] (K -1) + K(K -1) + (K -1)(p -K/2) + K 2 (K -1)/2 + 1 DLM [Σβ k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K(K -1)/2 + K DLM [Σβ] (K -1) + K(K -1) + (K -1)(p -K/2) + K(K -1)/2 + 1 DLM [α kj β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K 2 DLM [α kj β] (K -1) + K(K -1) + (K -1)(p -K/2) + K(K -1) + 1 DLM [α k β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + 2K DLM [α k β] (K -1) + K(K -1) + (K -1)(p -K/2) + K + 1 DLM [αj β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + (K -1) + K DLM [αj β] (K -1) + K(K -1) + (K -1)(p -K/2) + (K -1) + 1 DLM [αβ k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K + 1 DLM [αβ] (K -1) + K(K -1) + (K -1)(p -K/2) + 2 Full-GMM (K -1) + Kp + Kp(p + 1)/2 Com-GMM (K -1) + Kp + p(p + 1)/2 Diag-GMM (K -1) + Kp + Kp Sphe-GMM (K -1) + Kp + K MFA (K -1) + Kp + Kd[p -(d -1)/2)] + Kp (d = 3) Mixt-PPCA (K -1) + Kp + K[d(p -(d + 1)/2) + d + 1] + 1 (d = 3) PGMM-CUU (K -1) + Kp + d[p -(d + 1)/2] + Kp (d = 3) MCFA (K -1) + Kd + p + d[p -(d + 1)/2] + Kd(d + 1)/2 (d = 3) MCUFSA (K -1) + Kd + 1 + d[p -(d + 1)/2] + Kd (d = 3)
Table 1: Number of free parameters to estimate when d = K -1 for the DLM models and some classical models (see text for details).

and d = 3, then the number of parameters to estimate for the DLM [Σ k β k ] is 337 which is drastically less than in the case of the Full-GMM (20 603 parameters to estimate). For a comparison purpose, Table 1 presents also the complexity of other clustering methods, such as Mixt-PPCA [START_REF] Tipping | Mixtures of Probabilistic Principal Component Analysers[END_REF], MFA [START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF], PGMM [START_REF] Mcnicholas | Parsimonious Gaussian mixture models[END_REF], MCFA [START_REF] Baek | Mixtures of factor analyzers with common factor loadings: applications to the clustering and visualisation of high-dimensional data[END_REF] and MCUFSA [START_REF] Yoshida | Array cluster: an analytic tool for clustering, data visualization and model finder on gene expression profiles[END_REF] for which the complexity grows linearly with p as well.

The Fisher-EM algorithm

An estimation procedure, called the Fisher-EM algorithm, is also proposed in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] in order to estimate both the discriminative space and the parameters of the mixture model. This algorithm is based on the EM algorithm from which an additional step is introduced, between the E and the M-step. This additional step, named F-step, aims to compute the projection matrix U whose columns span the discriminative latent space. The Fisher-EM algorithm has therefore the following form, at iteration q:

The E-step This step computes the posterior probabilities t (q)

ik that the observations belong to the K groups using the following update formula:

t (q) ik = π(q-1) k φ(y i , θ(q-1) k )/ K ℓ=1 π(q-1) ℓ φ(y i , θ(q-1) ℓ ), (5)
with θk = {μ k , Σk , βk , Û}.

The F-step This step estimates, conditionally to the posterior probabilities, the orientation matrix U (q) of the discriminative latent space by maximizing the Fisher's criterion [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF][START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF] under orthonormality constraints:

Û(q) = max U trace (U t SU) -1 U t S (q) B U , w.r.t. U t U = I d , ( 6
)
where S stands for the covariance matrix of the whole dataset and S (q)

B , defined as follows:

S (q) B = 1 n K k=1 n (q) k (m (q) k -ȳ)(m (q) k -ȳ) t , (7)
denotes the soft between covariance matrix with n

(q) k = n i=1 t (q) ik , m (q) k = 1/n (q) k n i=1 t (q)
ik y i and ȳ = 1/n n i=1 y i . This optimization problem is solved in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] using the concept of orthonormal discriminant vector developed by [START_REF] Foley | An optimal set of discriminant vectors[END_REF] through a Gram-Schmidt procedure. Such a process enables to fit a discriminative and low-dimensional subspace conditionally to the current soft partition of the data while providing orthonormal discriminative axes. In addition, according to the rank of the matrix S (q) B , the dimensionality of the discriminative space d is strictly bounded by the number of clusters K.

The M-step This third step estimates the parameters of the mixture model in the latent subspace by maximizing the conditional expectation of the complete loglikelihood:

Q(θ) = - 1 2 K k=1 n (q) k -2 log(π k ) + trace(Σ -1 k Û (q)t C (q) k Û (q) ) + log(|Σ k |) + (p-d) log(β k ) + trace(C (q) k )-d j=1 û(q)t j C (q) k û(q) j β k + p log(2π) . ( 8
)
where

C (q) k = 1 n (q) k n i=1 t (q) ik (y i -m (q) k )(y i -m (q)
k ) t is the empirical covariance matrix of the kth group and û(q) j is the jth column vector of Û(q) , n

(q) k = n i=1 t (q)
ik . Hence, maximizing Q conditionally to Û(q) leads to the following update formula for the mixture parameters of the model DLM

[Σ k β k ] : π(q) k = n (q) k n , (9) μ(q
) k = 1 n (q) k n i=1 t (q) ik Û(q)t y i , (10)
Σ(q) k = Û(q)t C k Û(q) , (11)
β(q) k = trace(C k )-d j=1 û(q)t j C k û(q) j p -d . ( 12
)
The Fisher-EM procedure iteratively updates the parameters until the Aitken criterion is satisfied (see paragraph 4.5 of [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF]). The convergence properties of the Fisher-EM algorithm have been studied in [START_REF] Bouveyron | Theoretical and practical considerations on the convergence properties of the Fisher-EM algorithm[END_REF]. It is also proposed in this work to use a stopping criterion based on the Fisher criterion involved in the F-step to improve the clustering performance. Finally, since the latent subspace has a low dimension and common to all groups, the clustered data can be easily visualized by projecting them into the estimated latent subspace.

Sparse versions of the Fisher-EM algorithm

Even though the Fisher-EM algorithm turns out to be efficient both for modeling and clustering data, the interpretation of clustering results regarding the original variables remains difficult. In this section, we propose therefore three different ways to introduce sparsity into the loadings of the projection matrix estimated in the F-step of the Fisher-EM algorithm.

A two-step approach

In this first approach, we propose to proceed in two steps. First, at iteration q, the traditional F-step of the Fisher-EM algorithm computes an estimate Û(q) of the orientation matrix of the discriminative latent space conditionally to the posterior probabilities t (q) ik . Then, the matrix Û(q) is approximated by a sparse one Ũ(q) using the following result.

Proposition 3.1. The best sparse approximation Ũ(q) of Û(q) at the level λ is the solution of the following penalized regression problem:

min U X (q)t -Y t U 2 F + λ d j=1 |U j | 1 , where U = [U 1 , ..., U d ], U j ∈ R p is the jth column vector of U, . F is the Frobenius norm and X (q) = Û(q)t Y .
Proof. Let Û(q) be the orientation matrix of the discriminative latent space estimated by the F-step at iteration (q) and let us define X (q) = Û(q)t Y ∈ R d×n the matrix of the projected data into the subspace spanned by Û(q) , where Y ∈ R p×n denotes the original data matrix. Since X (q) is generated by Û(q) , then Û(q) is solution of the least square regression of X (q) on Y :

min U X (q)t -Y t U 2 F , where U = [U 1 , ..., U d ], U j ∈ R p is the jth column vector of U, . F is the Frobenius norm.
A penalized version of this regression problem can be obtained by adding a ℓ 1 -penalty term as follows:

min U X (q)t -Y t U 2 F + λ d j=1 |U j | 1 ,
and the solution of this penalized regression problem is therefore the best sparse approximation of Û(q) at the level λ.

The previous result allows to provide a sparse approximation Ũ(q) of Û(q) but we have no guarantee that the Ũ(q) is orthogonal as required by the DLM model. The following proposition solves this issue. Proposition 3.2. The best orthogonal approximation of Ũ(q) is Ū(q) = u (q) v (q)t where u (q) and v (q) are respectively the left and right singular vectors of the SVD of Ũ(q) . Proof. Let us consider the matrix Ũ(q) . Searching the best orthogonal approximation of the matrix Ũ(q) is equivalent to solving the following optimization problem:

min U Ũ(q) -U 2 F w.r.t. U t U = I d .
This problem is a nearest orthogonal Procrustes problem which can be solved by a singular value decomposition [START_REF] Gower | Procrustes Problems[END_REF]. Let u (q) Λ (q) v (q)t be the singular value decomposition of Ũ(q) , then u (q) v (q)t is the best orthogonal approximation of Ũ(q) . From an practical point of view, the penalized regression problem of Proposition 3.1 can be solved by alternatively regressing each column vector of the projected matrix Û(q) . The sparse and orthogonal approximation Ū(q) of Ũ(q) is obtained afterward through a SVD of Ũ(q) . The following algorithm summarizes these steps.

Algorithm 1 -F-step of the sparseFEM-1 algorithm 1. At iteration q, compute the matrix Û(q) by solving (6).

2. Compute X (q) = Û(q)t Y .

3. For j ∈ {1, . . . , d}, solve d independent penalized regression problems with the LARS algorithm [START_REF] Efron | Least angle regression[END_REF]:

Ũ(q) j = arg min U j x (q)t j -Y t U j 2 + λ |U j | 1 ,
4. Repeat step 3 several times until convergence.

5. Let Ũ(q) = [ Ũ(q) 1 , ..., Ũ(q) d ]
, compute the SVD of Ũ(q) = u (q) Λ (q) v (q)t and let Ū(q) = u (q) v (q)t .

Let us remark that this problem can be extended to a more general penalized regression by adding a ridge penalty term. This allows in particular to handle the n < p case which occurs frequently nowadays. In such a case, the elastic-net algorithm [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] has to be used instead of the LARS algorithm in Algorithm 1.

Nevertheless, a limitation of such a procedure may be the disconnection between the estimation of the discriminative subspace and the introduction of the sparsity in the loadings of the projection matrix. To avoid that, the two following approaches aim to propose penalized Fisher criteria for which the solutions fit directly a sparse and discriminative latent subspace.

A penalized regression criterion

We therefore propose here to reformulate the constrained Fisher criterion (6) involved in the F-step of the Fisher-EM algorithm as a penalized regression problem.

Consequently, the solution of this penalized regression problem will fit directly a sparse and discriminative latent subspace. To this end, let us introduce the soft matrices H (q)

W and H (q) B which will be computed, conditionally to the E-step, at each iteration q of the sparse F-step as follows: Definition 3.1. The soft matrices H (q) W ∈ R p×n and H (q) B ∈ R p×K are defined, conditionally to the posterior probabilities t (q) ik computed in the E-step at iteration q, as follows:

H (q) W = 1 √ n Y - K k=1 t (q) 1k m (q) k , . . . , Y - K k=1 t (q) nk m (q) k ∈ R p×n (13) H (q) B = 1 √ n n (q) 1 (m (q) 1 -ȳ), . . . , n (q) K (m (q) K -ȳ) ∈ R p×K , ( 14
)
where n

(q) k = n i=1 t (q)
ik and m

(q) k = 1 n n i=1 t (q)
ik y i is the soft mean vector of the cluster k.

According to these definitions, the matrices

H (q)
W and H (q)

B satisfy:

H (q) W H (q)t W = S (q) W
and

H (q) B H (q)t B = S (q) B , ( 15
)
where

S (q) W = 1/n K k=1 n (q)
k C k stands for the soft within covariance matrix computed at iteration q and S (q) B denotes the soft between covariance matrix defined in equation [START_REF] Bouveyron | Theoretical and practical considerations on the convergence properties of the Fisher-EM algorithm[END_REF]. A penalized version of the optimization problem ( 6) can be therefore formulated as a penalized regression-type problem: Proposition 3.3. The best sparse approximation Ũ(q) of the solution of ( 6) at the level λ is the solution B(q) of the following penalized regression problem:

min A,B K k=1 R (q)-t W H (q) B,k -AB t H (q) B,k 2 F + ρ d j=1 β t j S (q) W β j + λ d j=1 |β j | 1 , w.r.t. A t A = I d , where A = [α 1 , . . . , α d ] ∈ R p×d , B = [β 1 , . . . , β d ] ∈ R p×d , R (q)
W ∈ R p×p is a upper triangular matrix resulting from the Cholesky decomposition of S (q)

W , i.e. S (q) W = R (q)t W R (q) W , H (q) B,k is the kth column of H (q)
B and ρ > 0.

Proof. First, let us consider that the matrix A is fixed at iteration q . Then, optimizing :

min A,B K k=1 R (q)-t W H (q) B,k -AB t H (q) B,k 2 F + ρ d j=1 β t j S (q) W β j ( 16
)
conditionally to A leads to consider the following regularized regression problem:

min B d j=1 H (q)t B R (q)-t W α j -H (q)t B β j 2 F + ρβ t j S (q) W β j , with B = [β 1 , . . . , β d ].
Solving this problem is equivalent to solving d independent ridge regression problem and the solution B(q) is :

B(q) = S (q) B + ρS (q) W -1 S (q) B R (q)-1 W A. ( 17
)
By substituting B(q) in Equation ( 16), optimizing the objective function ( 16) over A, given A t A = I d and B(q) fixed, is equivalent to maximize the quantity:

max A trace B(q)t H (q) B H (q)t B R (q)-1 W A , w.r.t. A t A = I d .
According to Lemma 1 of [START_REF] Qiao | Sparse linear discriminant analysis with applications to high dimensional low sample size data[END_REF], this is a Procrustes problem [START_REF] Gower | Procrustes Problems[END_REF] which has an analytical solution by computing the singular value decomposition of the quantity:

R (q)-t W (H (q) B H (q)t B ) B(q) = u (q) Λ (q) v (q)t ,
where the column vectors of the p × d matrix u (q) are orthogonal and v (q) is a d × d orthogonal matrix. The solution is Â(q) = u (q) v (q)t . Substituting Â(q) into (17) gives:

B(q) = R (q)-1 W R (q)-t W S (q) B R (q)-1 W + ρI p -1 R (q)-t W S (q) B R (q)-1 W Â(q) = R (q)-1 W u (q) Λ (q) + ρI p -1 Λ (q) v (q)t .
By remarking that the d eigenvectors associated to the non-zero eigenvalues of the generalized eigenvalue problem ( 6) are the columns of R (q)-1 W u (q) , it follows that B(q) spans the same linear subspace than the solution Û(q) of ( 6). Therefore, the solution of the penalized optimization problem: min

A,B K k=1 R (q)-t W H (q) B,k -AB t H (q) B,k 2 F + ρ d j=1 β t j S (q) W β j + λ d j=1 |β j | 1 , w.r.t. A t A = I d ,
is the best sparse approximation of the solution of ( 6) at the level λ.

However and as in the previous case, the orthogonality of the column vectors of B(q) is not guaranteed but this issue can be tackled by Proposition 3.2. From a practical point of view, the optimization problem of Proposition 3.3 can be solved using the algorithm proposed by [START_REF] Qiao | Sparse linear discriminant analysis with applications to high dimensional low sample size data[END_REF] in the supervised case by optimizing alternatively over B with A fixed and over A with B fixed. This leads to the following algorithm in our case:

Algorithm 2 -F-step of the sparseFEM-2 algorithm 1. At iteration q, compute the matrices H (q)

B and H (q)

W from Equations ( 13) and [START_REF] Foley | An optimal set of discriminant vectors[END_REF]. Let S (q)

W = H (q) W H (q)t W and S (q) B = H (q) B H (q)t B .

Compute R (q)

W by using a Cholesky decomposition of S (q) W + γ/p trace(S

(q) W ) = R (q)t W R (q) W .

Initialization:

Let B (q) be the eigenvectors of

S -1 S (q) B . Compute the SVD R (q)-t W S (q)
B B (q) = u (q) Λ (q) v (q)t and let A (q) = u (q) v (q)t . 4. Solve d independent penalized regression problems. For j = 1, . . . , d:

β(q) j = arg min β j β t j W (q)t W (q) β j -2 Ỹ (q)t W (q) β j + λ 1 β j 1 ,
where

W (q) = H (q)t B √ ρR (q) W and Ỹ (q) = H (q)t B R (q)-1 W α (q) j O p . 5. Let B(q) = [ β1 , . . . , βd ]. Compute the SVD of R (q)-t W S (q) B
B(q) = u (q) Λ (q) v (q)t and let A (q) = u (q) v (q)t . 6. Compute the SVD of B(q) = u ′(q) Λ ′(q) v ′(q)t and let Ū(q) = u ′(q) v ′(q)t . 7. Repeat steps several times until convergence.

A penalized singular value decomposition

In this last approach, we reformulate the constrained Fisher criterion (6) involved in the F-step of the Fisher-EM algorithm as a regression problem which can be solved by doing the SVD of the matrix of interest in this regression problem. A sparse approximation of the solution of this regression problem will be obtained by doing a penalized SVD [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF] instead of the SVD. To that end, let us consider the following result.

Proposition 3.4. The solution of (6) is also solution of the following constrained optimization problem:

min U p ℓ=1 S (q) B,ℓ -UU t S (q) B,ℓ 2 w.r.t. U t U = I d ,
where S (q) B,ℓ is the ℓth column of the soft between covariance matrix S (q) B computed at iteration q.

Proof. Let us first prove that minimizing the quantity p ℓ=1 ||S

(q) B,ℓ -UU t S (q) B,ℓ || 2 is equivalent to maximize trace(U t S (q) B S (q)t B U).
To that end, we can write down the following equalities:

p ℓ=1 S (q) B,ℓ -UU t S (q) B,ℓ 2 = p ℓ=1 trace S (q)t B,ℓ (I p -UU t ) t (I p -UU t )S (q) B,ℓ = trace (I p -UU t ) t (I p -UU t ) p ℓ=1 S (q) B,ℓ S (q)t B,ℓ = trace S (q)t B (I p -UU t )S (q) B = trace(S (q)t B S (q) B ) -trace(U t S (q) B S (q)t B U).
Consequently, minimizing over U the quantity p ℓ=1 ||S

(q) B,ℓ -UU t S (q) B,ℓ || 2 is equivalent to maximize trace(U t S (q) B S (q)t
trace U t S B S t B U = trace U t (uΛ d v t )(uΛ d v t ) t U , = trace U t uΛ d Λ t d u t U , = d j=1 λ 2 j .
Consequently, the p × d orthogonal matrix Û such that p ℓ=1 ||S

(q) B,ℓ -UU t S (q) B,ℓ || 2 is minimized, is the matrix made of the d first left eigenvectors of S (q) B . Besides, since S (q)
B is symmetric and semi-definite positive, the matrix Û contains also the eigenvectors associated with the d largest eigenvalues of S (q)2 B and therefore the ones of S (q) B . Therefore, assuming without loss of generality that S = I p , Û is also solution of the constrained optimization problem (6) involved in the original F-step.

The optimization problem of Proposition (3.4) can be seen as looking for the projection matrix U such that the back-projection UU t S (q) B,ℓ is as close as possible to S (q) B,ℓ . In [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF], Witten et al. have considered such a problem with a constraint of sparsity on U. To solve this problem, they proposed an algorithm which performs a penalized SVD of the matrix of interest in the constrained optimization problem. Therefore, it is possible to obtain a sparse approximation Ũ(q) of the solution of ( 6) by doing a penalized SVD of S (q) B with the algorithm of [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF]. As previously, the orthogonality of the column vectors of Ũ(q) is not guaranteed but this issue can be again tackled by Proposition 3.2. From a practical point of view, this third approach can be implemented as follows:

Algorithm 3 -F-step of the sparseFEM-3 algorithm 1. Let M 1 = S (q)
B and d = rank (S B ) .

2. For j ∈ {1, . . . , d}:

(a) Solve û(q) j = arg max u j u t j M j v j w.r.t. u j 2 2 ≤ 1, v j 2 2 ≤ 1 and p ℓ=1 u (q)
jℓ ≤ λ 1 using the penalized SVD algorithm of [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF].

(b) Update M j+1 = M j -λ j u (q) j v t j . 3. Û(q) = [û (q) 1 , . . . , û(q) d ].
4. Compute the SVD of Û(q) = u (q) Λ (q) v (q)t and let Ũ(q) = u (q) v (q)t .

Practical aspects

The introduction of sparsity in the Fisher-EM algorithm presents several practical aspects among which the ability to interpret the discriminative axes. However, two questions remain: the choice of the hyper-parameter which determines the level of sparsity and the implementation strategy in the Fisher-EM algorithm. Both aspects are discussed below.

Choice of the tuning parameter

The choice of the threshold λ is an important problem since the number of zeros in the d discriminative axes depends directly on the degree of sparsity. In [START_REF] Zou | Sparse Principal Component Analysis[END_REF], Zou et al. chose the hyper-parameter of their sparse PCA with a criterion based on the explanation of the variance approximated by the sparse principal components. In [START_REF] Witten | A framework for feature selection in clustering[END_REF], Witten and Tibshirani proposed for their sparse-kmeans to base the choice of the tuning parameter on a permutation method closely related to the gap statistic previously proposed by Tibshirani et al. [START_REF] Tibshirani | Estimating the number of clusters in a dataset via the gap statistic[END_REF] for estimating the number of components in standard kmeans. Since our model is defined in a Gaussian mixture context, we propose to use the BIC criterion to select the threshold λ. According to the consistency results obtained by Zou et al. [START_REF] Zou | On the degrees of freedom of the Lasso[END_REF] and the fact that the sparsity constraint is applied on the projection matrix U, the effective number of parameters to estimate in the DLM [Σ k β k ] model is:

γ e = (K -1) + Kd + (d[p-(d + 1)/2] -d e ) + Kd(d + 1)/2 + K
where d e is the number of zeros in the loading matrix. In the same manner, this effective number of parameters to estimate can be declined for the 11 other submodels of the DLM family.

Implementation of the sparse Fisher-EM algorithm

We identified two different ways to implement the sparse versions of the Fisher-EM algorithm. First, it could be possible to replace the usual F-step of the Fisher-EM algorithm by a sparse F-step developed previously. The resulting algorithm would sparsify at each iteration the projection matrix U before estimating the model parameters. This can however leads to some drawbacks since an early introduction of the ℓ 1 penalty could penalize too much the loadings of the projection matrix, in particular if the initialization is far away from the optimal situation. An alternative implementation would be to, first, execute the traditional Fisher-EM algorithm until convergence and, then, initialize the sparse Fisher-EM algorithm with the result of the Fisher-EM algorithm. This strategy should combine the efficiency of the standard Fisher-EM algorithm with the advantage of having a sparse selection of discriminative variables. We therefore recommend this second implementation and it will be used in the experiments presented in the following sections.

Experimental comparison

This section presents comparisons with existing variable selection techniques on simulated and real-world data sets.

Comparison on simulated data

This first experiment aims to compare on simulated data the performances of the proposed sparseFEM algorithms (sparseFEM-1, sparseFEM-2, sparseFEM-3) to several competitors: Clustvarsel of Raftery and Dean [START_REF] Raftery | Variable selection for model-based clustering[END_REF], Selvarclust of Maugis et al. [START_REF] Maugis | Variable selection in model-based clustering: A general variable role modeling[END_REF] and sparse-kmeans of Witten and Tibshirani [START_REF] Witten | A framework for feature selection in clustering[END_REF]. For this experiment, we replicated the simulation proposed in Section 3.3 of [START_REF] Witten | A framework for feature selection in clustering[END_REF]. We simulated K = 3

Gaussian components of n observations in a 25-dimensional observation space whose components differ only on q = 5 features. The used parameters were Table 2: Clustering errors and numbers of non-zero variables averaged over 20 simulations for several clustering methods with p = 25 and q = 5. The results of Selvarclust are reported from [START_REF] Celeux | Letter to the editor[END_REF].

µ kj = µ × (1
ponents and σ 2 kj = 1 for the variance terms. Moreover, four different situations are considered: n = 30 or 300 and µ = 0.6 or 1.7. Each simulation was replicated 25 times.

Table 2 presents the means and standard deviations for both the clustering error and the number of non-zero variables for kmeans, sparse-kmeans, Clustvarsel, Selvarclust and the 3 procedures of sparseFEM. Note that the results of Selvarclust corresponds to clustering errors and non-zero variable rates found in [START_REF] Celeux | Letter to the editor[END_REF]. Moreover, the reported results concerning the 3 sparse Fisher-EM algorithms were obtained with the DLM [α k β] model for a sparsity level corresponding to the highest BIC value obtained at each trial.

Two main remarks can be done on the results presented in Table 2. First, by considering either the most difficult clustering cases (n = 30 and µ = 0.6) or the easiest one (n = 300 and µ = 0.6 or 1.7), all approaches present approximatively the same results in terms of clustering error rate. The methods differ however in the number of variables they retain to perform the clustering: Clustevarsel, sparseFEM-1, sparseFEM-2 and sparseFEM-3 turn out to select significantly less variables than sparse-kmeans and Clustvarsel. In particular, Clustevarsel and the sparseFEM algorithms select a number of useful variables consistent with the actual number of meaningful variables (q = 5). Second, in the situation where n = 30 and µ = 1.7, Selvarclust and sparse-kmeans present the lowest misclassification rate (0.08), even though the clustering error of sparseFEM-1 and kmeans remains relatively low (0.14). However, as previously, only Clustevarsel and the sparseFEM algorithms select a number of variables close to the right number of discriminative features.

Comparison on real data sets

Real-world data sets are now used to compare the efficiency of the sparseFEM algorithms to its competitors for both the clustering and variable selection tasks.

We considered 7 different benchmark data sets coming mostly from the UCI machine learning repository. We selected these data sets because they represent a wide range of situations in term of number of observations n, number of variables p and number of groups K. These characteristics are given in the top row of Table 3 and a detailed description of these data sets can be found in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF].

For this experiment, we used the 3 sparseFEM algorithms and the 3 sparse methods introduced previously (sparse-kmeans, Clustvarsel and Selvarclust). Since the evaluation of the clustering performance is a complex and very discussed problem, we chose to evaluate the clustering performance as the adequacy between the re-sulting partition of the data and the known labels for these data. For each data set, the sparseFEM algorithms were initialized with a common random partition drown from a multinomial distribution with equal prior probabilities. For Clustvarsel, Selvarclust and sparse-kmeans, the initialization was done with their own deterministic procedure. Moreover, for each method, the number K of groups has been fixed to the actual one. For Clustvarsel, Selvarclust and sparse-kmeans, the determination of the other free parameters was done according to the tools provided by each approach. For the sparseFEM algorithms, we used the penalized BIC criterion to select the model and the level of sparsity. More precisely, we first chose the model presenting the highest average BIC value on 20 replications. Then, given the selected model, we selected the level of sparsity associated with the highest BIC value.

Table 3 presents the average clustering accuracies and the associated standard deviations obtained for the 6 approaches. The average number of non-zero variables is also reported within brackets in the table. The results associated to the sparseFEM algorithms have been obtained by averaging over 20 trials with random initializations. The lack of standard deviations for Clustvarsel, Selvarclust and sparse-kmeans is due to the deterministic initializations they use. It first appears that the three sparse versions of the Fisher-EM algorithm perform rather similarly both in term of clustering and variable selection. It also appears clearly that the sparseFEM algorithms are competitive to existing methods regarding both the clustering performances and the selection of variables. Indeed, the sparseFEM algorithms obtain the best clustering accuracies on 4 of the 7 data sets whereas sparse-kmeans and Selvarclust obtain the best clustering accuracies on respectively 2 and 1 data sets.

The sparseFEM algorithms differ also from sparse-kmeans regarding the number of variables retained to perform the clustering. Indeed, sparse-kmeans turns out to frequently select a large number of variables whereas sparseFEM is usually rather sparse in the number of selected variables. Finally, Clustvarsel and Selvarclust turn out to select most of the time few variables, particularly in high-dimensional spaces, which seems to obstruct their clustering performance. To summarize, this experiment has shown that the sparseFEM algorithms seem to be good compromises between sparse-kmeans and Clustvarsel /Selvarclust in term of variable selection and, certainly thanks to this characteristic, they also provide good clustering results.

Comparison on the usps358 data set

We focus now on the usps358 dataset to stress the role of variable selection in the interpretation of clustering results. divided in 10 classes corresponding to the digits from 0 to 9. Each digit is a 16 × 16 gray level image represented as a 256-dimensional vector. For this experiment, we extracted a subset of the data (n = 1 756) corresponding to the digits 3, 5 and 8 which are the three most difficult digits to discriminate. This smaller dataset is hereafter called usps358. Figure 2 depicts the group mean images obtained from the true labels in the usps358 dataset. For this experiment, we used the three sparse-FEM algorithms with the model and the level of sparsity selected in the previous experiment for this data set. For Clustvarsel, Selvarclust and sparse-kmeans, the level of sparsity was again selected with their own selection procedure. First of all, it appears that Clustvarsel and Selvarclust select significantly fewer variables than both sparse-kmeans or the sparseFEM procedures. Furthermore, most of the selected variables by Clustvarsel and Selvarclust turn out to be irrelevant to discriminate the digit 3 from the digits 5 and 8. For instance, in Figures 3.b and 3.c, we can observe that the black pixels located in right bottom corner, do not correspond to any discriminative variable. This certainly explain the poor clustering performances (48.3% for Clustvarsel and 36.7% for Selvarclust) observed on this data set for these methods. On the contrary, sparse-kmeans turns out to perform well in term of clustering performance (74.7% of clustering accuracy). Nevertheless, the number of selected variables remains higher (213 selected variables amongst 256 original ones) than we would expect to ease the interpretation of results. Finally, sparseFEM-1 and sparseFEM-2 seem to answer quite well to both the clustering task and the task of feature selection. Indeed, on the one hand, the subset of selected pixels remains small for both algorithms: 6 and 15 pixels are selected amongst 256 for sparseFEM-1 and sparseFEM-2 respectively. Furthermore, the selected pixels appear to be relevant to discriminate the classes associated with the three digits. For instance, the darker pixel on the bottom right corner of Figure 4.b discriminates the digit 8 from the digits 3 and 5. On the other hand, and certainly due to this relevant selection of variables, both algorithms perform particularly well on this high-dimensional data set (84.7% for sparseFEM-1 and 82.8% for sparseFEM-2). However, on this data set, the sparseFEM-3 procedure shows a disappointing behavior regarding the variable selection even though its clustering performance remains satisfying. The fact that sparseFEM-3 succeeds in clustering the data set even with a bad selection of variables is certainly due to the nature of the DLM model which models also the non discriminative information through the parameter β k . Table 4 presents the computing time of the studied clustering methods (for a given model and with λ and K fixed) for clustering the usps358 data set. As we can remark, our procedures are much faster than the sparse-kmeans, Clustvarsel and Selvarclust algorithms. Consequently, the sparseFEM algorithms appear once again to be good compromises, in practice, to cluster high-dimensional data and select a set of discriminative variables in a reasonable time.

Application to the segmentation of hyperspectral images

Here, we propose to use sparseFEM to segment hyperspectral images of the Martian surface. Visible and near infrared imaging spectroscopy is a key remote sensing technique to study the system of the planets. Imaging spectrometers, which are inboard of an increasing number of satellites, provide high-dimensional hyperspectral images. In March 2004, the OMEGA instrument (Mars Express, ESA) [START_REF]Bibring and 42 co-authors. Mars Surface Diversity as Revealed by the OMEGA/Mars Express Observations[END_REF] has collected 310 Gbytes of raw images. The OMEGA imaging spectrometer has mapped the Martian surface with a spatial resolution varying between 300 to 3000 meters depending on the spacecraft altitude. It acquired for each resolved pixel the spectrum from 0.36 to 5.2 µm in 256 contiguous spectral channels. OMEGA is designed variables (wavelengths) as discriminative variables, i.e. the rows associated to these variables were non-zero in the loading matrix U. Looking closely at the selection, we indeed notice that the first selected variable (from left to right) discriminates the blue group from the others. The second selected variable discriminates the red and green groups from the black, blue and light blue groups whereas the third selected variable allows to discriminate the red, green and black groups from the blue and light blue groups. Similarly, the fourth and fifth selected variables discriminate the red and green groups from the black, blue and light blue groups whereas the sixth, seventh and eighth selected variable allows to discriminate the red, green and light blue groups from the blue and black groups.

A possible interest of such a selection could be the measurement of only a tens of wavelengths for future acquisitions instead of the 256 current ones for a result expected to be similar. This could in particular reduce the acquisition time for each pixel from a few tens of seconds to less than one second.

Conclusion

This article has focused on variable selection for clustering with the Fisher-EM algorithm which has been recently proposed in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF]. The aim of this work was to introduce sparsity in the Fisher-EM algorithm and thus select the discriminative variables among the set of original variables. We have proposed three different procedures based on a ℓ 1 -penalty term. Experiments on simulations and real data sets have shown that the three sparse versions of the Fisher-EM algorithm are highly competitive with existing approaches of the literature. In particular, the sparseFEM procedures present several assets regarding existing approaches. On the one hand, they tend to select an intermediate number of discriminative variables whereas existing approaches tend to select either too few (Clustvarsel and Selvarclust) or too much variables (sparse-kmeans). On the other hand, the sparseFEM procedures perform both the clustering and the variable selection in a reasonable time comparing to existing approaches in the case of high-dimensional data. The sparseFEM algorithms have been also applied with success to the segmentation of hyperspectral images of the planet Mars and relevant parts of the spectra which well discriminate the groups have been identified.

Among the possible extensions of this work, it may be first interesting to use different ℓ 1 -penalty values according to the relevance of each discriminative axis estimated in the Fisher-EM algorithm. Such an approach could identify different levels of relevancy among the original variables. Second, we used in this work a penalized BIC criterion to select the sparsity level by evaluating the model complexity in regard to the non-zero values as proposed by [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF]. Although Zou et al. [START_REF] Zou | On the degrees of freedom of the Lasso[END_REF] showed that the number of non-zero coefficients is an unbiased estimate of the degrees of freedom and is asymptotically consistent in the case of penalized regression problem, this result has no theoretical justification in the penalized GMM context. It would be therefore interesting to obtain theoretical guarantees of such a result in our context. Finally, since the ICL criterion [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] is also used to select the number of components, it would be a natural extension to consider a penalized ICL for selecting the sparsity level in the sparseFEM algorithms.

B

  U) according to U. Let us now consider the SVD of the n × p matrix S (q) B = uΛv t where u and v stands for respectively the left and right singular vectors of S (q) B and Λ is a diagonal matrix containing its associated singular values. Since the matrix S (q) B has a rank d at most equal to K -1 < p, with K the number of clusters, then only d singular values of the matrix S (q) B are non zeros, which enables us to write S (q) B = uΛ d v t , where Λ d = diag(λ 1 , . . . , λ d , 0, . . . , 0). Moreover, by letting U = u d the d first left eigenvectors of S B , then:

Figures 3

 3 Figures 3 illustrates, as images, the features selected respectively by sparsekmeans (Figure 3.a), Clustvarsel (Figure 3.b) and Selvarclust (Figure 3.c). In Figure 3.a, the weight assigned by sparse-kmeans to each feature is represented by gray levels: lighter is the pixel, weaker is the absolute value of the weight of the associated feature. For Clustvarsel and Selvarclust, only the selected variables are depicted and are associated to black pixels as it is illustrated in Figures 3.b and 3.c respectively. These representations are associated to the following clustering accuracies 74.7%, 48.3% and 36.7% for sparse-kmeans, Clustvarsel and Selvarclust respectively. For the 3 sparseFEM algorithms, we superimposed in a same figure the absolute values of the loadings of the two discriminative axes fitted by the sparseFEM-1, sparseFEM-2 and sparseFEM-3 procedures. The associated clustering accuracies are respectively 84.7%, 82.8% and 79.1%.

Figure 2 :

 2 Figure 2: Group means obtained from the true labels in the USP358 datasets.

Figure 3 :

 3 Figure 3: Variable selection obtained from (a) the sparse-kmeans algorithm, (b) the Clustvarsel approach and (c) the Selvarclust approach.

3 Figure 4 :

 34 Figure 4: Variable selection obtained from (a) the sparseFEM-1, (b) the sparseFEM-2 and (c) the sparseFEM-3 procedures with sparsity levels selected by the penalized BIC.

Figure 5 :

 5 Figure 5: Image of the studied zone of the Martian surface.

Figure 6 :

 6 Figure 6: Some of the 38 400 measured spectra described on 256 wavelengths (see text for details).

Figure 8 :

 8 Figure 8: Mean spectra of the 5 groups formed by sparseFEM and selection of the discriminative wavelengths (indicated by gray rectangles).

  where π k is the mixing proportion of the kth group and φ(.; m k , S k ) denotes the multivariate Gaussian density function parametrized by the mean vector m k = Uµ k and the covariance matrix S k = UΣ k U t +Ψ of the kth group. Furthermore, we define the p × p matrix W = [U, V ] such that W t W = W W t = I p , where the (pd) × p matrix V is an orthogonal complement of U. Finally, the noise covariance matrix Ψ is assumed to satisfy the conditions V ΨV t = βI p-d and UΨU t = 0 d , such that ∆ k = W t S k W has the following form:

  In particular, the complexity of the DLM [Σ k β k ] model mainly depends on the number of clusters K since the dimensionality of the discriminative subspace is such that d ≤ K -1. Notice that the complexity of the DLM [Σ k β k ] grows linearly with p contrary to the traditional Gaussian models in which the complexity increases with p 2 . As an illustration, if we consider the case where p = 100, K = 4

	[Σ k β] model by relaxing or adding constraints on model parameters. Firstly, it is possible to consider a more general case than the DLM [Σ k β] by relaxing the constraint on the variance term of the non discriminative information. Assuming that ε |Z=k ∼ N (0, Ψ k ) yields the DLM [Σ k β k ] model which can be useful in some practical cases. From this extended model, 10 parsimonious models can be obtained by constraining the parameters Σ k and β k to be common between and within the groups. For instance, the covariance matrices Σ 1 , . . . , Σ K in the latent space can be assumed to be common across the groups and this sub-model is referred to by DLM [Σβ k ] . Similarly, in each group, Σ k can be assumed to be diagonal, i.e. Σ k = diag(α k1 , . . . , α kd ). This sub-model is referred to by DLM [α kj β k ] . These sub-models can also be declined by considering that the parameter β is common to all classes (∀k, β k = β). A list of the 12 different DLM models is given by Table 1 and detailed descriptions can be found in [6]. Such a family yields very parsimonious models and allows, in the same time, to fit into various situations. Model K = 4 and Nb. of parameters p = 100

  The original dataset is made of 7 291 images

		iris	wine	chiro	zoo	glass	satimage	usps358
		(p=4,K=3) (p=13,K=3) (p=17,K=3) (p=16,K=7) (p=9,K=7) (p=36,K=6) (p=256,K=3)
	Approaches	(n=150)	(n=178)	(n=178)	(n=101)	(n=214)	(n=4435)	(n=1726)
	sparseFEM-1	96.5±0.3	97.8±0.2	84.2±11	71.4±8.5	50.2±1.9	69.6±0.1	84.7±3.2
		(2.0±0.0)	(2.0±0.0)	(2.3±0.5)	(13±2.5)	(6.0±1.0)	(36±0.0)	(5.5±0.7)
	sparseFEM-2	89.9±0.4	98.3±0.0	84.8±12	70.1±12.2	48.4±3.0	67.5±1.6	82.8±9.1
		(4.0±0.0)	(4.0±0.0)	(2.0±0.6)	(14±3.6)	(6.6±0.7)	(36±.0.0)	(15.5±16)
	sparseFEM-3	96.5±0.3	97.8±0.0	82.9±12	72.0±4.3	48.2±2.7	71.8±2.3	79.1±7.4
	sparse-kmeans	(2.0±0.3) 90.7	(2.0±0.0) 94.9	(2.0±0.0) 95.3	(10±2.8) 79.2	(7.0 ±0.0) 52.3	(36±0.0) 71.4	(6.0±1.3) 74.7
		(4.0)	(13.0)	(17.0)	(16.0)	(6.0)	(36.0)	(213)
	Clustvarsel	96.0	92.7	71.1	75.2	48.6	58.7	48.3
		(3.0)	(5.0)	(6.0)	(3.0)	(3.0)	(19.0)	(6.0)
	Selvarclust	96.0	94.4	92.6	92.1	43.0	56.4	36.7
		(3.0)	(5.0)	(8.0)	(5.0)	(6.0)	(22.0)	(5.0)

Table 3 :

 3 Clustering accuracies and their standard deviations (in percentage) on 7 UCI datasets (iris, wine, chironomus, zoo, glass, satimage, usps358) averaged on 20 trials. The average number of nonzero variables is reported in brackets. No standard deviation is reported for Clustvarsel/Selvarclust and sparse-kmeans since their initialization procedure is deterministic and always provides the same initial partition.

Table 4 :

 4 Computing times for the 3 versions of the sparseFEM algorithm, sparsekmeans, Clustvarsel and Selvarclust on the USPS358 data (for a given model and with λ and K fixed).

	Approaches: Procedure time (sec) Approaches:	Procedure time (sec)
	sparseFEM-1	729.04	sparse-kmeans	1 567.75
	sparseFEM-2	387.12	Clustvarsel	2 957.70
	sparseFEM-3	409.61	Selvarclust	9 257.10

Acknowledgments

The authors would like to thank Cathy Maugis for providing the results of Selvarclust on the zoo, glass, satimage and usps358 data sets.

Expert segmentation

to characterize the composition of surface materials, discriminating between various classes of silicates, hydrated minerals, oxides and carbonates, organic frosts and ices. For this experiment, a 300 × 128 image of the Martian surface is considered and a 256-dimensional spectral observation is therefore associated to each of the 38 400 pixels. Figure 5 presents an image of the studied zone and Figure 6 shows some of the 38 400 measured spectra. According to the experts, there are K = 5 mineralogical classes to identify.

The sparseFEM-1 algorithm was applied to this dataset using the model DLM [α kj β] and a sparsity ratio equals to 0.1 (it refers to the ratio of the ℓ 1 norm of the coefficient vector relative to the norm at the full least square solution). The sparseFEM algorithm was initialized with the results of the Fisher-EM algorithm and the whole segmentation process took 18 hours on a 2.6 Ghz computer. Figure 7 presents, on the right panel, the segmentation into 5 mineralogical classes of the studied zone with the sparseFEM algorithm. In comparison, the left panel of Figure 7 shows the segmentation obtained by experts of the domain using a physical model. It first appears that the two segmentations agree globally on the mineralogical nature of the surface of the studied zone (60.30% of agreement). We recall that both segmentations do not exploit the spatial information. When looking at the top-right quarter of the image, we can notice that sparseFEM seems to provide a finer segmentation than the segmentation based on the physical model. Indeed, sparseFEM segments better than the physical model the fine "rivers" which can be seen on Figure 5.

Finally, Figure 8 shows the mean spectra of the 5 groups formed by sparseFEM and the selection of the discriminative wavelengths. SparseFEM selected 8 original


Junjie Lai

André Seznec

Bound the Peak Performance of SGEMM on GPU with software-controlled fast memory

Keywords: GPGPU, Fermi GPU, Kepler GPU, Performance Upperbound, Assembly GPGPU, Fermi GPU, Kepler GPU, Performance Upperbound, Assembly

In this paper, we studied the NVIDIA GPU architecture characteristics concerning the SGEMM routine and the potential peak performance of SGEMM on Fermi GPU. Guiding by the analysis, our SGEMM routine achieved about 11% (NN), 4.5% (TN), 3% (NT) and 9% (TT) better performance than cublas in CUDA 4.1 package for large matrices on GTX580 Fermi Card. We also described how to use native assembly language directly in the CUDA runtime source code.

Introduction

It seems that Moore's law is still alive for GPU. Nowadays, vendors can still put more computing cores per die and the GPU architecture is evolving rapidly. Fermi GPU is the present generation of NVIDIA GPU. We can find many articles about optimizing specific kernels on Fermi. However, since the architecture is changing very quickly, within a few months (the Kepler GPU is about to come), we may need to repeat the same work again. It would be interesting to generalize our present experience, find the critical architecture parameters, and bound the potential peak performance of those kernels. Hence, with these knowledge, it would be easier for us to move to the newer architecture.

In the report, we chose the GEMM kernel as our example (since we can only access the GTX580 Fermi card, so just SGEMM is evaluated). There are already a few articles about optimizing GEMM kernels on Fermi GPU [START_REF] Nath | An improved magma gemm for fermi gpus[END_REF] [6], and also auto-tuning framework is also presented [START_REF] Kurzak | Autotuning gemm kernels for the fermi gpu[END_REF]. We want to generalize the bounding factors for this computation-bound kernel and hope this knowledge can profit future optimization work on NVIDIA GPUs. Using the native assembly language to implement the device kernels, the four peak SGEMM performances (4800x4800 matrices) that we have achieved are 1165 GFLOPs (NN), 1157 GFLOPs (TN), 1173 GFLOPs (NT) and 1163 GLOPs (TT), which are 11% (NN), 4.5% (TN), 3% (NT) and 9% (TT) better cublas from CUDA 4.1 (released in Jan. 2012) respectively. Although the result is better than cublas, we have to see that the highest performance of 1173 GFLOPs (NT) is still only about 73.2% of the theoretical peak performance of our GTX580 card (1566KHz shader clock and 1604GLops theoretical peak performance).

GPU Architecture Characteristics Concerning SGEMM

General Architecture Difference between Fermi and GT200 GPU

The most obvious improvement from GT200 to Fermi GPU is that the SP number doubles, from 240 (GTX280, 65nm) to 512 (GTX580, 40nm). Each SM in Fermi GPU consists of 32 SPs instead of 8 SPs. For GTX280, each SM has 16KB shared memory and 16K 32bit registers. In GTX580, shared memory per SM increases to 48KB and the 32bit register number is 32K. However, if we consider the resource per SP, the on-die storage per SP actually decreases. Also for Fermi instruction set, there is a hard limit of maximum 63 registers per thread (for GT200 generation the limit is 127 registers per thread) since in the instruction coding, only 6 bits are left for one register. To reduce the effects of register spilling, Fermi GPU introduces L1 cache. Local writes are written back to L1. Global stores bypass L1 cache since multiple L1 caches are not coherent for global data.

L2 cache is also introduced in Fermi and could reduce the penalty of some irregular global memory accesses. A comparison between GT200 GPU and Fermi GPU is listed in Table 1.

Necessity of Register Blocking

In this work, we define a warp instruction as the same instruction shared by 32 threads in the same warp, and a thread instruction as the instruction executed by one thread. So a warp instruction launches 32 operations or consists of 32 thread instructions.

In GT200 generation, the scheduler can issue one warp instruction per core cycle and since there are 8 SPs per SM, SPs need 4 shader cycles, or 2 core cycles to process one warp instruction. Apparently, the issuing throughput is higher than the SP's processing throughput. That is why NVIDIA introduces the concept of dual-issue which means that the scheduler can use the However, in Fermi GPUs, SM are redesigned with 2 warp scheduler equipped with a dispatch unit, and 32 SPs. Each warp scheduler issue instructions to 16 SPs. With an issue rate of one warp instruction per core cycle, the scheduler's or dispatcher unit's ability could be fully utilized by SPs. The theoretical peak performance for math instructions comes from the SPs' performance. Also, the percentage of other instructions becomes an issue since if there are too many other kinds of auxiliary instructions like LD/ST, the scheduler doesn't have enough cycles to issue useful instructions like FFMA.

NVIDIA just launched the new Kepler GPU with 192 SPs on new designed SM or SMX. Now each SMX has 4 warp scheduler, each of which has 2 dispatch unit. And there is no shader cycle, all the units on die run at the core cycle. The hardware details are barely disclosed by now. Let us suppose that the dispatch unit's ability is the same as before, that is one warp instruction per core cycle. Then each core cycle, the dispatch units can supply 256 thread instructions. Clearly, if the 192 SPs all run at twice the core clock, there won't be enough instructions to feed them.

A comparison of the three generations of NVIDIA GPUs is listed in Table 1 and the data with a ? is not obtained or speculated.

For Fermi GPUs, according to the disassembly output of cuobjdump, data from shared memory cannot be used as operands of arithmetic instructions like FFMA. The instruction LDS is needed to first load data from shared memory to register and then math instructions can operate register operands. In the worst case, without register reuse, 2 LDS instructions are needed to fetch data for 1 FFMA instruction. So around 1/3 of the instructions devote to FLOPs. Also because the issuing power is not high enough, we can achieve less than 1/3 of the theoretical Blocking is a well known technique to better utilize memory hierarchy for scientific programs [3] [4]. With different blocking factors of registers, we illustrate the percentage of FFMA instructions in Figure 1. For example, if 6*6 register blocking is used on Fermi GPU, 75% of the instructions in the main loop of SGEMM are FFMA instructions. In this case, 6*6 registers are needed for C sub matrix per thread and 7 or more registers are needed for one line of A and B matrices per thread. Detailed analysis is in Section 4.1. In our analysis, the L1 and L2 cache don't devote to the peak performance. L1 cache is not coherent for different SMs and just reduces the latency of accessing some local data. For L2 cache, since the executing sequence of different C sub matrices cannot be controlled by software and if we consider that after some cycles, the blocks executing on different SMs are computing C sub matrices from random positions, there will be little chance for different SMs getting a hit in L2 cache.

Bound performance upperbound

DO { A d d
To make sure that there are enough threads per SM, equation 1 needs to be satisfied. Equation 4is the loose condition for register locking factor. B 2 R is the register needed to hold C sub matrix, B R is one line of A or B sub matrix, and at least another register is needed to load from one line of B or A sub matrix. Equation 7 is the tense condition for register blocking factor. If we want to prefetch data, at least 2 *

√ T B * B R * L T B
registers are needed. And we also need a few registers to store the address of matrices in global memory and shared memory.

T B * B 2 R * L is the size of a prefetched sub matrix for A or B. Equation 5 guarantees that it is easy to program the prefetching.

T SM ≥ 32 * #P ipeline_latency #Issue_latency (1)
R T ≤ R M ax (2)
Inria T SM * R T ≤ R SM (3)
B 2 R + B R + 1 < R T (4)
( T B * B R * L)%T B = 0 (5) Blk * 2 * T B * B R * L ≤ Sh SM (6) B 2 R + 2 * √ T B * B R * L T B + B R + 1 + R index ≤ R T (7)
B Sh = T B * B 2 R (8) #T hroughput_F actor = f (B R , #Issue_T P, #SP _T P, #LDS_T P ) (9)
P M emBound = #GlobalM em_bandwidth * B Sh (10)

P SM Bound = B 2 R B 2 R + B R * 2 * #T hroughput_F actor * P theoretical (11)
P potential = min(P M emBound , P SM Bound ) (12)
The #Throughput_Factor is a function of register blocking factor, throughput of SPs, LD/ST units and dispatch units. The function f for Fermi GPU is illustrated in Figure 1 and Figure 2.

For Fermi GPU, because of the hard limit of 63 registers per thread, according to the loose condition 4, the maximum register blocking factor is 7. Considering prefetching and using the tense condition 7, the maximum factor is only 6. The detailed register allocation is illustrated in Section 4.1. With this register blocking factor, the register resource per SM can support up to 512 threads. Using equation 5 and 8, we choose 256 threads per block to maximize shared memory blocking factor.

With a 6 register blocking factor, mixing LDS or LDS.64 with FFMA instructions, can both achieve close to 32 thread instructions per shader clock per shader. While using a LDS.64 instruction can increase the FFMA instruction percentage to 85.7% from 75%. Though LDS.128 instruction can provide higher percentage of 'useful' instructions, the instruction processing throughput is too low.

According to equation 10, 11 and 12, the performance is bounded by SM's processing throughput, and the potential peak is about 85% of the theoretical peak performance for SGEMM. And the main limitation factor comes from the nature of Fermi instruction set and the limited issuing throughput of schedulers.

Kepler GPU (GTX680) has already been released. The underlying architecture is different from Fermi GPU. However, our analysis method is still applicable to the new architecture. Since we do not own a new Kepler GPU card, we can only speculate some characteristics of Kepler. From unconfirmed information, the Kepler's instruction set is similar to Fermi GPU instruction set. That means that the limit of 63 registers per thread still exists. Thus, 6 register blocking is also applicable. And the register resource can support 1024 active threads per SM. The shared memory size can guarantee a blocking size around 96. Thus the main memory bandwidth is not bound. From section 2.2, we speculate that the issuing power of the dispatch units can support both SPs and LD/ST units running at full speed. However, since the FFMA/LDS.X ratio is 3:1 using LDS instructions, SPs may not work in full power (32 * 3 = 96). And using LDS.64 or LDS.128 instruction is not helping since they cannot increase the throughput to shared memory. To sum up all these factors, we speculate that the limitation factor for Kepler GPU to implement SGEMM would be the limited number of LD/ST units and SGEMM's potential peak performance is the only half of theoretical peak performance. So less than 50% of the theoretical peak performance on Kepler GPU is achievalbe for SGEMM and the key is to reduce the accesses to shared memory. Apparently, Kepler GPU is born for gaming instead of general computing. Other optimization methods are very similar to Fermi GPU.

Using Native Assembly Code in CUDA Rumtime Api Source Code

According to public materials, cuda binary file cannot be directly used in a project built with CUDA runtime api. The .cubin file can only be used with CUDA driver api. We find that loading .cubin file at runtime may degrade the performance. (We still need to look into that.) Besides, a lot of projects are programmed with CUDA runtime api. That restricts the usage of the code written in assembly language. We managed to integrate our .cubin code into a CUDA runtime project. Since NVIDIA never promises to support this, the method is not guaranteed to be working on future CUDA toolkit. The CUDA binary is generated using an assembly tool for Fermi GPU called Asfermi [1];

1. Add -v -keep to nvcc options, so that all the intermediate files are saved and we can have all the compiling steps.

2. Write one .cu files, eg. kernel.cu, which has one dummy kernel function with the same device function name as your CUDA binary code.

3. Add kernel.cu into the project.

4. Make the project and collect all the compiling command line.

5. Replace the generated kernel.cubin with the one that you implement. Regenerate the kernel.fatbin.c file, like fatbin -key="0de11f308b8b4b9d" "-image=profile=sm_20,file=kernel.cubin" -embedded-fatbin="kernel.fatbin.c" 6. Regenerate the kernel.cu.cpp and then kernel.cu.o according to the original command line information.

7. Now we can build the whole project.

LDS.64 & LDS.128

According to our analysis, it is essential for achieving better performance to minimize auxiliary instructions, especially LDS instruction. The assembly code for CUDA sm_20 provides LDS.64 and LDS.128 instructions to load 64bit and 128bit data at one time per thread. To use these two instructions, the start address in shared memory should be 64bit and 128bit aligned. Also, the index of registers needs to be 2 and 4 aligned respectively. In the assembly language level, it is possible for us to carefully design the data layout and register allocation to satisfy these requirements. 2, with the ratio of 0, it is the processing throughput for LDS.X instructions. Clearly, LDS.128 has only the throughput of 2 thread instructions per shader cycle. That means, the LD/ST units need 16 shader cycles to process one LDS.128 instruction.

If 6 register blocking is used, the FFMA/LDS.X ratios are 3:1, 6:1, and 12:1 for LDS, LDS.64 and LDS.128 respectively. According to our experiments, the processing throughput for one SM in these cases are 31.3, 30.4 and 24.5 thread instructions per shader clock. Apparently, if the accesses to shared memory are implemented all by LDS.128, in the best case we can only achieve around 76.5% of SMs' throughput.

Apparently, it is not worth the programming efforts to mix FFMA with LDS.128. Especially, a lot of padding has to be used since we are using 6 register blocking.

Assembly Level Optimization

Eliminate Register Spilling

The register resource is 32K 32-bit registers per SM for Fermi GPU and each thread uses maximum 63 registers. The register R1 normally is occupied as stack pointer. According to our analysis, the per-thread register number with prefetching is

B 2 R + 2 * √ T B * B R * L T B + B R + 1 + R index .
With the register blocking factor of 6 for Fermi GPU, the register allocation of our implementation is as the following. Also, we used 32bit address to save address registers. So, in all 36 + 12 + 6 + 2 + 2 + 1 + 2 + 2 = 63 registers are used. Since we don't need thread stack, R1 is used to store the loop end condition in our code. Thus, register spilling can be fully eliminated.

Instruction Reordering

1. In the main loop, between the 2 barriers are all shared memory accesses. By moving address calculation from start of the loop to mix with the shared memory accesses, we can achieve better performance.

2. By interleaving prefetching from global memory with FFMA and LDS instructions can benefit performance too.

Optimization Result

As in Figure 3, it is the performance comparison between our implementation (alf_sgemm), cublas from CUDA 4.1 and MAGMA library [START_REF] Nath | An improved magma gemm for fermi gpus[END_REF].

Conclusion

In this work, we analyzed the potential peak performance bound on Fermi GPU for SGEMM. We showed that the nature of Fermi instruction set and the limited issuing throughput of schedulers are the main limitation factors for approaching the theoretical peak performance. Guiding by the analysis, we further optimized the four SGEMM kernels. The 4 kernels achieved 11% (NN), 4.5% (TN), 3% (NT) and 9% (TT) better performance than cublas from CUDA 4.1. We also described how to use assembly code directly in a CUDA runtime source project and discussed the possible optimization methods for SGEMM on Kepler GPU and the possible peak performance that we can achieve.

Figure 1 :

 1 Figure 1: FFMA instruction percentage in SGEMM main-loop with different Blocking Factors

Figure 3 :

 3 Figure 3: SGEMM Performance

Table 2 :

 2 Architecture and Algorithm Parameters

	Prefetch_A_GlobalMem ( ) ;
	Prefetch_B_GlobalMem ( ) ;
	LDS_A_SharedMem ( ) ;
	LDS_B_SharedMem ( ) ;
	FFMA_C_A_B( ) ;
	__syncthreads ( ) ;
	STS_Prefetched_A ( ) ;
	STS_Prefetched_B ( ) ;
	__syncthreads ( ) ;
	}
	While (ADDR_B < LOOP_END)
	Listing 1: GEMM Main Loop Pseudocode
	The main loop of the SGEMM routine with prefetching is illustrated in Listing 1. At each
	iteration, first addresses of A, B are recalculated according to the stride; second, start to load the

r e s s _ C a l c u l a t i o n ( ) ; A, B data for next iteration into registers; third, load data from shared memory and accumulate to C; after a barrier, store the prefetched data into shared memory.

  Figure 2: Instruction Processing Rate Mixing FFMA and LDS However, LDS.128 normally leads to 2-way shared-memroy bank conflict. As in Figure

		32		LDS (3:1, 31.3)
		28			LDS.64 (6:1, 30.4)
	Processing Rate	12 16 20 24						LDS.128 (12:1, 24.5)						LDS LDS.64 LDS.128
		8
		4
		0
		0	2	4	6	8	10	12	14	16	18	20	22	24	26	28	30	32
									FFMA/LDS.X ratio
																		Inria

RR n°7923

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Abderrahmane Bendali

† M'

Barek Fares

Sophie Laurens

Sébastien Tordeux

Numerical study of acoustic multiperforated plates

It is rather classical to model multiperforated plates by approximate impedance boundary conditions. In this article we would like to compare an instance of such boundary conditions obtained through a matched asymptotic expansions technique to direct numerical computations based on a boundary element formulation in the case of linear acoustic.

Introduction

Multiperforated boundaries are involved in several industrial devices (liners [START_REF] Lidoine | Numerical prediction of SDOF-Perforated Plate Acoustic Treatment Impedance. Part 1: Linear domain[END_REF], mufflers [START_REF] Bonnet-Bendhia | Mathematical analysis of the acoustic diffraction by a muffler containing perforated ducts[END_REF], combustion chambers [START_REF] Mendez | Acoustic modeling of perforated plates with bias flow for Large-Eddy Simulations[END_REF], noise barriers ...). The occurrence of a high density of holes of small size has a significant impact on the noise generated by the devices.

Until recently, it was almost impossible to predict the effect of a multiperforated plate by a direct numerical computation. To overcome this bottleneck, many authors have proposed to replace the perforated wall by an effective impedance boundary condition of transmission type, see [START_REF] Tuck | Matching problems involving flow through small holes[END_REF] for the case of linear acoustic. This approach is suitable for numerical computation since it requires no mesh refinement in the neighborhood of the holes. The rapid increase of numerical computational abilities allows for the validation of this 1 method by means of high performance computing. In this article, we would like to address such an issue using the parallel boundary element library CESC of CERFACS.

In section 1, we describe the considered geometry and state the boundaryvalue problem obtained by assuming that the equations of linearized acoustics are valid everywhere. In section 2, we briefly review some effective boundary conditions previously used for modeling multiperforated plates. Numerical experiments are presented in section 3.

Statement of the Problem 1.Description of the propagation domain

We consider an infinite waveguide T with rectangular section

T = x ∈ R 3 : x 1 ∈]0, h 1 [ and x 2 ∈]0, h 2 [ .
This domain is split into two parts by a multiperforated plate with zero thickness

P = x ∈ R 3 : x 1 ∈]0, h 1 [, x 2 ∈]0, h 2 [ and x 3 = 0 .
The wall P is perforated by a grating of characteristic lengths η 1 and η 2 of holes Θ m of radius δ > 0

           x m = η 1 2 , η 2 2 , 0 + m 1 + m 2 α η 1 , m 2 η 2 , 0 with α ∈ [0, 1[ Θ m = {x ∈ R 3 : x -x m < δ and x 3 = 0} with m ∈ S = (m 1 , m 2 ) ∈ Z 2 : Θ m ⊂ P . (1)
η 2 η 1 αη 1 h 2 h 1 x 1 x 2 Figure 1: The multiperforated plate
The spacing between two holes is large enough so that short range interactions can be neglected, i.e., δ/min(η 1 , η 2 ) = o(1). We can thus define the propagation domain:

Ω = T \ P perf with P perf = P \ m∈S Θ m .

Linear acoustics modeling

The time variable and the speed of sound are respectively denoted by t and c. Let us recall that in the context of linear acoustic the pressure field p : Ω × R + → R is governed by a scalar wave equation

∂ 2 p ∂t 2 (x, t) -c 2 ∆p(x, t) = 0 in Ω.
The disturbed velocity vector field v : Ω × R + → R 3 is linked to the pressure by linearized Euler's equations

ρ 0 ∂v ∂t (x, t) + ∇p(x, t) = 0 in Ω,
where ρ 0 is the density of the fluid at rest. The boundary of the physical domain is supposed to be a perfectly rigid wall

v • n = 0 on ∂Ω = ∂T \ P \ m∈S Θ m ,
with n the normal unit vector outwardly directed to Ω. Considering a wave with pulsation ω and wave number k = ω/c, we denote also by p and v the phasors of respectively the pressure and the velocity p(x, t) = ℜ p(x) exp(-iωt) and v(x, t) = ℜ v(x) exp(-iωt) .

In the sequel we consider the well-posed boundary-value problem

                     Find p : Ω -→ C satisfying χp ∈ H 1 (Ω) for all χ ∈ D(R 3 ), ∆p(x) + k 2 p(x) = 0 in Ω, ∂p ∂n (x) = 0 on ∂Ω, p(x) -exp(ikx 3
) is outgoing, see [START_REF] Goldstein | A finite element method for solving Helmholtz type equations in waveguides and other unbounded domains[END_REF] or [START_REF] Lenoir | The localized finite element method and its application to the two-dimensional sea-keeping problem[END_REF] for example,

which provides a model for the time-harmonic waves standing in this acoustic system when it is fed by the propagating mode coming from x 3 = -∞.

Remark 1. For a waveguide with small height and width (kh 1 < π and kh 2 < π) the harmonic pressure p can be expressed as

   p(x, t) = exp(ikx 3 ) + R exp(-ikx 3 ) + w + (x) for x 3 < 0, p(x, t) = T exp(ikx 3 ) + w -(x) for x 3 > 0.
Parameters R and T are respectively the reflection and transmission coefficients. Both functions w -and w + are evanescent in the x 3 direction (exponential decay):

w -(x) = o x 3 →-∞
(1) and w + (x) = o

x 3 →+∞
(1).

The approximate transmission model for multiperforated plates

For a large number of holes, even when resorting to high performance computing, the direct numerical solution of problem ( 2) is generally out of reach.

To overcome this difficulty a solution consists in modeling the multiperforated plate by an equivalent transmission condition (also called impedance condition) which requires no mesh refinement. These approximate conditions have been well-studied by the acoustic community both from the theoretical and experimental point of view. It can either be achieved using heuristic approaches based on a physical reasonning [START_REF] Tuck | Matching problems involving flow through small holes[END_REF] or by means of asymptotic analysis, see [START_REF] Bonnet-Bendhia | Simulation of muffler's transmission losses by a homogenized finite element method[END_REF], [START_REF] Bonnet-Bendhia | Mathematical analysis of the acoustic diffraction by a muffler containing perforated ducts[END_REF] or [START_REF] Chechkin | On boundary-value problems for the Laplacian in bounded and in unbounded domains with perforated boundaries[END_REF] and references therein. The last approach benefits from a rigorous framework based on the matching of asymptotic expansions techniques [START_REF] Il'in | Matching of asymptotic expansions of solutions of boundary value problems, volume 102 of Translations of Mathematical Monographs[END_REF].

The conclusion of these studies can be physically interpreted as follows: the multiperforated plate can be modeled as an impenetrable moving membrane of surface mass ρ S . This leads to a system of partial differential equations posed on the domain Ω = T \ P governing the pressure field of the fluid p, the flow velocity vector field of the fluid v and the surface velocity of the membrane V = (0, 0,

V 3 )                                           
Acoustic equations in the fluid:

∂ 2 p ∂t 2 (x, t) -c 2 ∆ p(x, t) = 0 in Ω, ρ 0 ∂ v ∂t (x, t) + ∇ p(x, t) = 0 in Ω, v(x, t) • n = 0 on ∂T.
Impenetrability of the membrane:

v 3 (x 1 , x 2 , 0 -, t) = v 3 (x 1 , x 2 , 0 + , t) = v 3 (x 1 , x 2 , t).
Second Newton's law for the membrane:

ρ S ∂ v 1 ∂t (x 1 , x 2 , t) = p(x 1 , x 2 , 0 -, t) -p(x 1 , x 2 , 0 + , t).
In the case of an incident harmonic wave coming from x 3 = -∞, this system simplifies into a problem involving only p :

Ω -→ C                        χ p ∈ H 1 ( Ω) for all χ ∈ D(R 3 ), ∆ p(x) + k 2 p(x) = 0 in Ω, ∂ p ∂n (x) = 0 on ∂T, ∂ p ∂x 3 (x 1 , x 2 , 0 -or 0 + ) = K R A p(x 1 , x 2 , 0 + ) -p(x 1 , x 2 , 0 -) , p(x) -exp(ikx 3 ) is outgoing, (3)
expressed in terms of the so-called Rayleigh conductivity K R of the hole and the surface of an unit-cell of the network A = η 1 η 2 . These quantities are related to the surface mass of the membrane by:

ρ S = ρ 0 A K R .
Remark 2 (Physical interpretation of the transmission condition). In the theory of acoustics, it is well known that the effect of a small perforation in a rigid plate can be modeled by two acoustic point sources located on each side of the plate, see [START_REF] Tuck | Matching problems involving flow through small holes[END_REF] for example. The algebraic fluxes of fluid coming from these sources Φ -and Φ + are opposite (due to the conservation of mass) and related to the difference of pressure on the two sides of the plate p + -p -by an analogous relation to the second Newton's law:

ρ 0 dΦ - dt = -ρ 0 dΦ + dt = K R p + -p -.
The multiplicative constant K R is traditionally called the Rayleigh conductivity, see [START_REF] Rayleigh | The theory of sound[END_REF]. It is given by

K R = 2δ
in the case of linear acoustics for a circular hole of radius δ in a plate of negligible thickness.

A fluid crossing a grating of holes has a complex structure in the neighborhood of the perforations but can be seen as a homogeneous fluid at large scale. The velocity of fluid is then determined by distributing over one unit cell of the network the flux coming from one hole :

v 3 (x 1 , x 2 , 0 -, t) = - Φ - A and v 3 (x 1 , x 2 , 0 + , t) = Φ + A .
This leads to the transmission conditions

             v 3 (x 1 , x 2 , 0 -, t) = v 3 (x 1 , x 2 , 0 + , t), A ∂ v 3 ∂t (x 1 , x 2 , 0 -, t) = ρ 0 K R p(x 1 , x 2 , 0 -, t) -p(x 1 , x 2 , 0 + , t) , A ∂ v 3 ∂t (x 1 , x 2 , 0 + , t) = ρ 0 K R p(x 1 , x 2 , 0 -, t) -p(x 1 , x 2 , 0 + , t) . (4)
Remark 3. The solution of the approximate model can be analytically computed. Its solution depends only on the x 3 coordinate and is given by

p -(x) = exp(ikx 3 ) + R exp(-ikx 3 ) for x 3 < 0, p + (x) = T exp(ikx 3 ) for x 3 > 0. ( 5
)
Inserting (5) in the transmission conditions (4) leads to the reflection and transmission coefficients:

R = - ikA K R 2 - ikA K R and T = 2 2 - ikA K R .

A numerical experiment

We consider two multiperforated plates with almost the same geometrical characteristic lengths h 1 = 5 cm, h 2 = 5 cm, η 1 = 5 mm, η 2 = 5 mm, δ = 0.225 mm.

In the first configuration, the grating is straight, α = 0 (see [START_REF] Bonnet-Bendhia | Simulation of muffler's transmission losses by a homogenized finite element method[END_REF] and Figure 1 for the definition of α) and tilted in the second one, α = 1/2, see Figure 2. The meshes, see Figures 3, were generated with the proprietary licensed software Altair HyperMesh library (http://www.hyperworks.fr) and its scripting language to ensure a local refinement in the neighborhood of each hole, see Figure 4. The outgoing condition is approximated by truncating the domain at h 3 = 2cm from the multiperforated plate, see Figure 5. On the two sections x 3 = -h 3 and x 3 = h 3 is prescribed the following accurate absorbing boundary conditions :

       ∂p ∂x 3 (x 1 , x 2 , -h 3 ) + ikp(x 1 , x 2 , -h 3 ) = 2ik exp(ikh 3 ), ∂p ∂x 3 (x 1 , x 2 , h 3 ) -ikp(x 1 , x 2 , h 3 ) = 0.
These conditions are rather accurate since they are obtained by neglecting the evanescent modes which are exponentially decaying. The numerical experiments have been achieved with the boundary element library CESC. We have chosen to use a double layer formulation [START_REF] Sauter | Boundary element methods[END_REF] to discretize the system of equations ( 2) and numerically compute its reflection coefficient R. In order to check the convergence of the method, we have compared the numerical results with approximately 4 × 10 4 degrees of freedom (d.o.f.) to the results with 1.6 × 10 5 d.o.f. Both numerical approximations of the reflection coefficients coincide up to the third digit.

x 3 = 0 h 2 h 1 h 3 h 3
We have considered 26 frequencies (from 500 Hz to 3000 Hz with a 100 Hz step). The determination of the transmission and reflection coefficients at each frequency requires to form and invert a linear system with a dense matrix of approximately 40,000 unknowns. It was carried on a parallel platform having 400 processors. It required a rather large computational time (5.2 × 10 6 s i.e. 1.3 × 10 4 s • cpu -1 ).

The modulus |R| and the phase angle of the reflection coefficients respectively obtained from the direct numerical solution and the homogenized approximate model are depicted in Figure 6 for the straight case (α = 0) and in Figure 7 for the tilted case (α = 1/2).

The conclusions that can be drawn from these plots are the following:

• As predicted by the theory, there is almost no difference between the straight case and the tilted case.

• Even if the overall behavior of the reflection coefficient is well predicted, an error of 10 % on its modulus and of 5 % on its phase angle can be observed.

Conclusion

A numerical validation of the approximate transmission condition modeling the reflection and the transmission of acoustic waves by a multiperforated plate has been carried out. The modulus and the phase angle for the reflection coefficients obtained through a direct numerical solution and the use of the homogenized transmission condition have exactly the same overall behavior even if a discrepancy between the two results gives rise to an error of about 10 %. To reduce this discrepancy, it would be interesting in our opinion to carry out further investigations.

First, the homogenized transmission conditions can be improved by taking into account the short rang interactions. These objectives can be achieved by performing a multi-scale in place of a two-scale asymptotic analysis.

Moreover, in most applications (combustion chamber, liners), linear acoustics is not generally used for modeling the behavior of the acoustic field in the proximity of the holes. Reproducing this approach for these more complicated models seem to be rather challenging. One can refer to [START_REF] Sanchez-Hubert | Acoustic fluid flow through holes and permeability of perforated walls[END_REF] for viscosity and to [START_REF] Howe | On the theory of unsteady high reynolds number flow through a circular aperture[END_REF][START_REF] Howe | Acoustics of fluid-structure interactions[END_REF][START_REF] Hughes | The absorption of sound by perforated linings[END_REF] for aeroacoustic effect.

The last issue is the choice of the numerical method. We have chosen to use a direct solution by a plain boundary element method without accelerating it by a Fast Multipole Method (FMM) (see, e.g., [START_REF] Gumerov | Fast Multipole Methods for the Helmholtz Eqaution in Three Dimensions[END_REF]). This choice is motivated by the fact that the problem being considered here is in the low frequencies range and at the best of authors knowledge no completely satisfactory FMM is available in this case.

Figure 2 :

 2 Figure 2: Multiperforated plates with α = 0 and α = 1 2

Figure 3 : 2 Figure 4 :

 324 Figure 3: The meshes of the multiperforated plates for α = 0 and α = 1/2

Figure 5 :

 5 Figure 5: The truncated domain of computation

Figure 6 :Figure 7 :

 67 Figure 6: Comparison of the reflection coefficient R = |R| exp(iϕ) of the direct computations (green) and of the approximate model (blue) for straight multiperforated plate

Acknowledgments.

This work was supported by the French National Research Agency under grant no. ANR-08-SYSC-001.