Charles Bouveyron

Camille Brunet

Model-Based Clustering of High-Dimensional Data: A review

Keywords: Model-based clustering, high-dimensional data, dimension reduction, regularization, parsimonious models, subspace clustering, variable selection, softwares, R packages

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Clustering is a data analysis tool which aims to group data into several homogeneous groups. The clustering problem has been studied for years and usually occurs in applications for which a partition of the data is necessary. In particular, more and more scientific fields require to cluster data in the aim to understand or interpret the studied phenomenon. Earliest approaches were based on heuristic or geometric procedures. They relied on dissimilarity measures between pairs of observations. A popular dissimilarity measure is based on the distance between groups, previously introduced by Ward [START_REF] Ward | Hierarchical groupings to optimize an objective function[END_REF] for hierarchical clustering. In the same way, the k-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] is perhaps the most popular clustering algorithm among the geometric procedures. Clustering was also defined in a probabilistic framework, allowing to formalize the notion of clusters through their probability distribution. One of the main advantages of this probabilistic approach is in the fact that the obtained partition can be interpreted from a statistical point of view. The first works on finite mixture models were from Wolfe [START_REF] Wolfe | Object cluster analysis of social areas[END_REF], Scott et al. [START_REF] Scott | Clustering methods based on likelihood ratio criteria[END_REF] and Duda et al. [START_REF] Duda | Pattern classification[END_REF]. Since then, these models have been extensively studied and, thanks to works such as those of McLachlan et al. [START_REF] Mclachlan | Mixture models : inference and applications to clustering[END_REF][START_REF] Mclachlan | Finite Mixture Models[END_REF], Banfield & Raftery [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF] or Fraley & Raftery [START_REF] Fraley | Algorithms for model-based Gaussian Hierarchical Clustering[END_REF], [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF], model-based clustering has become a popular and reference technique.

Nowadays, the measured observations in many scientific domains are frequently highdimensional and clustering such data is a challenging problem for model-based methods. Indeed, model-based 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 over-parametrized in high-dimensional spaces.

Furthermore, in several applications, such as mass spectrometry or genomics, the number of available observations is small compared to the number of variables and such a situation increases the problem difficulty. Since the dimension of observed data is usually higher than their intrinsic dimension, it is theoretically possible to reduce the dimension of the original space without loosing any information. For this reason, dimension reduction methods are frequently used in practice to reduce the dimension of the data before the clustering step. Feature extraction methods, such as principal component analysis (PCA), or feature selection methods are very popular. However, dimension reduction usually does not consider the classification task and provide a sub-optimal data representation for the clustering step. Indeed, dimension reduction methods imply an information loss which could have been discriminative.

To avoid the drawbacks of dimension reduction, several approaches have been proposed to allow model-based methods to efficiently cluster high-dimensional data. This work proposes to review the alternatives to dimension reduction for dealing with high-dimensional data in the context of model-based clustering. Earliest approaches include constrained and parsimonious models or regularization. More recently, subspace clustering techniques and variable selection techniques have been proposed to overcome the limitations of previous approaches. Subspace clustering techniques are based on probabilistic versions of the factor analysis model. This modeling allows to cluster the data in low-dimensional subspaces without reducing the dimension. Conversely, variable selection techniques do reduce the dimension of the data but select the variables to retain regarding the clustering task. Both techniques turn out to be very efficient and their practical use will be discussed as well in this article. This article is organized as follows. Section 2 briefly recalls the bases of mixture modeling and its inference with the EM algorithm. Section 3 introduces the curse of dimensionality in model-based clustering. Approaches based on dimension reduction, regularization and parsimonious models are reviewed in Section 4. Then, Section 5 and 6 present respectively the approaches based on subspace clustering and variable selection. Existing softwares for model-based clustering of high-dimensional data are also reviewed in Section 7 and their practical use is discussed in Section 8. Finally, some concluding remarks are made in Section 9.

The mixture model and the EM algorithm

This section first recalls the bases of mixture modeling and its inference with the expectation-maximization (EM) algorithm.

The mixture model

Let us consider a data set of n observations {y 1 , . . . , y n } ∈ R p that one wants to divide into K homogeneous groups. The aim of clustering is to determine, for each observation y i , the value of its unobserved label z i such that z i = k if the observation y i belongs to the kth cluster. To do so, model-based clustering [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF][START_REF] Mclachlan | Finite Mixture Models[END_REF] considers the overall population as a mixture of the groups and each component of this mixture is modeled through its conditional probability distribution. In this context, the observations {y 1 , . . . , y n } ∈ R p are assumed to be independent realizations of a random vector Y ∈ R p whereas the unobserved labels {z 1 , ..., z n } are assumed to be independent realizations of a random variable Z ∈ {1, ..., K}. The set of pairs {(y i , z i )} n i=1 is usually referred to as the complete data set. By denoting by g the probabilistic density function of Y , the finite mixture model is:

g(y) = K k=1 π k f k (y), (1)
where π k (with the constraint K k=1 π k = 1) and f k respectively represent the mixture proportion and the conditional density function of the kth mixture component. Furthermore, the clusters are often modeled by the same parametric density function in which case the finite mixture model is:

g(y) = K k=1 π k f (y; θ k ), (2)
where θ k is the parameter vector for the kth mixture component. For a set of observations y = {y 1 , ..., y n }, the log-likelihood of this mixture model is then:

ℓ (θ; y) = n i=1 log K k=1 π k f (y i ; θ k ) . (3)
However, the inference of this model cannot be directly done through the maximization of the likelihood since the group labels {z 1 , ..., z n } of the observations are unknown. Indeed, due to the exponential number of solutions to explore, the maximization of equation ( 3) is unfortunately intractable, even for limited numbers of observations and groups. Before introducing the most popular inference algorithm used in this context, let us introduce the complete log-likelihood:

ℓ c (θ; y, z) = n i=1 K k=1 z ik log (π k f (y i ; θ k )) ,
where z ik = 1 if the ith observation belongs to the kth cluster and z ik = 0 otherwise.

The EM algorithm

Although it is not specifically dedicated to mixture models, the expectation-maximization (EM) algorithm, proposed by Dempster et al. [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], is certainly the most popular technique for inferring mixture models. The EM algorithm iteratively maximizes the conditional expectation of the complete log-likelihood:

E [ℓ c (θ; y, z) |θ * ] = K k=1 n i=1 t ik log (π k f (y i ; θ k )) ,
where t ik = E [z = k|y i , θ * ] and θ * is a given set of the mixture parameters. From an initial solution θ (0) , the EM algorithm alternates two steps: the E-step and the M-step. First, the expectation step (E-step) computes the expectation of the complete log-likelihood E ℓ c (θ; y, z) |θ (q) conditionally to the current value of the parameter set θ (q) . Then, the maximization step (M-step) maximizes E ℓ c (θ; y, z) |θ (q) over θ to provide an update for the parameter set. This algorithm therefore forms a sequence θ (q) q≥1 which satisfies, for each q ≥ 1:

θ (q+1) = arg max θ E ℓ c (θ; y, z) |θ (q) .
One of the most outstanding properties of the EM algorithm is that it guarantees an improvement of the likelihood function at each iteration. Each update of the parameter set, resulting from an E-step followed by an M-step, is guaranteed to increase the log-likelihood function. In particular, Wu [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF] proved that the sequence of θ (q) q≥1 converges to a local optimum of the likelihood. For further details on the EM algorithm, the reader may refer to [START_REF] Mclachlan | The EM algorithm and extensions[END_REF].

The two steps of the EM algorithm are iteratively applied until a stopping criterion is satisfied. The stopping criterion may be simply |ℓ(θ (q) ; y)ℓ(θ (q-1) ; y)| < ε where ε is a positive value to provide. It would be also possible to use the Aitken's acceleration criterion [START_REF] Lindsay | Mixture models: Theory, geometry and applications[END_REF] which estimates the asymptotic maximum of the likelihood and allows to detect in advance the algorithm convergence. Once the EM algorithm has converged, the partition {ẑ 1 , . . . , ẑK } of the data can be deduced from the posterior probabilities t ik = P (Z = k|y i , θ) by using the maximum a posteriori (MAP) rule which assigns the observation y i to the group with the highest posterior probability.

The Gaussian mixture model

Among the possible probability distributions for the mixture components, the Gaussian distribution is certainly the most used for both theoretical and computational reasons. Let us however notice that several recent works focused on different distributions such as the skew normal [START_REF] Lin | Finite mixture modelling using the skew normal distribution[END_REF], asymmetric Laplace [START_REF] Browne-R | Mixtures of shifted asymmetric Laplace distributions[END_REF] or t-distributions [START_REF] Andrews | Model-based clustering, classification, and discriminant analysis via mixtures of multivariate t-distributions[END_REF][START_REF] Lee | Finite mixtures of multivariate skew t-distributions: some recent and new results[END_REF][START_REF] Mclachlan | Robust cluster analysis via mixtures of multivariate t-distributions[END_REF]. Nevertheless, the density function f (y; θ k ) is most commonly assumed to be a multivariate Gaussian density φ(y; θ k ) parametrized by its mean µ k and its covariance matrix Σ k , such that the density function of Y can be written as:

g(y; θ) = K k=1 π k φ(y; θ k ), (4)
where:

φ(y; θ k ) = 1 (2π) p/2 |Σ k | 1/2 exp - 1 2 (y -µ k ) t Σ -1 k (y -µ k )
is the multivariate Gaussian density with parameter θ k = (µ k , Σ k ). This specific mixture model is usually referred in the literature as the Gaussian mixture model (GMM). In this case, the EM algorithm has the following form, at iteration q: E-step. This step aims to compute the expectation of the complete log-likelihood conditionally to the current value of the parameter θ (q-1) . In practice, it reduces to the com-

putation of t (q) ik = E[z i = k|y i , θ (q-1)
]. Let us also recall that t

(q)
ik is as well the posterior probability P (z i = k|y i , θ (q-1) ) that the observation y i belongs to the kth component of the mixture under the current model. Using Bayes' theorem, the posterior probabilities t (q) ik , i = 1, ..., n, k = 1, ..., K, can be expressed as follows:

t (q) ik = π (q-1) k φ(y i , θ (q-1) k ) K l=1 π (q-1) l φ(y i |θ (q-1) l ) , (5)
where π

(q-1) k and θ (q-1) k = µ (q-1) k , Σ (q-1) k
are the parameters of the kth mixture component estimated at the previous iteration.

M-step.

This step updates the model parameters by maximizing the conditional expectation of the complete log-likelihood. The maximization of E ℓ c (θ; y, z) |θ (q-1) conduces to an update of the mixture proportions π k , the means µ k and the covariance matrices Σ k as follows, for k = 1, ..., K:

π(q) k = n (q) k n , (6) μ
(q) k = 1 n (q) k n i=1 t (q) ik y i , (7)
Σ(q) k = 1 n (q) k n i=1 t (q) ik y i - μ(q) k y i - μ(q) k t , ( 8
)
where n

(q) k = n i=1 t (q)
ik .

The curse of dimensionality in model-based clustering

Before to present classical and recent methods for high-dimensional data clustering, we focus in this section on the causes of the curse of dimensionality in model-based clustering.

The origins: the Bellman' sentence

When reading research articles or books related to high-dimensional data, it is very likely to find the term "curse of dimensionality" to refer to all problems caused by the analysis of high-dimensional data. This term was introduced by R. Bellman in the preface of his book [START_REF] Bellman | Dynamic Programming[END_REF] promoting dynamic programming. To illustrate the difficulty to work in high-dimensional spaces, Bellman recalled that if one considers a regular grid with step 0.1 on the unit cube in a 10-dimensional space, the grid is made of 10 10 points. Consequently, the search for the optimum of a given function in this unit cube requires 10 This paragraph will indeed show that the Bellman's thought was corrected since, at least for clustering, high dimensions have nice properties which do allow to obtain significant results.

The curse of dimensionality in model-based clustering

In the context of model-based clustering, the curse of dimensionality takes a particular form. Indeed, model-based clustering methods require the estimation of a number of parameters which directly depends on the dimension of the observed space. If we consider the classical Gaussian mixture model for instance, the total number of parameters to estimate is equal to:

ν = (K -1) + Kp + Kp(p -1)/2,
where (K -1), Kp and Kp(p -1)/2 are respectively the numbers of free parameters for the proportions, the means and the covariance matrices. It turns out that the number of parameters to estimate is therefore a quadratic function of p in the case of the Gaussian mixture model and a large number of observations will be necessary to correctly estimate those model parameters. Furthermore, a more serious problem occurs in the EM algorithm when computing the posterior probabilities

t ik = E[Z = k|y i , θ] which depend, in the GMM context, on the quantity H k (x) = -2 log(π k φ(x; µ k , Σ k ))
. Indeed, the computation of H k , which can be rewritten as:

H k (x) = (x -µ k ) t Σ -1 k (x -µ k ) + log(det Σ k ) -2 log(π k ) + p log(2π), (9)
requires the inversion of the covariance matrices Σ k , k = 1, ..., K. Consequently, if the number of observations n is small compared to ν, the estimated covariance matrices Σk are ill-conditioned and their inversions conduce to unstable classification functions. In the worst case where n < p, the estimated covariance matrices Σk are singular and modelbased clustering methods cannot be used at all. Unfortunately, this kind of situation tends to occur more and more frequently in Biology (DNA sequences, genotype analysis) or in computer vision (face recognition) for instance.

The curse of dimensionality has been also exhibited by [START_REF] Pavlenko | On feature selection, curse of dimensionality and error probability in discriminant analysis[END_REF][START_REF] Pavlenko | Effect of dimensionality on discrimination[END_REF] in the Gaussian case with the estimation point of view. Let us consider the estimation of the normalized trace τ (Σ) = tr(Σ -1 )/p of the inverse covariance matrix Σ of a multivariate Gaussian distribution N (µ, Σ). The estimation of τ from a sample of n observations {x 1 , ..., x n } conduces to:

τ (Σ) = τ ( Σ) = 1 p tr( Σ-1 ),
and its expectation is:

E[ τ (Σ)] = 1 - p n -1 -1 τ (Σ). Consequently, if the ratio p/n → 0 when n → +∞, then E[ τ (Σ)] → τ (Σ). However, if the dimension p is comparable with n, then E[ τ (Σ)] → cτ (Σ)
when n → +∞, where c = lim n→+∞ p/n. We refer to [START_REF] Pavlenko | On feature selection, curse of dimensionality and error probability in discriminant analysis[END_REF][START_REF] Pavlenko | Effect of dimensionality on discrimination[END_REF] for further details on the effect of the dimensionality on classification in the asymptotic framework, i.e. p and n → +∞.

The blessing of dimensionality in clustering

Hopefully, as expected by Bellman, high-dimensional spaces have specific features which could facilitate their exploration. Several authors, such as [START_REF] Hall | Geometric representation of high dimension, low sample size data[END_REF][START_REF] Murtagh | The remarkable simplicity of very high dimensional data: application of model-based clustering[END_REF], have shown that, in the context of clustering, high-dimensional spaces do have useful characteristics which ease the classification of data in those spaces. In particular, Scott and Thompson [START_REF] Scott | Probability density estimation in higher dimensions[END_REF] showed that high-dimensional spaces are mostly empty. The experiment suggested by Huber [START_REF] Huber | Projection pursuit[END_REF] consists in drawing realizations of a p-dimensional random vector X with uniform probability distribution on the hypersphere of radius 1. The probability that a realization x i of this experiment belongs to the shell between the hypersphere of radius 0.9 and the unit hypersphere is therefore:

P (x i ∈ S 0.9 (p)) = 1 -0.9 p .
In particular, the probability that x i belongs to the shell between the hypersphere of radius 0.9 and the unit hypersphere in a 20-dimensional space is roughly equals to 0.88. Therefore, most of the realizations of the random vector X live near a p -1 dimensional subspace and the remaining of the space is mostly empty. This suggests that clustering methods should model the groups in low-dimensional subspaces instead to model them in the whole observation space. Furthermore, it seems reasonable to expect that different groups live in different subspaces and this may be a useful property for discriminating the groups. Subspace clustering methods, presented in Section 6, exploit this specific characteristic of high-dimensional spaces.

Earliest approaches

Earliest approaches to deal with the clustering of high-dimensional data can be split into three families: dimension reduction methods, regularization methods and parsimonious methods.

Dimension reduction

Approaches based on dimension reduction assume that the number p of measured variables is too large and, implicitly, that the data at hand live in a space of lower dimension, let us say d < p. Once the data projected in a low-dimensional space, it is then possible to apply the EM algorithm on the projected observations to obtain a partition of the original data.

The most popular linear method used for dimension reduction is certainly principal component analysis (PCA). It was introduced by Pearson [82] who defines PCA as a linear projection that minimizes the average projection cost. Later, Hotelling [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF] proposed another definition for PCA which reduces the dimension of the data by keeping as much as possible the variation of the data set. In other words, this method aims to find an orthogonal projection of the data set in a low-dimensional linear subspace, such that the variance of the projected data is maximum. This leads to the classical result where the principal axes {u 1 , ..., u d } are the eigenvectors associated with the largest eigenvalues of the empirical covariance matrix S of the data. Several decades after, Tipping and Bishop [START_REF] Tipping | Probabilistic principal component analysis[END_REF] proposed a probabilistic view of PCA by assuming that the observations are independent realizations of a random variable Y ∈ R p which is linked to a latent variable X ∈ R d through the linear relation:

Y = Λ t X + ε.
It is further assume that X ∼ N (µ, I d ) and ε ∼ N (0, σ 2 I p ), such that the marginal distribution of Y is N (Λµ, Λ t Λ + σ 2 I p ). The estimation of the parameters µ, Λ and b by maximum likelihood conduces in particular to estimate Λ by the eigenvectors associated with the largest eigenvalues of the empirical covariance matrix S of the data.

Factor analysis (FA) is an other way to deal with dimension reduction. This approach is as old as PCA since its origins are relative to Spearman [START_REF] Spearman | The proof and measurement of association between two things[END_REF] and there is an important literature on this subject too (see for example [12, chap. 12]). The basic idea of factor analysis is to both reduce the dimensionality of the space and to keep the observed covariance structure of the data. It turns out that the probabilistic PCA (PPCA) model is in fact a particular case of the factor analysis model. Indeed, the FA model makes the same assumption as the PPCA model except regarding the distribution of ε which is assumed to be N (0, Ψ), where Ψ is a diagonal covariance matrix. However, conversely to the PPCA model, the estimation of model parameters by maximum likelihood does not conduce to closed-form estimators.

Regularization

It is also possible to see the curse of dimensionality in clustering as a numerical problem in the inversion of the covariance matrices Σ k in Equation [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF]. From this point of view, a way to tackle the curse of dimensionality is to numerically regularize the estimates of the covariance matrices Σ k before their inversion. As we will see, most of the regularization techniques have been proposed in the supervised classification framework, but they can be easily used for clustering as well. A simple way to regularize the estimation of Σ k is to consider a ridge regularization which adds a positive quantity σ k to the diagonal of the matrix:

Σk = Σk + σ k I p .
Notice that this regularization is often implicitly used in statistical softwares, such as R [START_REF] Venables | Modern Applied Statistics with S[END_REF] for performing a linear discriminant analysis (LDA) where, for instance, the lda() function spheres the data before analysis. A more general regularization has been proposed by Hastie et al. [START_REF] Hastie | Penalized discriminant analysis[END_REF]:

Σk = Σk + σ k Ω,
where Ω is a p × p regularization matrix. This penalization differs from the previous one by the fact that it penalizes also correlations between the predictors. This regularization is used in particular in the supervised classification method penalized discriminant analysis (PDA). Friedman [START_REF] Friedman | Regularized discriminant analysis[END_REF] also proposed, for his famous regularized discriminant analysis (RDA), the following regularization:

Σk (λ, γ) = (1 -γ) Σk (λ) + γ tr( Σk (λ)) p I p ,
where :

Σk (λ) = (1 -λ)(n k -1) Σk + λ(n -K) Σ (1 -λ)(n k -1) + λ(n -K) .
Thus, the parameter γ controls the ridge regularization whereas λ controls the contribution of the estimators Σk and Σ, where Σ estimates the within covariance matrix. Finally, it is also possible to use the Moore-Penrose pseudo-inverse of Σ in Equation ( 9) instead of the usual inverse Σ-1 . The reader can also refer to [START_REF] Mkhadri | Regularization in discriminant analysis: a survey[END_REF] which provides a comprehensive overview of regularization techniques in classification.

More recently, with the exponential growth of high-dimensional data and since many areas in statistical analysis require the inversion or the estimation of covariance matrices, several works were focused on the estimation of large covariance matrices. According to the type of variables, different approaches have been developed. On the one hand, approaches have been proposed to deal with variables having a natural ordering such as longitudinal data or time series (see the works of [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] for example). On the other hand, works were focused on variables for which there is no notion of distances, such as gene expression arrays. For such applications, sparsity is therefore necessary: the methods proposed are mainly based on the use of lasso penalties (see [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] in particular) or on thresholding the covariance matrix [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF][START_REF] Karoui | Operator norm consistent estimation of large dimensional sparse covariance matrices[END_REF]. The reader can refer to the introduction of [START_REF] Bickel | Covariance regularization by thresholding[END_REF] for an overview of existing methods in this field.

Constrained and parsimonious models

A third way to look at the the curse of dimensionality in clustering is to consider it as a problem of over-parameterized modeling. Indeed, as we discussed before in Section 3.2, the Gaussian model turns out to be highly parameterized which naturally yields inference problems in high-dimensional spaces. Consequently, the use of constrained or parsimonious models is another solution to avoid the curse of dimensionality in model-based clustering.

Constrained Gaussian models. A traditional way to reduce the number of free parameters of Gaussian models is to add constraints on the model through their parameters.

Let us recall that the unconstrained Gaussian model (Full-GMM hereafter) is a highly parametrized model and requires the estimation of 20603 parameters when the number of components is K = 4 and the number of variables is p = 100. A first possible constraint for reducing the number of parameters to estimate is to constraint the K covariance matrices to be the same across all mixture components, i.e. Σ k = Σ, ∀k. This model will be denoted to by Com-GMM in the sequel. Notice that this model yields the famous linear discriminant analysis (LDA) [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF] method in the supervised classification case. For the sake of comparison, Table 1 lists the most used constrained models which can be obtained from a Gaussian mixture model with K components in a p-dimensional space. The number of free parameters to estimate, given in the central column, can be decomposed in the number of parameters to estimate for the proportions (K -1), for the means (Kp) and for the covariance matrices (last terms). As one can see in Table 1, it is also possible to assume that the variables are conditionally independent. Such an assumption implies that the covariance matrices are diagonal, i.e. Σ k = diag σ 2 k1 , ..., σ 2 kp with σ 2 k. ∈ R + , and the associated model (Diag-GMM) has a very low number of free parameters. Finally, Sphe-GMM refers to the Gaussian mixture model for which Σ k = σ 2 k I p . Two other intermediate models are presented in Table 1. The Com-Diag-GMM which supposes diagonal common covariances such as Σ k = Σ = diag σ 2 1 , . . . , σ 2 p or the Com-Sphe-GMM model which assumes that the covariance matrices of each class are equal and spherical such that Σ k = Σ = σ 2 I p , for k = 1, ..., K and with σ 2 ∈ R. One can see in Table 1 that, whereas the Full-GMM model is a highly parametrized model, the Sphe-GMM and Com-Sphe-GMM models are Model Name Nb. of parameters

K = 4 p = 100 [λ k D k A k D t k ] VVV (K -1) + Kp + Kp(p + 1)/2 20603 [λD k A k D t k ] EVV* (K -1) + Kp + Kp(p + 1)/2 -(K -1) 20600 [λ k D k AD t k ] VEV (K -1) + Kp + Kp(p + 1)/2 -(K -1)(p -1) 20306 [λD k AD t k ] EEV (K -1) + Kp + Kp(p + 1)/2 -(K -1)p 20303 [λ k DA k D t ] VVE* (K -1) + Kp + p(p + 1)/2 + (K -1)p 5753 [λDA k D t ] EVE* (K -1) + Kp + p(p + 1)/2 + (K -1)(p -1) 5750 [λ k DAD t ] VEE* (K -1) + Kp + p(p + 1)/2 + (K -1) 5456 [λDAD t ] EEE (K -1) + Kp + p(p + 1)/2 5453 [λ k B k ] VVI (K -1) + Kp + Kp 803 [λB k ] EVI (K -1) + Kp + Kp -(K -1) 800 [λ k B] VEI (K -1) + Kp + p + (K -1) 506 [λB] EEI (K -1) + Kp + p 503 [λ k I p ] VII (K -1) + Kp + K 407 [λI p ] EII (K -1) + Kp + 1 404
Table 2: Number of free parameters to estimate for parsimonious Gaussian mixture models with K components and p variables. The models flagged with a star are unfortunately not available in the mclust software (see Section 8.1).

conversely very parsimonious models. They indeed respectively require the estimation of only 407 and 404 parameters when K = 4 and p = 100. Those models however make a strong assumption on the independence of the variables which may be unrealistic in several situations. Finally, the Com-GMM model which requires the estimation of an intermediate number of parameters ( 5453) is known to be a useful model in practical situations. This model is known as well to be efficient for the classification of data sets for which the Gaussian assumption does not hold.

Parsimonious Gaussian models. In a similar spirit than constrained models, Banfield &

Raftery [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF] and Celeux & Govaert [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF] proposed, almost simultaneously, a parameterization of the Gaussian mixture model which yields a family of parsimonious models. To this end, they parametrize the covariance matrices from their eigenvalue decomposition:

Σ k = λ k D k A k D t k ,
where D k is the matrix of eigenvectors which determines the orientation of the cluster, A k is a diagonal matrix proportional to the eigenvalues which explains its shape, and λ k is a scalar which controls its volume. This model is referred to by the [ [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF] and to by VVV in [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF]. By constraining the parameters λ k , D k and A k within and across the groups, 14 different parsimonious models can be enumerated. This family of 14 models are listed in Table 2 in First of all, we can observe that this family of models can be divided in three levels of parsimony, as shown in Table 1. Among the 14 models, 4 models are highly parametrized as the Full-GMM model, 4 models have an intermediate level of parsimony as the Com-GMM model and, finally, 6 models are very parsimonious and are in the same order as the Diag-GMM and Sphe-GMM models. Besides, this reformulation of the covariance matrices can be viewed as a generalization of the constrained models, presented previously. For example, the Com-GMM model is equivalent to the model [λDAD t ]. The model proposed by [START_REF] Murtagh | Fitting straight lines to point patterns[END_REF], which uses the equal shape (λ k = λ, ∀k) and equal volume (A k = A, ∀k), turns out to The reader can refer to [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF] for more details on these models.

λ k D k A k D t k ] model in

Discussion on classical approaches and related works

Firstly, regarding the dimension reduction solution, we would like to caution the reader that reducing the dimension without taking into consideration the clustering goal may be dangerous. Indeed, such a dimension reduction may yield a loss of information which could have been useful for discriminating the groups. In particular, when PCA is used for reducing the data dimensionality, only the components associated with the largest eigenvalues are kept. Such a practice is disproved mathematically and practically by Chang [START_REF] Chang | On using principal component before separating a mixture of two multivariate normal distributions[END_REF] who shows that the first components do not necessary contain more discriminative information than the others. In addition, reducing the dimension of the data may not be a good idea since, as discussed in Section 3, it is easier to discriminate groups in high-dimensional spaces than in lower dimensional spaces, assuming that one can build a good classifier in high-dimensional spaces. With this point of view, subspace clustering methods are good alternatives to dimension reduction approaches. The solution based on regularization does not have the same drawbacks than dimension reduction and can be used with less fear. However, all regularization techniques require the tuning of a parameter which is difficult to tune in the unsupervised context, although this can be done easily in the supervised context using cross validation. Finally, the solution which introduces parsimony in the models is clearly a better solution in the context of model-based clustering since it proposes a trade-off between the perfect modeling and what one can correctly estimate in practice.

We will see in the next sections that recent solutions for high-dimensional clustering are partially based on the idea of parsimonious modeling.

Subspace clustering methods

Conversely to previous solutions, subspace clustering methods exploit the "empty space" phenomenon to ease the discrimination between groups of points. To do so, they model the data in low-dimensional subspaces and introduce some restrictions while keeping all dimensions. Subspace clustering methods can be split into two categories: heuristic and model-based methods. Heuristic methods use algorithms to search for subspaces of [START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF]). high density within the original space. The Clique algorithm [START_REF] Agrawal | Automatic subspace clustering of high-dimensional data for data mining application[END_REF] is one of the heuristic algorithms and remains a reference in this family of methods. A review on heuristic methods is given by [START_REF] Parsons | Subspace clustering for high-dimensional data: a review[END_REF]. On the other hand, model-based subspace clustering methods are mostly related to the factor analysis [START_REF] Rubin | EM algorithms for ML factor analysis[END_REF] model which assumes that the observation space is linked to a latent space through a linear relationship. Figure 1 illustrates the typical modeling of subspace clustering methods.

X Y µ, Λ ε Ψ (a) FA model X Y Z π µ k , Λ k ε Ψ (b) G-MFA model.

Mixture of factor analyzers (MFA)

Mixture of factor analyzers (MFA) [START_REF] Ghahramani | The EM algorithm for factor analyzers[END_REF][START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF] The G-MFA model is an extension of the factor analysis model to a mixture of K factor analyzers. Let {y 1 , . . . , y n } be independent observed realizations of a random vector Y ∈ R p . Let us also consider that Y can be expressed from an unobserved random vector X ∈ R d , named the factor and described in a lower dimensional space of dimension d < p.

Moreover, the unobserved labels {z 1 , . . . , z n } are assumed to be independent unobserved realizations of a random vector Z ∈ {1, . . . , K} where z i = k indicates that y i is generated by the kth factor analyzer. The relationship between these two spaces is finally assumed, conditionally to Z, to be linear:

Y |Z=k = Λ k X + µ k + ε, ( 10
)
where Λ k is a p×d matrix and µ k ∈ R p is the mean vector of the kth factor analyzer. Moreover ε ∈ R p is assumed to be a centered Gaussian noise term with a diagonal covariance matrix Ψ, common to all factors:

ε ∼ N (0, Ψ). (11)
Besides, as in the FA model, the factor X ∈ R d is assumed to be distributed according to a Gaussian density function such as X ∼ N (0, I d ). This implies that the conditional distribution of Y is also Gaussian:

Y |X, Z = k ∼ N (Λ k X + µ k , Ψ). ( 12
)
The marginal density of Y is thus a Gaussian mixture model such as f (y) = K k=1 π k φ(y|θ k ), where π k stands for the mixture proportion and

θ k = {µ k , Λ k Λ t k + Ψ}.
K = 4,d = 3 p = 100 M-MFA S k = Λ k Λ t k + Ψ k (K -1) + Kp + Kd[p -(d -1)/2] + Kp 1991 G-MFA S k = Λ k Λ t k + Ψ (K -1) + Kp + Kd[p -(d -1)/2] + p 1691 MCFA S k = AΩ k A t + Ψ (K -1) + Kd + p + d[p -(d + 1)/2] + Kd(d + 1)/2 433 HFMA S k = V Ω k V t + Ψ (K -1) + (K -1)d + p + d[p -(d -1)/2] + (K -1)d(d + 1)/2 427 MCUFSA S k = A∆ k A t + λI p (K -1) + Kd + 1 + d[p -(d + 1)/2] + Kd 322 A is defined such as A t A = I d , V such as V Ψ -1 V t is
diagonal with decreasing order and ∆ k is a diagonal matrix.

Y |X, Z = k ∼ N (Λ k X + µ k , Ψ k ).
In this case, since there are K covariance matrices of noise to compute in comparison to the G-MFA model, the model complexity increases and takes the following expression for the M-MFA model:

γ M -M F A = (K -1) + Kp + Kd(p -(d -1)/2) + Kp.
Regarding the model inference, Ghahramani & Hinton proposed in [START_REF] Ghahramani | The EM algorithm for factor analyzers[END_REF] an exact EM algorithm for their G-MFA model. In the case of the M-MFA model, McLachlan et al. [START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF] proposed to make use of the alternating expectation-conditional maximization (AECM) [START_REF] Meng | The EM algorithm -an old folk song sung to a fast new tune[END_REF] for parameter estimation. [START_REF] Baek | Mixtures of Factor Analyzers with Common Factor Loadings: Applications to the Clustering and Visualisation of High-Dimensional Data[END_REF] provided an alternative approach which aims to lower the complexity of the MFA model by proposing a more parsimonious modeling. To that end, they re-parametrized the mixture model with restrictions on the means, such as

Extensions of the MFA model

More recently, McLachlan & Baek

µ k = Aρ k , where A is a p × d orthonormal matrix (A t A = I d ) and ρ k is a d-dimensional
vector, and on the covariance matrix S k = AΩ k A t + Ψ, where Ω k is a d × d positive definite symmetric matrix and Ψ a diagonal p × p matrix. This model is referred to by the mixture of factor analyzers with common factor loadings (MCFA) by its authors, as the matrix A is common to the factors. According to the MCFA assumptions, there are only Kd means parameters to estimate instead of Kp in the MFA model. Moreover, since the matrix A is constrained to have orthonormal columns and to be common to all classes, then only pdd(d + 1)/2 loadings are required to estimates it. Finally, according to the restriction on the matrices Ω k , the number of parameters to estimate for these K matrices is Kd(d + 1)/2. Consequently, the complexity of the MCFA model is:

γ M CF A = (K -1) + Kd + p + (pd -d(d + 1)/2)) + Kd(d + 1
)/2 and, for our recurrent practical case, this complexity is equal to 433 which is much more parsimonious than the previous MFA models. Besides, this MCFA approach is a special case of the MFA model but has the main advantage to allow the data to be displayed in a common low-dimensional plot.

The MCFA approach is also a generalization of the works of Yoshida et al. [START_REF] Yoshida | A mixed factor model for dimension reduction and extraction of a group structure in gene expression data[END_REF][START_REF] Yoshida | Array cluster: an analytic tool for clustering, data visualization and model finder on gene expression profiles[END_REF] since they constrained the covariance of the noise term to be spherical (Ψ = λI p ) and the component-covariance matrices of the factors to be diagonal (Ω k = ∆ k = diag(σ 2 k1 , . . . , σ 2 kd )). This approach, called mixtures of common uncorrelated factor spherical-error analyzers

Model name

Cov. structure Nb. of parameters

K = 4,d = 3 p = 100 UUUU -UUU S k = Λ k Λ t k + Ψ k (K -1) + Kp + Kd[p -(d -1)/2] + Kp UUCU - S k = Λ k Λ t k + ω k ∆ k (K -1) + Kp + Kd[p -(d -1)/2] + [1 + K(p -1)] UCUU - S k = Λ k Λ t k + ω k ∆ (K -1) + Kp + Kd[p -(d -1)/2] + [K + (p -1)] UCCU -UCU S k = Λ k Λ t k + Ψ (K -1) + Kp + Kd[p -(d -1)/2] + p UCUC -UUC S k = Λ k Λ t k + ψ k I p (K -1) + Kp + Kd[p -(d -1)/2] + K UCCC -UCC S k = Λ k Λ t k + ψI p (K -1) + Kp + Kd[p -(d -1)/2] + 1 CUUU -CUU S k = ΛΛ t + Ψ k (K -1) + Kp + d[p -(d -1)/2] + Kp CUCU - S k = ΛΛ t + ω∆ k (K -1) + Kp + d[p -(d -1)/2] + [1 + K(p -1)] CCUU - S k = ΛΛ t + ω k ∆ (K -1) + Kp + d[p -(d -1)/2] + [K + (p -1)] 803 CCCU -CCU S k = ΛΛ t + Ψ (K -1) + Kp + d[p -(d -1)/2] + p 800 CCUC -CUC S k = ΛΛ t + ψ k I p (K -1) + Kp + d[p -(d -1)/2] + K 704 CCCC -CCC S k = ΛΛ t + ψI p (K -1) + Kp + d[p -(d -1)/2] + 1 701 where ω k ∈ R + and |∆ k | = 1.
Table 4: Nomenclature of the members of the PGMM and EPGMM families and the corresponding covariance structure.

(MCUFSA), is therefore more parsimonious than MCFA according to the additional assumptions done on the parameters of the MFA model. Finally, Montanari & Viroli [START_REF] Montanari | Heteroscedastic factor mixture analysis[END_REF] presented an approach called heteroscedastic mixture factor model (HMFA) which is very similar to the model described in MCFA. Their model differs from the MCFA approach only on the definition of the common loadings matrix A which does not need to have orthonormal columns. However, to obtain a unique solution for the matrix A, Montanari & Viroli added restrictions on this matrix such as A t Ψ -1 A is diagonal with elements in decreasing order. The links and differences between all these MCFA models are summarized in Table 3 which presents both the covariance structure and the model complexity of each approach.

The inference of both the MCFA and MCUFSA models is done using the AECM algorithm whereas a simple EM algorithm is used for the HMFA model.

Mixture of parsimonious Gaussian mixture models (PGMM)

A general framework for the MFA model was also proposed by McNicholas & Murphy [START_REF] Mcnicholas | Parsimonious Gaussian mixture models[END_REF] which, in particular, includes the previous works of Ghahramani and Hinton and of McLachlan et al. [START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF]. By considering the previous framework, defined by Equations [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] and [START_REF] Bishop | Pattern recognition and machine learning[END_REF], McNicholas & Murphy [START_REF] Mcnicholas | Model-based clustering of microarray expression data via latent gaussian mixture models[END_REF] proposed a family of models known as the expanded parsimonious Gaussian mixture model (EPGMM) family. They decline 12 EPGMM models by either constraining the terms of the covariance matrix to be equal or not, considering an isotropic variance for the noise term, or re-parametrizing the factor analysis covariance structure. The nomenclature of both PGMM and EPGMM is given by Table 4 in which the covariance structure of each model is detailed as well. In particular, the terminology of the PGMM family is as follows: the first letter refers to the loading matrix which is constrained to be common between groups (C..) or not (U..), the second term indicates whether the noise variance is common between factors (.C.) or not (.U.), and the last term refers to the covariance structure which can be either isotropic (..C) or not (..U). Thus, or not (...U). The table also gives the maximum number of free parameters to estimate according to K, p and d for the 12 models. Once again, this number of free parameters to estimate can be decomposed in the number of parameters to estimate for the proportions (K -1), for the means (Kp) and for the covariance matrices (last terms).

X Y Z π k µ k , d k a k1 , ..., a kd k Q k ε b k (a) HD-GMM model X Y Z π k µ k Σ k U ε Ψ (b) DLM model
According to this family of 12 models, the previous approaches developed by [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] Ghahramani | The EM algorithm for factor analyzers[END_REF][START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF][START_REF] Mcnicholas | Parsimonious Gaussian mixture models[END_REF][START_REF] Tipping | Mixtures of Probabilistic Principal Component Analysers[END_REF] then become submodels of the EPGMM approach. For example, by constraining only the noise variance to be isotropic on each class (Ψ k = σ 2 k I p ), which by the way corresponds to the UUC and UUUC models, it produces the famous mixture of prababilistic PCA (Mixt-PPCA) of Tipping & Bishop [START_REF] Tipping | Mixtures of Probabilistic Principal Component Analysers[END_REF]. In the same way, by considering the covariance structure of the UCU and UCCU models such that Ψ k = Ψ and Λ k , then we obtain the mixture of factor analyzers model developed by Ghahramani In the case of the EPGMM models, McNicholas & Murphy [START_REF] Mcnicholas | Model-based clustering of microarray expression data via latent gaussian mixture models[END_REF] proposed for the model inference to make use of the AECM algorithm to speed up the convergence. It is worth to notice that the AECM could be used for inferring most of the MFA-based models.

Mixture of high-dimensional Gaussian mixture models (HD-GMM)

In a slightly different context, Bouveyron et al. [START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF][START_REF] Bouveyron | High Dimensional Discriminant Analysis[END_REF] proposed a family of 28 parsimonious and flexible Gaussian models to deal with high-dimensional data. Conversely to the previous approaches, this family of GMM was directly proposed in both supervised and unsupervised classification contexts. In order to ease the designation of this family, we propose to refer to these Gaussian models for high-dimensional data by the acronym HD-GMM. Bouveyron et al. [START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF] proposed to constraint the GMM model through the eigen-decomposition of the covariance matrix Σ k of the kth group:

Σ k = Q k Λ k Q t k ,
where Q k is a p × p orthogonal matrix which contains the eigenvectors of Σ k and Λ k is a p × p diagonal matrix containing the associated eigenvalues (sorted in decreasing order).

The Such an approach can be viewed in two different ways: on the one hand, these models enable to regularize the models in high-dimension. In particular, by constraining d k such that d k = p -1 for k = 1, ..., K, the proposed approach can be viewed as a generalization of the works of [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF][START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. Indeed, the model [ [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF]. In the same manner, the model [

a kj b k Q k (p -1)] is equivalent to the Full-GMM model or the [λ k D k A k D k ] model in
a kj b k Q (p -1)]
is equivalent to the Diag-GMM and the [a j bQ (p -1)] is also the Com-Diag-GMM. On the other hand, this approach can also be viewed as an extension of the mixture of principal component analyzer (Mixt-PPCA) model [START_REF] Tipping | Mixtures of Probabilistic Principal Component Analysers[END_REF] since it relaxes the equality assumption on d k made in [START_REF] Tipping | Mixtures of Probabilistic Principal Component Analysers[END_REF] and the model [a kj b k Q k d] is therefore equivalent to the Mixt-PPCA model.

In the case of the 16 HD-GMM models listed in Table 5, the inference can be done easily using the EM algorithm since update formula for mixture parameters are closedform. We refer to [START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF] regarding the inference of the other 13 models which require an iterative M step. Finally, the estimation of the intrinsic dimensions d k , k = 1, ..., K, relies on the scree test of Cattell [START_REF] Cattell | The scree test for the number of factors[END_REF] which looks for a break in the eigenvalue scree of the Model name Nb. of parameters p = 100

K = 4, d = 3 [a kj b k Q k d k ] (K -1) + Kp + K k=1 d k [p -(d k + 1)/2] + K k=1 d k + 2K 1599 [a kj bQ k d k ] (K -1) + Kp + K k=1 d k [p -(d k + 1)/2] + K k=1 d k + 1 + K 1596 [a k b k Q k d k ] (K -1) + Kp + K k=1 d k [p -(d k + 1)/2] + 3K 1591 [ab k Q k d k ] (K -1) + Kp + K k=1 d k [p -(d k + 1)/2] + 1 + 2K 1588 [a k bQ k d k ] (K -1) + Kp + K k=1 d k [p -(d k + 1)/2] + 1 + 2K 1588 [abQ k d k ] (K -1) + Kp + K k=1 d k [p -(d k + 1)/2] + 2 + K 1585 [a kj b k Q k d] (K -1) + Kp + Kd[p -(d + 1)/2] + Kd + K + 1 1596 [a j b k Q k d] (K -1) + Kp + Kd[p -(d + 1)/2] + d + K + 1 1587 [a kj bQ k d] (K -1) + Kp + Kd[p -(d + 1)/2] + Kd + 2 1593 [a j bQ k d] (K -1) + Kp + Kd[p -(d + 1)/2] + d + 2 1584 [a k b k Q k d] (K -1) + Kp + Kd[p -(d + 1)/2] + 2K + 1 1588 [ab k Q k d] (K -1) + Kp + Kd[p -(d + 1)/2] + K + 2 1585 [a k bQ k d] (K -1) + Kp + Kd[p -(d + 1)/2] + K + 2 1585 [abQ k d] (K -1) + Kp + Kd[p -(d + 1)/2] + 3 1582 [a j bQd] (K -1) + Kp + d[p -(d + 1)/2] + d + 2 702 [abQd] (K -1) + Kp + d[p -(d + 1)/2] + 3 700
Table 5: Nomenclature for the members of the Hd-GMM family and the number of parameters to estimate. For the numerical example, the intrinsic dimension of the clusters has been fixed to

d k = d = 3, ∀k = 1, . . . , K.
empirical covariance matrix of each group. Let us finally notice that Bouveyron et al. [START_REF] Bouveyron | Intrinsic Dimension Estimation by Maximum Likelihood in Isotropic Probabilistic PCA[END_REF] have demonstrated the surprising result that the maximum likelihood estimator of the intrinsic dimensions d k is asymptotically consistent in the case of the model [

a k b k Q k d k ].

The discriminative latent mixture (DLM) models

Recently, Bouveyron & Brunet [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] proposed a family of mixture models which fit the data into a common and discriminative subspace. This mixture model, called the discriminative latent mixture (DLM) model, differs from the FA-based models in the fact that the latent subspace is common to all groups and is assumed to be the most discriminative subspace of dimension d. This latter feature of the DLM model makes it significantly different from the other FA-based models. Indeed, roughly speaking, the FA-based models choose the latent subspace(s) maximizing the projected variance whereas the DLM model chooses the latent subspace which maximizes the separation between the groups.

Let us nevertheless start with a FA-like modeling. Let Y ∈ R p be the observed random vector and let Z ∈ {1, . . . , K} be once again the unobserved random variable to predict. The DLM model then assumes that Y is linked to a latent random vector X ∈ E through a linear relationship of the form Y = U X + ε, where E ⊂ R p is assumed to be the most discriminative subspace of dimension d ≤ K -1 such that 0 ∈ E, K < p, U is a p × d orthonormal matrix common to the K groups and satisfying U t U = I d , and ε ∼ N (0, Ψ) models the non discriminative information. Besides, within the latent space and conditionally to Z = k, X is assumed to be distributed as

X|Z = k ∼ N (µ k , Σ k )
where µ k ∈ R d and Σ k ∈ R d×d are respectively the mean vector and the covariance matrix Model Nb. of parameters p = 100 and

K = 4 DLM [Σ k β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K 2 (K -1)/2 + K 337 DLM [Σ k β] (K -1) + K(K -1) + (K -1)(p -K/2) + K 2 (K -1)/2 + 1 334 DLM [Σβ k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K(K -1)/2 + K 319 DLM [Σβ] (K -1) + K(K -1) + (K -1)(p -K/2) + K(K -1)/2 + 1 316 DLM [α kj β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K 2 325 DLM [α kj β] (K -1) + K(K -1) + (K -1)(p -K/2) + K(K -1) + 1 322 DLM [α k β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + 2K 317 DLM [α k β] (K -1) + K(K -1) + (K -1)(p -K/2) + K + 1 314 DLM [αj β k ] (K -1) + K(K -1) + (K -1)(p -K/2) + (K -1) + K 316 DLM [αj β] (K -1) + K(K -1) + (K -1)(p -K/2) + (K -1) + 1 313 DLM [αβ k ] (K -1) + K(K -1) + (K -1)(p -K/2) + K + 1 314 DLM [αβ] (K -1) + K(K -1) + (K -1)(p -K/2) + 2 311
Table 6: Number of free parameters to estimate when d = K -1 for the DLM models and some classical models (see text for details).

of the kth group. Given these distribution assumptions, the marginal distribution of Y is once again a mixture of Gaussians, i.e. g(y)

= K k=1 π k φ(y; m k , S k ), where m k = U µ k and S k = U Σ k U t + Ψ. Let W = [U, V ] be the p × p matrix such that W t W = W W t = I p and
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 is block-diagonal: ∆ k = diag (Σ k , βI 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]. Following the classical strategy, several other models can be obtained from the DLM [Σ 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. The list of the 12 different DLM models is given by Table 6 which also provides the number of free parameters to estimate. As we can see, the DLM family yields very parsimonious models and allows, in the same time, to fit into various situations. 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. Let us finally notice that a model similar to the DLM [αβ] model has been considered in [START_REF] Sanguinetti | Dimensionality reduction of clustered datasets[END_REF].

Conversely to most of the MFA-based models, the inference of the DLM models cannot be directly done using the EM algorithm because of the specific features of its latent subspace. To overcome this problem, an estimation procedure, called the Fisher-EM al-gorithm, is also proposed in [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] for estimating both the discriminative subspace 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. This step estimates at iteration q, 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 and conditionally to the posterior probabilities:

Û (q) = max U trace (U t SU ) -1 U t S (q) B U , w.r.t. U t U = I d , ( 13
)
where S stands for the empirical covariance matrix 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 , ( 14
)
denotes the soft between covariance matrix with n

(q) k = n i=1 t 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 and can be set to K -1 in practice.

Two additional procedures are proposed in [START_REF] Bouveyron | On the estimation of the latent discriminative subspace in the Fisher-EM algorithm[END_REF] for the estimation of the latent subspace orientation. The convergence properties of the Fisher-EM algorithm were also studied in [START_REF] Bouveyron | Theoretical and practical considerations on the convergence properties of the Fisher-EM algorithm[END_REF] from both the theoretical and the practical points of view.

Discussion on subspace clustering methods and related works

The subspace clustering methods which have been presented in this section belong to a huge family of Gaussian mixture models and several links exist between these approaches. Indeed, in the last paragraph, we saw that some constrained HD-GMM models are equivalent to traditional parsimonious models such that Com-Diag-GMM or Diag-GMM and consequently to the models proposed by Raftery & Fraley or by Celeux & Govaert. Besides, it also appears that the model [a kj b k Q k d] is equivalent to the Mixt-PPCA model. In the same manner, few models which belong to the EPGMM family of [START_REF] Mcnicholas | Model-based clustering of microarray expression data via latent gaussian mixture models[END_REF] are also included in the HD-GMM family (and vice-versa). In particular, the UCUC model of [START_REF] Mcnicholas | Model-based clustering of microarray expression data via latent gaussian mixture models[END_REF] [START_REF] Baek | Mixtures of Factor Analyzers with Common Factor Loadings: Applications to the Clustering and Visualisation of High-Dimensional Data[END_REF]. Let us highlight that the hypothesis of HD-GMM models are more restrictive than the MFA models since, for example, the subspace of each class is spanned by orthogonal vectors, whereas it is not a necessary condition in MFA, even if such a situation can occur sometimes as in the case of the model UCUC (CUC). The links between the different families of parsimonious models are presented in Figure 4. However, it is worth noting that, despite their efficiency to cluster high-dimensional data, these probabilistic methods based on subspace clustering present several limitations. Indeed, since most models assume that each cluster lives in a specific subspace, the visualization of the clustering results is not easy. Only the MCFA, MCUFSA, HMFA and DLM, for which the latent subspace is assumed to be common, can provide low-dimensional plots of the clustered data.

[a kj b k Q k d k ] [a kj bQ k d k ] [a k b k Q k d k ] [ab k Q k d k ] CCUU [abQ k d k ] UCUU [a j b k Q k d] CUCU [a kj bQ k d] UUCU [a j bQ k d] [λ k B k ] CUUU CUU [a k b k Q k d] [λB k ] CCUC CUC [ab k Q k d] [λ k B] UUCC UUC [a k bQ k d] [λB] CCCC CCC [abQ k d] [λD k AD t k ] EEV CCCU CCU M-MFA [a j bQd] [λ k D k AD t k ] VEV UUUU UUU G-MFA [abQd] [λ k A k ] VVI UUUC UUC Mixt-PPCA [a kj b k Q k d] [λA k ] EVI [a kj b k I p (p -1)] Diag-GMM [λ k A] VEI [a kj b k Q k (p -1)] Full-GMM [λ k D k A k D t k ] VVV [a j bQ (p -1)] Com-GMM [λDAD t ] EEE Sphe-GMM [λ k I p ] VII Com-Sphe-GMM [λI p ] EII
Among the related works, we may cite the robust versions [START_REF] Andrews | Extending mixtures of multivariate t-factor analyzers[END_REF][START_REF] Mclachlan | Extension of the mixture of factor analyzers model to incorporate the multivariate t-distribution[END_REF] of the MFA models which relies on t-distributions and the GMMDR technique [START_REF] Scrucca | Dimension Reduction for Model-Based Clustering[END_REF] which iteratively looks for the subspace containing the most relevant information for clustering. Let us finally notice that some subspace clustering methods have been modified for performing variable selection through ℓ 1 penalization. We will discuss in the next section the sparse extensions of the HMFA and DLM models.

Variable selection for clustering

Conversely to the approaches of the previous section, 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 removing 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. 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 cluster-specific 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.

Variable selection as a model selection problem

The underlying idea of the works of Law et al. [START_REF] Law | Simultaneous Feature Selection and Clustering Using Mixture Models[END_REF], Raftery and Dean [START_REF] Raftery | Variable selection for model-based clustering[END_REF] and Maugis et al. [START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF] is to find the variables which are relevant for the clustering task. The determination of the role of each variable is in particular apprehended in [START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF][START_REF] Raftery | Variable selection for model-based clustering[END_REF] and Celeux & Govaert [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF], based on a specific decomposition of the mixture component variance matrix (see Section 4 for more details).

In the Raftery & Dean's approach, the authors define two different sets of variables: S which denotes the set of relevant variables and S c which is the set containing the irrelevant variables. An interesting aspect of their approach is that they do not assume that the irrelevant variables are independent of the clustering variables conversely to Law et al. [START_REF] Law | Simultaneous Feature Selection and Clustering Using Mixture Models[END_REF]. In particular, they define the irrelevant variables as those which are independent of the clustering but which remain dependent of the set of relevant variables according to a linear relationship. The models in competition are compared with the integrated loglikelihood via a BIC approximation. Thus, the selected model maximizes the following quantity:

K, m, Ŝ = arg max (K,m,r,ℓ,V ) BIC clust (y S |K, m) + BIC reg (y Sc |y S ) , ( 15
)
where K is the number of clusters and m ∈ M is a model which belongs to the family of parsimonious models available in the mclust [START_REF] Fraley | MCLUST: Software for Model-Based Cluster Analysis[END_REF] software. Note that the quantity to be maximized in expression [START_REF] Bouveyron | On the estimation of the latent discriminative subspace in the Fisher-EM algorithm[END_REF], can be decomposed into two parts: the first term corresponds to the Gaussian mixture model of K components on the subset of relevant variables S, whereas the second one is relative to the regression of irrelevant variables in S c on the set of all clustering variables in S. However, the dependence assumption which defines the irrelevant set of variables according to all the relevant ones remains debatable. Indeed, on the one hand, considering only the case where the irrelevant variables are independent on both the clustering and the relevant partition, as it was considered in the work of Law et al. [START_REF] Law | Simultaneous Feature Selection and Clustering Using Mixture Models[END_REF], seems to be unrealistic. On the other hand, considering that all the irrelevant variables depends on the relevant variables by a linear relationship seems to be as well a strong hypothesis which may be not valid in certain practical cases. An other limitation of the Raftery & Dean's procedure is linked to their variable selection algorithm. Indeed, they proposed in [START_REF] Raftery | Variable selection for model-based clustering[END_REF] a forward-stepwise algorithm which considers only few variables at the beginning and which prevents from taking into account the block interactions between variables.

To overcome these limitations, Maugis et al. [START_REF] Celeux | Letter to the editor[END_REF][START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF][START_REF] Maugis | Variable selection in modelbased clustering: A general variable role modeling[END_REF] K, m, r, ĥ, V = arg max

(K,m,r,h,V ) BIC clust (y S |K, m) + BIC reg (y U |r, y R ) + BIC ind (y W |h) , ( 16
)
where V = (S, R, U, W) stands for the variable partition. The first term of this expression, called BIC clust , corresponds to the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF] for a Gaussian mixture of K components on the relevant subset of variables S. The model m belongs here to a collection of 28 parsimonious models which are available in the Mixmod software [START_REF] Biernacki | Model-based cluster and discriminant analysis with the mixmod software[END_REF] and include the GMM family introduced by Celeux & Govaert [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF]. The second term denoted by BIC reg , is linked to the BIC criterion for a linear regression of the irrelevant variables U on a subset of clustering variables R. Note that the index r stands for the structure of the covariance matrix which can be assumed to be spherical, diagonal or non-constraint. Finally, the last term depicts the BIC criterion for a Gaussian density on the variable subset W independent of the clustering variables. This Gaussian marginal distribution is characterized by a variance matrix σ which is constrained to be either diagonal or spherical and is specified by the index h in the expression [START_REF] Bouveyron | Discriminative variable selection for clustering with the sparse Fisher-EM algorithm[END_REF].

The identifiability of the SRUW model and the consistency of their variable selection problem has been studied in [START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF]. Regarding the implementation, they propose an algorithm based on a backward-stepwise selection. It implies that all the variables are considered at the beginning of the procedure and only a block of variables is either included or excluded of the clustering relevant set of features at each iteration. Such an approach enables them to take into account variable block interactions, if they exist. Then a second algorithm is executed to select both the model and the number of components for the mixture model.

Variable selection by likelihood penalization

An other way to combine variable selection and clustering is to penalize the clustering criteria in order to yield sparsity in the features. This technique has been used, in particular, by penalizing the log-likelihood function to optimize. A general form for the penalized log-likelihood function is:

L p (θ) = ℓ(θ) -p λ (θ) (17)
where ℓ(θ) stands for the log-likelihood function and p λ (θ) is the penalty function. In GMM context, Pan & Shen [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] proposed a penalized log-likelihood criterion by assuming a Gaussian mixture model with common diagonal covariance matrices, meaning that ∀k ∈ {1, . . . , K},

Σ k = Σ = diag(σ 2 1 , . . . , σ 2 j , . . . , σ 2 p ) where σ 2 j ∈ R.
The penalty function is focused on the means of K clusters (m 1k , . . . , m pk , ∀k ∈ {1, . . . , K}) and has the following form:

p λ (θ) = λ 1 K k=1 p j=1 |m kj | , ( 18
)
where m kj denotes the mean of the jth variable in the component k and λ 1 an hyperparameter which stands for the desired level of sparsity. Thus, since the observations are standardized, if the means of a variable j on each component are equal i.e. m 1j = • • • = m Kj = 0, then this variable is irrelevant and can be removed from the clustering variables.

Therefore, a variable selection is realized when some m kj 's can be shrunken toward 0. This situation occurs for an ℓ 1 penalty term large enough. In the same spirit, Wang & Zhou [START_REF] Wang | Variable selection for model-based high dimensional clustering and its application to microarray data[END_REF] proposed two other penalty terms. The first one is based on ℓ ∞ -norm:

p λ (θ) = λ ∞ p j=1 max k∈{1,...,K} |m kj | , ( 19
)
which has the advantage to incorporate group information. Thus, this penalty tends to shrink all the m kj 's toward 0 as soon as the jth variable is non informative. However, such a penalty tends to shrink the m kj 's in the same magnitude and thus does no take into account the situation where a variable is different from 0 on only one component. To that end, Wang & Zhou proposed a second penalty function based on hierarchical penalties.

These three penalized log-likelihood functions have been developed with the restriction of the diagonal common covariance matrix of each cluster in the mixture model. Xie et al. [START_REF] Xie | Penalized model-based clustering with cluster-specific diagonal covariance matrices and grouped variables[END_REF] extended the model of Pan & Shen [START_REF] Pan | Penalized model-based clustering with application to variable selection[END_REF] by relaxing the equality constraint on the covariance matrices. Indeed, they proposed an approach dealing with the case of clusterspecific diagonal covariance matrices (∀k ∈ {1, . . . , K}, S k = diag(σ 2 k1 , . . . , σ 2 kp )) leading to the following penalty function:

p λ (θ) = λ 1 K k=1 p j=1 |m kj | + λ 2 K k=1 p j=1 σ 2 kj -1 . ( 20
)
In this case, a second regularized term is added and holds on the variance of the variable j of the kth component σ 2 kj which can be shrunk towards 0. As previously, the hyperparameters λ 1 and λ 2 are selected through a modified BIC criterion, which takes into account the level of sparsity in the model complexity term. Finally, Zhan et al. [START_REF] Zhang | A flexible and efficient algorithm for regularized fisher discriminant analysis[END_REF] recently proposed a penalization in the case of Com-GMM model (S k = S, ∀k ∈ {1, . . . K}):

p λ (θ) = λ 1 K k=1 p j=1 |m kj | + λ 2 p ℓ=1 p j=1 |C jℓ | ,
where {C jℓ } p j,ℓ=1 are the elements of the inverse covariance matrix S -1 . In the same work, Zhan et al. also proposed an estimation procedure to deal with the n < p case.

Variable selection by penalization of the loadings

An alternative approach for selecting the relevant variables through penalization is to directly apply the lasso penalty on the loading matrix of a MFA-based model. This has been achieved in particular in [START_REF] Bouveyron | Discriminative variable selection for clustering with the sparse Fisher-EM algorithm[END_REF][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]. In the case of the MFA model, Galimberti et al. [START_REF] Galimberti | Penalized factor mixture analysis for variable selection in clustered data[END_REF] introduced an ℓ 1 -penalty on the factor loadings in the log-likelihood function such as:

p λ (θ) = λ 2 d ℓ=1 p j=1 |b ℓj | (21)
where b ℓj stands for the factor loadings. In a very recent work, Xie et al. [START_REF] Xie | Penalized mixtures of factor analyzers with application to clustering high-dimensional microarray data[END_REF] proposed a penalized MFA approach from the model introduced by Gharamani & Hinton where the covariance matrix of the noise term is diagonal and common to all factors. The penalty function has the following form:

p λ (θ) = λ 1 K k=1 p j=1 |m kj | + λ 2 K k=1 p j=1 b 2 kj , ( 22
)
where b kj stands for the factor loading of the kth factor. As in the previous approaches, the first term based on the ℓ 1 -norm is used to shrink the means m kj to be exactly equal to 0 while the second term serves as a grouped variable penalty. Indeed, this last penalty aims to shrink the estimates of factor loadings b kj which are close to 0 to be exactly equal to 0. Consequently, if a variable has a common mean equal to 0, a common variance on each factor across the clusters and is independent with all other cluster such as b kj = 0 ∀k, then this variable is irrelevant and does not contribute in the clustering task.

In the context of Fisher-EM, the direct penalization of the loading matrix U makes particularly sense since it is not estimated by likelihood maximization. The matrix U is indeed estimated in the F-step of the Fisher-EM algorithm by maximizing the Fisher criterion conditionally to the current partition of the data. To achieve the penalization of U, two solutions are proposed in [START_REF] Bouveyron | Discriminative variable selection for clustering with the sparse Fisher-EM algorithm[END_REF]. The first solution is a two stage approach which first estimate U , at each iteration, with the F-step and then looks for its best sparse approximation as follows:

min U X (q)t -Y t U 2 F + λ d j=1 u j 1 ,
where u j is the jth column vector of U , X (q) = Û (q)t Y and F refers to the Frobenius norm. The solution of this penalized regression problem can be computed through the LARS algorithm [START_REF] Efron | Least angle regression[END_REF]. The second solution directly recasts the maximization of the Fisher criterion as a regression problem and provides a sparse loading matrix by solving the lasso problem associated to this regression problem. Let us define, conditionally to the posterior probabilities t (q) ik , the matrices 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 and H (q) B = 1 √ n n (q) 1 (m (q) 1 -ȳ), . . . , n (q) K (m (q) K -ȳ) ∈ R p×K , where n (q) k = n i=1 t (q) ik and m (q) k = 1 n n i=1 t (q)
ik y i . According to these definitions, the matrices

H (q)
W and

H (q) B sat- isfy H (q) W H (q)t W = S (q) W and H (q) B H (q)t B = S (q)
B where S (q) W and S (q) B are respectively the soft within and between covariance matrices computed at iteration q. Then, the best sparse approximation of the solution of the Fisher criterion 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. However, solving this lasso problem is not direct in this case and requires the use of an iterative algorithm. Regarding the implementation details, it is proposed in [START_REF] Bouveyron | Discriminative variable selection for clustering with the sparse Fisher-EM algorithm[END_REF] to initialize the sparseFEM algorithm with the result of the Fisher-EM algorithm and to determine the value of λ by model selection through a modified BIC criterion.

Discussion on variable selection methods and related works

As we have seen, two main approaches are nowadays used for achieving variable selection in the context of model-based clustering. On the one hand, the approaches based on model selection present the clear advantage of being fully automatic since they automatically select the number of variables to retain. Their current implementations, based on forward-backward strategies, can however prevent their application to very highdimensional data. On the other hand, the approaches based on lasso-type penalizations are usually fast and work well even on very high-dimensional data. However, the selection of the sparsity parameter λ is still an open issue. Solutions on model selection criteria are indeed relatively unstable and require further studies.

Among the related works, we shall cite the recent work of Witten & Tibshirani [START_REF] Witten | A framework for feature selection in clustering[END_REF], even though it is actually not a model-based approach. They proposed a general nonprobabilistic framework for variable selection in clustering, based on a penalized criterion which governs both variable selection and clustering. Several clustering methods can be reformulated from such a criterion and Witten & Tibshirani apply in particular their criterion to the k-means and hierarchical clustering methods. Let us also notice that a recent work of Galimberti & Soffritti [START_REF] Galimberti | Using conditional independence for parsimonious model-based Gaussian clustering[END_REF] proposes an intermediate solution between parsimonious modeling and variable selection by constraining the group covariance matrices to be block-diagonal.

Softwares for model-based clustering

This paragraph aims to list the existing softwares which provide model-based clustering methods for high-dimensional data.

Sofwares with parsimonious models

A first and basic function for the Matlab software, called emgm [START_REF] Mo | emgm: EM algorithm for Gaussian mixture model[END_REF], allows to fit data with a Gaussian mixture model, without allowing to impose any restrictions on the covariance matrices. The statlearn toolbox [START_REF] Bouchard | The statlearn toolbox: statistical learning tools for Matlab[END_REF] for Matlab contains conversely many supervised and unsupervised learning algorithms among which several model-based methods for clustering. For the R software, the reference package mclust [START_REF] Fraley | MCLUST: Software for Model-Based Cluster Analysis[END_REF] allows to fit a mixture of Gaussians with different volumes and shapes of the group covariance matrices. The modeling choices can be controlled by the user (see next paragraph) or chosen according to the BIC criterion. An EM algorithm fits the parameters of the considered model and allows agglomerative hierarchical clustering based on maximum likelihood. In order to deal with large data sets, the package pmclust [START_REF] Chen | Parallel Model-Based Clustering[END_REF], available for the R software, performs a parallel version of the EM algorithm for mixtures of Gaussians and this allows to fit ultra large data. The standalone EMMIX software [START_REF] Mclachlan | The EMMIX software[END_REF][START_REF] Mclachlan | The emmix software for the fitting of mixtures of normal t-components[END_REF] provides a robust approach based on tdistributions and on the ECM algorithm [START_REF] Mclachlan | Robust cluster analysis via mixtures of multivariate t-distributions[END_REF]. With this t-mixture model-based approach, the normal distribution for each component in the mixture is embedded in a wider class of elliptically symmetric distributions. The MIXMOD [START_REF] Biernacki | Model-based cluster and discriminant analysis with the mixmod software[END_REF] software, written in C++ and with interfaces with the Matlab, Scilab and R softwares, is a general tool for model-based supervised and unsupervised classification. The main asset of this software remains in the fact that it deals with both quantitative and qualitative data. To that end, different models are proposed: 14 parsimonious Gaussian models are available for quantitative data and 5 multinomial mixture models can be used for qualitative data. This software allows as well to deal with high-dimensional data. Moreover, likelihood maximization can be done using the CEM or SEM algorithms in addition to the traditional EM algorithm. The model and the number of clusters can be chosen by different criteria according to the required task: the BIC, ICL [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] and NEC [START_REF] Bouchard | Model selection in supervised classification[END_REF] criteria can be used for clustering and the cross-validation is also available in the supervised classification context. Finally, an R version of MIXMOD has been released recently and the package is called Rmixmod.

Subspace clustering methods

Different softwares and R packages provide subspace clustering methods in the Gaussian mixture model context. Firstly, the package HDclassif [START_REF] Bergé | HDclassif : an R Package for Model-Based Clustering and Discriminant Analysis of High-Dimensional Data[END_REF] for the R software provides routines for model-based clustering and classification of high-dimensional data. In particular, the hddc function implements the subspace clustering method of [START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF]. Among the 28 original HD-GMM models, 14 models can be used in the function: 12 models with group-specific orientation matrices and 2 models with common covariance matrices (see Table 5). The hddc function allows to choose the most appropriate model for the data at hand according to the BIC criterion and, for each group, the intrinsic dimension is chosen according the Cattell scree-test or BIC. The EM algorithm is used for fitting the model parameters. Let us notice that HDclassif provides also the hdda function for executing a discriminant analysis on high-dimensional data. Both classifiers are also available in the MIXMOD software. Amongst the mixture of factor analyzers, the R package pgmm [START_REF] Mcnicholas | Parsimonious Gaussian mixture models[END_REF][START_REF] Mcnicholas | Model-based clustering of microarray expression data via latent gaussian mixture models[END_REF] provides an implementation of model-based clustering and model-based classification using several parsimonious Gaussian mixture models. The Fortran EMMIX-MFA software [START_REF] Mclachlan | The EMMIX-MFA software[END_REF] implements the mixture of factor analyzers [START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF]. The software allows the user to choose the number of factor analyzers to fit and the shape of covariance matrices to use. Finally, the mfma function [START_REF] Viroli | The mmfa function for the R software[END_REF] for the R software extends and combines the mixture of factor analyzers [START_REF] Mclachlan | Modelling high-dimensional data by mixtures of factor analyzers[END_REF] and the factor mixture analyzers [START_REF] Montanari | Heteroscedastic factor mixture analysis[END_REF]. More precisely, it provides a unified class of dimensionally reduced mixture models and offers an interesting tool for modeling non-Gaussian latent variables. Some of these subspace clustering approaches model the data in a low and common subspace which allows the data to be displayed in a low-dimensional plot. The mcfa function [START_REF] Mclachlan | The mcfa function for the R software[END_REF] for R, which implements the method of [START_REF] Baek | Mixtures of Factor Analyzers with Common Factor Loadings: Applications to the Clustering and Visualisation of High-Dimensional Data[END_REF], is one of them. It enables model-based density estimation to be undertaken for high-dimensional data, where the number of observations is not very large relative to their dimension. Still for the R software, the hmfa function [START_REF] Viroli | The hmfa function for the R software[END_REF] implements the approach proposed in [START_REF] Montanari | Heteroscedastic factor mixture analysis[END_REF], also based on the mixture of factor analyzers. Finally, the FisherEM package for the R software implements the eponymous algorithm [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF] and provides an efficient tool to model and cluster the data at hand in a discriminative and low-dimensional latent subspace.

Variable selection and clustering in a GMM context

The clustvarsel package for the R software implements the approach proposed by [START_REF] Raftery | Variable selection for model-based clustering[END_REF]. The package allows the user to simultaneously partition the data in the GMM context and select the relevant variables using an approximate Bayes factor. The C++ software selvarclust [START_REF] Maugis | The selvarclust software[END_REF] implements the extension proposed by [START_REF] Maugis | Variable selection for Clustering with Gaussian Mixture Models[END_REF][START_REF] Maugis | Variable selection in modelbased clustering: A general variable role modeling[END_REF] which embeds the Raftery & Dean's approach. This software allows to cluster data where individuals are described by quantitative block variables. It returns a data clustering and the variable partition. A sparse version of the FisherEM algorithm is also available in FisherEM package for R. Performing a selection of discriminative variables within the FisherEM algorithm is made possible by the use of the option method='sparse' in the fem function. Finally, the R package bclust [START_REF] Partovi Nia | High-Dimensional Bayesian Clustering with Variable Selection: The R Package bclust[END_REF] is useful for clustering high-dimensional continuous data. The package uses a parametric spike-and-slab Bayesian model to down-weight the effect of noise variables and to quantify the importance of each variable in agglomerative clustering. Unlike HDclassif, the bclust package implements a Bayesian approach to the clustering, with priors for model parameters and for the allocation of subjects to groups. The model and its priors are chosen so that the marginal posterior is analytically tractable, providing a fast algorithm.

Numerical experiments

This section now considers some specific R packages and data sets in order to illustrate the practical use of model-based clustering techniques on high-dimensional data.

Why it is important not to reduce the dimension before the clustering ?

Before to move further, we would like to convince the reader that reducing the data dimension without taking into account the clustering goal conduces to suboptimal results. Indeed, as we discussed earlier, the disconnection between the dimension reduction and the clustering steps can lead to a loss of information which could have been discriminative for the clustering. Such a problem is illustrated in this paragraph on the Crabs data set [START_REF] Campbell | A multivariate study of variation in two species of rock crabs of genus Leptograpsus[END_REF] because, even though it is not high-dimensional, it is known to be difficult to cluster and allows nice visualizations. The data consist in 200 Leptograpsus variegatus crabs, collected at Fremantle, Australia. Each individual is described by 5 morphological measurements: frontal lobe size, rear with, carapace length and width, and finally body depth. The data set is made of 4 groups of 50 crabs each, grouped together according to their color (blue or orange) and their sex (female or male). The main difficulty in clustering these data comes from the existing linear correlation between the 5 variables which implies a "size effect" when a PCA is executed. Therefore, if a PCA is run without any special care, then the clustering task can lead to poor results. Let us first illustrate such a remark by running in first a PCA on the Crabs data set before clustering the data using the mclust package for R. The scree plot is presented on Figure 5a and shows that the 2 first principal components seem to be sufficient to represent the data. Indeed, almost 99% of the total variance is explained by the 2 first principal axes. However, the quality of the partition obtained by running Mclust on the 2 first principal components is disappointing since almost 51% of the data are misclassified with regard to the true partition. In particular, as we can observe in Table 7a, only 2 different groups can be distinguished in this case and the gender is not taking into account. This can be explained by the task of dimension reduction and the size effect. We then switched the 1st principal component with the third one since, as one can observe in Figure 5b, the 2nd and 3rd principal components appear to be discriminant for the 4 groups.

res <-Mclust(X[,2:3],G=4)

The resulted partition is now clearly better since the correct classification rate reaches 82% of agreement with the know partition for this data set. In particular, the orange males (OM) and females (OF) form now 2 distinct groups but we have still a mixture of genders for the blue crabs. This is illustrated in the clustered data obtained with 3 subspace clustering methods: MCFA, HMFA and Fish-erEM. We use for this experiment the chironomus data set which is made of 148 larvae belonging to 3 different species: cloacalis, frommeri and staegeri. For each larva, 17 characters of the larval head capsule have been measured. A detailed description of this data set can be found in [START_REF] Montanari | Heteroscedastic factor mixture analysis[END_REF]. In first, for each method, the default settings have been used. In particular, the FisherEM algorithm fixes itself the intrinsic dimension of the discriminative subspace to d = K -1. For the MCFA and HMFA approaches, the number of factor analyzers is determined using the BIC criterion. The first results shows that, except for the FisherEM algorithm for which the discriminative space is 2-dimensional, the BIC criterion selects only one factor analyzer for both MCFA and HMFA approaches. In order to visualize the fitted subspaces of 1 dimension, we have projected the data in the 2 first principal components of PCA for the MCFA and HMFA approaches. They are represented by a straight line in Figures 7a and7b. As we can observe, 2 of the 3 groups seem to be difficult to separate linearly. This difficulty is illustrated with HMFA which merges 2 groups. This might be explained by the fact that only one factor analyzer is not sufficient to discriminate 3 groups. Conversely, the subspaces built by MCFA and HMFA are such that the variance of the projected data is maximum.

Variable selection through penalized approaches

We focus now on variable selection for clustering with penalized approaches. We consider here the USPS358 data set which is a subset of the USPS data set from the UCI repository. The original data set is made of 7291 images divided in 10 classes corresponding to handwritten 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 = 1756) corresponding to the digits 3, 5 and 8 which are the most difficult to discriminate. Figure 9 shows a subset of images from the USPS358 data set and Figures 10a-10d show the group means and the global mean of the data. As we can observe on these images, there are 3 The level of sparsity is: 0.2 with a penalized bic: 2997.332

Figure 10e shows the absolute values of the loadings of the two discriminative axes estimated by sparseFEM. The weights assigned by sparseFEM to each feature are represented by gray levels: lighter is the pixel, weaker is the absolute value of the loadings.

We can observe in Figure 10e that the subset of selected pixels is small compared to the average number of pixels used to draw the digits 3, 5 and 8. Only 15 pixels are selected by sparseFEM among the 256 original ones. Furthermore, the selected pixels appear to be relevant since, for instance, the darker pixel on the bottom right corner of Figure 10e discriminates the digit 8 from the digits 3 and 5. The variables 2, 9 and 13 are selected and the resulting partition on these 3 variables for the zoo data set is resumed by Table 9a. We can observe that 4 among 7 known groups are correctly recognized. For the selvarclust algorithm, we used the C++ code 9b that 5 groups are well recognized.

Variable selection as a model selection problem

Conclusion

This article has therefore presented a review of existing solutions for model-based clustering of high-dimensional data. We hope that the reader is now convinced that it is better to use parsimonious models, subspace clustering methods or variable selection methods designed for clustering instead of preprocessing the data with dimension reduction. The short software review and the few practical examples may also help the practitioner in applying recent model-based techniques to his own data.

  which the first column stands for the model names used by Celeux & Govaert and the second one corresponds to the nomenclature used by Raftery & Fraley.

  be equivalent to the model [λD k AD t k ]. It is worth to notice that the work of Celeux & Govaert widens the family of parsimonious models since they add unusual models which allow different volumes for the clusters such as the [λ k DAD t ], [λ k DA k D t ] and [λ k D k AD t k ] models. Furthermore, by assuming that the covariance matrix Σ k are diagonal matrices, Celeux & Govaert proposed a new parametrization Σ k = λ k B k where |B k | = 1. Such a parametrization leads to 4 additional models listed at the bottom of Table 2. Finally, by considering the spherical shape, it leads to 2 other models: the [λ k I p ] and [λI p ] models.

Figure 1 :

 1 Figure1: Illustration of the modeling of subspace clustering methods in the case of two groups (illustration from[START_REF] Bouveyron | High-Dimensional Data Clustering[END_REF]).

Figure 2 :

 2 Figure 2: Graphical summary of factor analysis (FA) model (a) and mixture of factor analyzers (G-MFA) model of Ghahramani and Hinton (b).

  may be considered at the earliest subspace clustering method which both clusters the data and reduces locally the dimensionality of each cluster. The MFA model differs from the FA model by the fact that it allows to have different local factor models, in different regions of the input space, whereas the standard FA assumes a common factor model. The MFA model was firstly introduced by Ghahramani & Hinton [43], and then extended by McLachlan et al. [65]. In this first paragraph, the original work of Ghahramani & Hinton on MFA is introduced before to develop the work of Baek & McLachlan on this subject and other recent works [4, 74, 103, 104]. To make a distinction between both approaches, we will refer to the model of Ghahramani & Hinton by G-MFA and to the model of McLachlan et al. by M-MFA.

  Figures 2a and 2b summarize respectively the FA and G-MFA models. The complexity of the G-MFA model can be computed according to the number of parameters to estimate. Since the G-MFA model is in a Gaussian mixture model of K components, there are (K -1) parameters for the proportions and Kp for the means. Moreover, Kd(p -(d -1)/2) + p parameters are required to estimate the component covariance matrices, since the covariances matrices are defined in a factor representation such as S k = Λ k Λ t k + Ψ. The model complexity is then γ G-M F A = (K -1) + Kp + Kd(p -(d -1)/2) + p and, by considering the practical case where p = 100, K = 4 and d = 3, then 1691 parameters have to be estimated for this G-MFA model. This approach, introduced by Ghahramani & Hinton, was generalized a few years later by McLachlan et al. [65] who removed in particular the constraint on the variance of the noise. Therefore, the conditional distribution of the noise term becomes ε|Z = k ∼ N (0, Ψ k ) where Ψ k stands for the diagonal covariance matrix of the cluster k. The new Model name Cov. structure Nb. of parameters

Table 3 :

 3 Nomenclature of the MFA models developed by Ghahramani and Hinton (G-MFA), MacLachlan et al. (M-MFA), and several MCFA models with their corresponding covariance structure. conditional distribution of Y is then:

Figure 3 :

 3 Figure 3: Graphical representations of the HD-GMM model [a kj b k Q k d k ] (a) and of the DLM [Σ k β] model (b).

  & Hinton. The UUUU model is as well equivalent to the MFA model proposed by McLachlan et al. in [65]. Finally, by parameterizing the covariance structure as Ψ k = ω k ∆ k , where ∆ k is a diagonal matrix and |∆ k | = 1, McNicholas & Murphy proposed 4 additional models to their previous work [69].

  key idea of the work of Bouveyron et al. is to reparametrize the matrix Λ k , such as Σ k has only d k + 1 different eigenvalues: Λ k = diag (a k1 , . . . , a kd k , b k , . . . , b k ) , where the d k first values a k1 , . . . , a kd k parametrize the variance in the group-specific subspace and the pd k last terms, the b k 's model the variance of the noise and d k < p.With this parametrization, these parsimonious models assume that, conditionally to the groups, the noise variance of each cluster k is isotropic and is contained in a subspace which is orthogonal to the subspace of the kth group. Following the classical parsimony strategy, the authors proposed a family of parsimonious models from a very general model,the model [a kj b k Q k d k ],to very simple models. Table5lists the 16 HD-GMM models with closed-form estimators and reports their complexity for comparison purposes. Once again, the first quantity (K -1) + Kp stands for the number of parameters for the means and the mixture proportions of K clusters. Then, there are K k=1 d k [p(d k + 1)/2] loadings to estimate for the K orientation matrices Q k . Finally, the last terms represent the number of free parameters for the covariance matrices in the latent and in the noise subspaces of K clusters and their intrinsic dimension.

  corresponds to the HD-GMM model [a kj b k Q k d]. Moreover, the EPGMM family proposed by McNicholas & Murphy includes individual works on MFA such as the works of Ghahramani [43], Tipping and Bishop [92], McLachlan [65], McNicholas and Murphy [70] and Baek et al.

Figure 4 :

 4 Figure 4: Links between MFA-bsed models for subspace clustering and some parsimonious models.

  as a model selection problem in the GMM context. In particular, Raftery & Dean and Maugis et al. consider a collection of parsimonious and interpretable models, developed by Banfield & Raftery[START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF]

  relax such restrictions and propose a more general variable role modeling. They define two subsets of variables: on the one hand, the relevant ones, which are grouped in S and, on the other hand, its complementary S c , which is formed by the irrelevant variables. Maugis et al. consider two types of behaviors among these irrelevant variables: a subset U of irrelevant variables which can be explained by a linear regression from a subset R of the clustering variables and a subset W of irrelevant variables which are totally independent of all relevant variables. Such a variable partition allows to both consider the approaches developed by Law et al.[START_REF] Law | Simultaneous Feature Selection and Clustering Using Mixture Models[END_REF] and by Raftery & Dean[START_REF] Raftery | Variable selection for model-based clustering[END_REF]. It is referred to by the model collection SRUW. From this characterization, the authors also recast the variable selection problem into a model selection problem through an approximation of the integrated log-likelihood. Then the selected model satisfies:

Figure 5 :

 5 Figure 5: Scree plot obtained by executing a PCA on the Crabs'data set (a) and scatterplots of projected data in the 3 first principal components (b).

Figure 7 :

 7 Figure 7: Projection of the clustered chironomus data in the 2 first principal components of PCA. The fitted subspace is represented by parallel lines centered on the means of each cluster for the MCFA and HMFA approaches.

  (scale(Y),3,model='all') The best model is: AkjBk and the BIC is: 59.33106 library(mcfa) resultat = pairlist() bic = rep(NA,6)for (i in 1:6){resultat[[i]] <-mcfa(scale(Y),g=3,q=i,maxinit=5) for (i in 1:6){bic[i]<-resultat[[i]]$BICval} max_bic = max(bic,na.rm=T) for(i in 1:6){ if(max_bic==bic[i]) {res2=resultat[[i]] for(i in 1:6){ if(max_bic==bic[i]) {q = i} }The number of factor analyzers is: 1 and the BIC value is: 4694.725 source(L1mfa.r)resultat = pairlist() bic = rep(NA,6) for (i in 1:6){resultat[[i]] <-L1mfa.em(scale(Y),K,i) for (i in 1:6){bic[i]<-resultat[[i]]$bic} max_bic = max(bic,na.rm=T) for(i in 1:6){ if(max_bic==bic[i]) {res3=resultat[[i]] for(i in 1:6){ if(max_bic==bic[i]) {q = i} }The number of factor analyzers is: 1 and the BIC value is: 4698.284

Figure 8 :

 8 Figure 8: Projection of the clustered chironomus data in the 2-dimensional subspaces fitted by MCFA, HMFA and FisherEM.

Figure 9 :

 9 Figure 9: Samples from the usps358 data set.

Finally, another way

  to combine variable selection and clustering is to consider variable selection as a model selection problem in the mixture model context. We consider here the clustvarsel [83] and selvarclust [58] algorithms and applied them on the zoo data set coming from the UCI repository. This data set is made of 7 families of 101 animals characterized by 16 attributes. For the clustvarsel algorithm, a package for the R software is available: library(clustvarsel) res <-clustvarsel(Y,G=7) colnames(res$sel.var) The variables selected are: 2 9 13 result <-Mclust(res$sel.var,1:3)

( a )

 a Group mean of digit 3 (b) Group mean of digit 5 (c) Group mean of digit 8 (d) Overall mean of the data set (e) Variable selection by sparseFEM

Figure 10 :Table 9 :

 109 Figure 10: Pixels averaging on the 1756 images (a) and discriminative pixels selected by the sparseFEM algorithm (b).

  10 evaluations of this function which was unfordable problem in the 60's (it remains a difficult problem nowadays). Although the term "curse of dimensionality" used by Bellman is of course rather pessimistic, the paragraph of the preface in which the term first appeared is in fact

more optimistic:

All this [the problems linked to high dimension] may be subsumed under the heading « the curse of dimensionality ». Since this is a curse, [...], there is no need to feel discouraged about the possibility of obtaining significant results despite it.

Table 1 :

 1 Number of free parameters to estimate for constrained Gaussian models.

	Model	Nb. of parameters	K = 4 and p = 100
	Full-GMM	(K -1) + Kp + Kp(p + 1)/2 20603
	Com-GMM	(K -1) + Kp + p(p + 1)/2	5453
	Diag-GMM	(K -1) + Kp + Kp	803
	Com-Diag-GMM (K -1) + Kp + p	503
	Sphe-GMM	(K -1) + Kp + K	407
	Com-Sphe-GMM (K -1) + Kp + 1	404

Table 7b .

 7b The dimension reduction task done independently of the clustering task is in fact prejudicial. Let us now use a subspace clustering method which allows to simultaneously reduce the dimension and find a partition of data. The clustering of the Crabs data set is done with the pgmmEM function from the pgmm package (all the default settings are kept):We now consider subspace clustering methods with group-specific subspaces. In this experiment, the HDclassif and PGMM packages are used to cluster the Diagnostic Wisconsin Breast Cancer data set coming from the UCI repository[START_REF] Frank | UCI Machine Learning Repository[END_REF]. This data set consists of 30 quantitative features which are computed from a digitized image of a fine needle aspirate of a breast mass. They describe characteristics of the cell nuclei present in the image. The number of instances is 569 whose are divided on two groups of nuclei: 357 are benign and 212 are malignant. Of course the class label is not used for the clustering task but, for both methods, the number of clusters have been fixed to 2 and the BIC criterion is used for selecting the model and the intrinsic dimension of each group.The model selected by hddc is the [a kj b k Q k d] model, which is formally equivalent to the Mixt-PPCA model, with all intrinsic dimensions equal to 11. Table8ashows the contingency table which allows to compare the known partition for the data with the one returned by the hddc function. One can observed that in this case the two groups are well-identified. Regarding the PGMM approach, the model selected by the BIC criterion is the CUU model. This choice of model implies that the loading matrix is common to the groups but the error variance remains different. Let us notice that this model is much more parsimonious than the one chosen by hddc since the intrinsic dimension for both groups is equal to 2 here. Consequently, the complexity of the CUU model is 180 whereas the one of the [a kj b k Q k d] model reaches 614 parameters to estimate. The counterpart of such a parsimony is perhaps a little degradation of the quality of the partition obtained by the PGMM approach on this data set, as we can observe in Table8b.

	8.2. Subspace clustering methods with group-specific subspaces
		label			label
	cluster Malign Benign	cluster Malign Benign
	1	172	16	1	145	52
	2	40	341	2	67	305
	(a) HDDC		(b) PGMM
	library(HDclassif)
	res1 <-hddc(as.matrix(Y),K=K,model="ALL",d='BIC')
	SELECTED: model AKJBKQKD with 2 clusters, BIC=16579.28.
	res1$d
	Intrinsic dimensions of the classes:
	1 2
	dim: 11 11
	library(pgmm)
	res2 <-pgmmEM(Y,Gmin=2,Gmax=2)
	Based on k-means starting values, the best model (BIC) for the range of
	factors and components used is a CUU model with q = 2 and G = 2. The BIC for
	this model is 23779.74.
	8.2.1. Subspace clustering methods with common subspaces
	library(pgmm)
	res <-pgmmEM(crabs[,5:8],Gmin=4,Gmax=4)

Several subspace clustering approaches model the data in a low-dimensional and common subspace. This allows in particular the clustered data to be displayed on a lowdimensional plot. In this experiment, we aim to compare the different visualizations of

Table 8 :

 8 Contingency tables for the HDDC and PGMM algorithms for the subset of the breast cancer data.

Acknowledgments

The authors would like to greatly thank the associate editor and two referees for their helpful remarks and comments on the manuscript.

UCC model with q = 1 and G = 4

The BIC for this model is -2589.95.

The quality of the resulted partition is once again really improved and only 9.5% of Crabs are now misclassified: while the orange females (OF) and blue males (BM) are now almost totally well classified, the orange males (OM) still overlap with the orange females (OF), as resumed in Table 7c. As the estimated dimension of the intrinsic subspace of each class is estimated by the pgmmEM function to be equal to 1, it allows an easy representation of the group subspaces using line segments. Figure 6a shows these group subspaces as line segments with the projection of the clustered data on the 2 first principal components.

This practical example has therefore make the demonstration that, when groups are localized in different subspaces, a global dimension reduction is not adapted for the clustering task. A efficient way to overcome such a problem is to use subspace clustering methods which model and cluster the data in group-specific low-dimensional subspaces.