Charles Bouveyron

Mathieu Fauvel

Stéphane Girard

Kernel discriminant analysis and clustering with parsimonious Gaussian process models

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Classification is an important and useful statistical tool in all scientific fields where decisions have to be made. Depending on the availability of a learning data set, two main situations may happen: supervised classification (also known as discriminant analysis) and unsupervised classification (also known as clustering). Discriminant analysis aims to build a classifier (or a decision rule) able to assign an observation x in an arbitrary space E with unknown class membership to one of k known classes C 1 , ..., C k . For building this supervised classifier, a learning dataset {(x 1 , z 1 ), ..., (x n , z n )} is used, where the observation x ℓ ∈ E and z ℓ ∈ {1, ..., k} indicates the class belonging of the observation x ℓ .

In a slightly different context, clustering aims to directly partition an incomplete dataset {x 1 , ..., x n } into k homogeneous groups without any other information, i.e., assign to each observation x ℓ ∈ E its group membership z ℓ ∈ {1, ..., k}. Several intermediate situations exist, such as semi-supervised or weakly-supervised classifications [START_REF] Chapelle | Semi-Supervised Learning[END_REF], but they will not be considered here.

Since the pioneer work of Fisher [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF], a huge number of supervised and unsupervised classification methods have been proposed in order to deal with different types of data. Indeed, there exist a wide variety of data such as quantitative, categorical and binary data but also texts, functions, sequences, images and more recently networks. As a practical example, biologists are frequently interested in classifying biological sequences (DNA sequences, protein sequences), natural language expressions (abstracts, gene mentioning), networks (gene interactions, gene co-expression), images (cell imaging, tissue classification) or structured data (gene structures, patient information). The observation space E can be therefore R p if quantitative data are considered, L 2 ([0, 1]) if functional data are considered (time series for example) or A p , where A is a finite alphabet, if the data at hand are categorical (DNA sequences for example). Furthermore, the data to classify can be a mixture of different data types: categorical and quantitative data or categorical and network data for instance.

Classification methods can be split into two main families: generative and discriminative techniques. On the one hand, generative or model-based techniques model the data of each class with a probability distribution and deduce the classification rule from this modeling. On the other hand, discriminative techniques directly build the classification rule from the learning dataset. Among the discriminative methods, kernel methods and Gaussian process classification are the most popular. Section 2 briefly reviews these techniques.

In this work, we propose to adapt model-based methods for the classification of any kind of data by working in a feature space of high or even infinite dimensional space. To this end, we propose a family of parsimonious Gaussian process models which allow to build, from a finite sample, a model-based classifier in a infinite dimensional space. It will be demonstrated that the building of the classifier can be directly done from the observation space through the so-called "kernel trick". The proposed classification method will be thus able to classify data of various types (categorical data, mixed data, functional data, networks, etc). The methodology is as well extended to the unsupervised classification case (clustering).

The paper is organized as follows. Section 2 reviews generative and discriminative techniques for classification. Section 3 then presents the context of our study and introduces the family of parsimonious Gaussian process models. The inference aspects are addressed in Section 4. It is also demonstrated in this section that the proposed method can work directly from the observation space through a kernel. Section 5 is dedicated to some special cases and to the extension to the unsupervised framework through the derivation of an EM algorithm. Experimental comparisons with state-of-the-art kernel methods are presented in Section 6 as well as applications of the proposed methodologies to various types of data including functional, categorical and mixed data. Some concluding remarks are given in Section 7 and proofs are postponed to the appendix.

Related work

This section briefly reviews model-based and discriminative techniques for classification.

Model-based classification

Model-based discriminant analysis assumes that {x 1 , ..., x n } are independent realizations of a random vector X on E and that the class conditional distribution of X is parametric, i.e. f (x|z = i) = f i (x; θ i ).

When E = R p , among the possible parametric distributions for f i , the Gaussian distribution is often preferred and, in this case, the marginal distribution of X is therefore a mixture of Gaussian distributions:

f (x) = k i=1 π i φ(x; µ i , Σ i ),
where φ is the Gaussian density, π i is the prior probability of the ith class, µ i is the mean of the ith class and Σ i is its covariance matrix. The optimal decision rule is called the maximum a posteriori (MAP) rule which assigns a new observation x to the class with the largest posterior probability. Introducing the classification function D i (x) = log |Σ i | + (xµ i ) t Σ -1 i (xµ i ) -2 log(π i ), which can be rewritten as:

D i (x) = p j=1 1 λ ij < x -µ i , q ij > 2 R p + p j=1 log(λ ij ) -2 log(π i ), (1)
where q ij and λ ij are respectively the jth eigenvector and eigenvalue of Σ i , it can be easily shown that the MAP rule reduces to finding the label i ∈ {1, . . . , k} for which D i (x) is the smallest. Estimation of model parameters is usually done by maximum likelihood. This method is known as the quadratic discriminant analysis (QDA), and, under the additional assumption that Σ i = Σ for all i ∈ {1, . . . , k}, it corresponds to the linear discriminant analysis (LDA). A detailed overview on this topic can be found in [START_REF] Mclachlan | Discriminant Analysis and Statistical Pattern Recognition[END_REF]. Recent extensions allowing to deal with high-dimensional data include [START_REF] Bouveyron | Simultaneous model-based clustering and visualization in the Fisher discriminative subspace[END_REF]Bouveyron et al., 2007b,a;[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] Montanari | Heteroscedastic factor mixture analysis[END_REF][START_REF] Murphy | Variable selection and updating in model-based discriminant analysis for high dimensional data with food authenticity applications[END_REF] and an extensive review of those techniques is given in [START_REF] Bouveyron | Model-based clustering of high-dimensional data: A review[END_REF]. Although model-based classification is usually enjoyed for its multiple advantages, model-based discriminant analysis methods is mostly limited to quantitative data. Some extensions exist to handle categorical data using multinomial [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF] or Dirichlet [START_REF] Bouguila | Novel mixtures based on the Dirichlet distribution: application to data and image classification[END_REF] distributions for instance. In addition, even in the case of quantitative data, the Gaussian assumption may not be well-suited for the data at hand. Recent works focused on different distributions such as the skew normal [START_REF] Lin | Robust mixture modeling using the skew t distribution[END_REF], asymmetric Laplace [START_REF] Franczak | 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] Lin | Robust mixture modeling using multivariate skew t distribution[END_REF][START_REF] Forbes | A new family of multivariate heavy-tailed distributions with variable marginal amounts of tailweight: application to robust clustering[END_REF].

Kernel methods for classification

Among discriminative classification methods, kernel methods [START_REF] Hofmann | Kernel methods in machine learning[END_REF] are probably the most popular and overcome some of the shortcomings of generative techniques. Kernel methods are non-parametric techniques and can be applied to any data for which a kernel function can be defined. A kernel K : E × E → R is a positive definite function such as every evaluation can be written as K(x i , x j ) =< ϕ(x i ), ϕ(x j ) > H , with x i , x j ∈ E, ϕ a mapping function (called the feature map), H a finite or infinite dimensional reproducing kernel Hilbert space (the feature space) and < •, • > H the dot product in H. An advantage of using kernels is the possibility of computing the dot product in the feature space from the original input space without explicitly knowing ϕ (kernel trick) [START_REF] Hofmann | Kernel methods in machine learning[END_REF]. Turning conventional learning algorithms into kernel learning algorithms can be done if the algorithms operate on the data only in terms of dot product. For instance, kernels can be defined on strings [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF][START_REF] Cuturi | The context-tree kernel for strings[END_REF], graphs [START_REF] Smola | Kernels and regularization on graphs[END_REF][START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF], trees (Schölkopf et al., 2004, Chap. 5) and vector-valued functions [START_REF] Caponnetto | Universal multi-task kernels[END_REF][START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF][START_REF] Kadri | Multiple Operator-valued Kernel Learning[END_REF]. Many conventional linear algorithms have been turned to non-linear algorithms thanks to kernels [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]: KPCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF], KFD [START_REF] Mika | Fisher discriminant analysis with kernels[END_REF], KGMM [START_REF] Dundar | Toward an optimal supervised classifier for the analysis of hyperspectral data[END_REF], mixture of Gamma distributions [START_REF] Murua | Kernel-based Mixture Models for Classification[END_REF] and others. Let us however highlight that those kernel versions often involve the inversion of a kernel matrix, i.e., a n × n matrix estimated with only n samples.

Usually, the kernel matrix is ill-conditioned and regularization is needed, while sometimes a simplified model is required too. Thus, it may limit the effectiveness of the kernel version. In addition, and conversely to model-based techniques, classification results provided by kernel methods are difficult to interpret which would be useful in many application domains.

Classification with Gaussian processes

A recent approach in discriminative classification is to make use of Gaussian processes (GP) in a Bayesian view of logistic regression. In this context, Gaussian processes are used as prior distribution for the regression function which links the predictive variable (Z) to the descriptive vector (X). The conditional

distribution of Z is then given through the relation logit (p(Z = 1|X = x)) = f (x)
where the prior distribution on the regression function f is f ∼ GP (0, Σ) with Σ a given covariance function. In the regression context, the inference is usually tractable since the Gaussian process prior is combined with a Gaussian likelihood. Unfortunately, in the classification case, the discrete nature of the predictive variable makes the exact inference infeasible. To overcome this difficulty, several approximation techniques have been proposed lately, such as the Laplace approximation [START_REF] Kuss | Assessing approximate inference for binary Gaussian process classification[END_REF] or through the expectation-propagation algorithm (Minka, 2001). Let us notice that the latter approximation is usually more accurate in term of classification performance than the former in practical situations (Rasmussen and Williams, 2006b). Furthermore, a limiting point of Gaussian process classification is, once again, the need to numerically invert a n × n covariance matrix. Thus, the inference of a Gaussian process predictor, even in a feasible case, is computationally demanding (O(n 3 ) where n is the number of learning observations). We refer to (Rasmussen and Williams, 2006b, Chapter 3) for an extensive overview of Bayesian classification with Gaussian processes.

Classification with parsimonious Gaussian process models

In this section, we first define the context of our approach and exhibit the associated computational problems. Then, a parsimonious parameterization of Gaussian processes is proposed in order to overcome the highlighted computational issues. Before to move forward and in order to avoid any misunderstanding, we would like to emphasize that we do not consider in this work Gaussian processes in a Bayesian regression framework. In our approach, the Gaussian processes are used as conditional distributions of a latent process and not as prior distributions for Bayesian logistic regression.

Classification with Gaussian processes

Let us consider a learning set {(x 1 , z 1 ), ..., (x n , z n )} where {x 1 , ..., x n } ⊂ E are assumed to be independent realizations of a, possibly non-quantitative and non-Gaussian, random variable X. Let us assume that there exists a non-linear mapping ϕ such that the latent process Y = ϕ(X) is, conditionally on Z = i, a Gaussian process on J ⊂ R with mean µ i and continuous covariance function Σ i . Hence, the use of the transformation ϕ allows to embed a possible non-quantitative vector into a quantitative space, as in kernel methods. More specifically, one has µ i (t) = E(Y (t)|Z = i) and [START_REF] Shorack | Empirical Processes with Applications to Statistics[END_REF] that, for all i = 1, . . . , k, there exist positive eigenvalues (sorted in decreasing order) {λ ij } j≥1 , together with eigenvector functions {q ij (.)} j≥1 continuous on J, such that

Σ i (s, t) = E(Y (s)Y (t)|Z = i) -µ i (t)µ i (s). It is then well-known
Σ i (s, t) = ∞ j=1 λ ij q ij (s)q ij (t),
where the series is uniformly convergent on J × J. Moreover, the eigenvector functions are orthonormal in L 2 (J) for the dot product < f, g > L2 = ´J f (t)g(t)dt. It is then easily seen, that, for all r ≥ 1 and i ∈ {1, . . . , k}, the random vector on R r defined by {< Y, q ij > L2 } j=1,...,r is, conditionally on Z = i, Gaussian with mean {< µ i , q ij >} j=1,...,r and covariance matrix diag(λ i1 , . . . , λ ir ). To classify a new observation x, we therefore propose to apply the Gaussian classification function (1) to ϕ(x):

D i (ϕ(x)) = r j=1 1 λ ij < ϕ(x) -µ i , q ij > 2 L2 + r j=1 log(λ ij ) -2 log(π i ).
From a theoretical point of view, if the Gaussian process is non degenerated, one should use r = +∞.

In practice, r has to be large in order not to loose to much information on the Gaussian process.

Unfortunately, in this case, the above quantities cannot be estimated from a finite sample set. Indeed, only a part of the classification function can be actually computed from a finite sample set:

D i (ϕ(x)) = ri j=1 1 λ ij < ϕ(x) -µ i , q ij > 2 L2 + ri j=1 log(λ ij ) -2 log(π i ) computable quantity + r j=ri+1 1 λ ij < ϕ(x) -µ i , q ij > 2 L2 + r j=ri+1 log(λ ij ) non computable quantity
, where r i = min(n i , r) and n i = Card (C i ) is the cardinality of the ith class. Consequently, the Gaussian model cannot be used directly in the feature space to classify data if r > n i for i = 1, ..., k.

A parsimonious Gaussian process model

To overcome the computation problem highlighted above, it is proposed here to use in the feature space a parsimonious model for the Gaussian process modeling each class. Following the idea of (Bouveyron et al., 2007a), we constrain the eigen-decomposition of the Gaussian processes as follows.

Definition 1. A parsimonious Gaussian process (pgp) model is a Gaussian process Y for which, conditionally to Z = i, the eigen-decomposition of its covariance operator Σ i is such that:

(A1) it exists a dimension r < +∞ such that λ ij = 0 for j ≥ r and for all i = 1, ..., k, (A2) it exists a dimension d i < min (r, n i ) such that λ ij = λ for d i < j < r and for all i = 1, ..., k.

It is worth noticing that r can be as large as it is desired, as long it is finite, and in particular r can be much larger than n i , for any i = 1, ..., k. From a practical point of view, this modeling can be viewed as assuming that the data of each class live in a specific subspace of the feature space. The variance of the actual data of the ith group is modeled by the parameters λ i1 , ..., λ idi and the variance of the noise is modeled by λ. This assumption amounts to supposing that the noise is homoscedastic and its variance is common to all the classes. The dimension d i can be also considered as the intrinsic dimension of the latent subspace of the ith group in the feature space. This model is referred to as

D i (ϕ(x)) = di j=1 1 λ ij - 1 λ < ϕ(x) -µ i , q ij > 2 L2 + 1 λ ||ϕ(x) -µ i || 2 L2 + di j=1 log(λ ij ) + (d max -d i ) log(λ) -2 log(π i ) + γ, ( 2
)
where γ is a constant term which does not depend on the index i of the class.

At this point, it is important to notice that the classification function D i depends only on the eigenvectors associated with the d i largest eigenvalues of Σ i . This estimation is now possible due to the inequality d i < n i for i = 1, ..., k. Furthermore, the computation of the classification function does not depend any more on the parameter r. As shown in Section 4, it is possible to reformulate the classification function such that it does not depend either on the mapping function ϕ.

Submodels of the parsimonious model

By fixing some parameters to be common within or between classes, it is possible to obtain particular models which correspond to different regularizations. Table 1 presents the 8 additional models which can be obtained by constraining the parameters of model M 0 . For instance, fixing the dimensions d i to be common between the classes yields the model M 1 . Similarly, fixing the first d i eigenvalues to be common within each class, we obtain the more restricted model M 2 . It is also possible to constrain the first d i eigenvalues to be common between the classes (models M 4 and M 7 ), and within and between the classes (models M 5 , M 6 and M 8 ). This family of 9 parsimonious models should allow the proposed classification method to fit into various situations.

Model inference and classification with a kernel

This section focuses on the inference of the parsimonious models proposed above and on the classification of new observations through a kernel. Model inference is only presented for the model M 0 since inference for the other parsimonious models is similar. Estimation of intrinsic dimensions and visualization in the feature subspaces are also discussed.

Estimation of model parameters

In the model-based classification context, parameters are usually estimated by their empirical counterparts [START_REF] Mclachlan | Discriminant Analysis and Statistical Pattern Recognition[END_REF] which leads, in the present case, to estimate the proportions π i by πi = n i /n and the mean function

µ i by μi (t) = 1 n i xj ∈Ci ϕ(x j )(t)
, where n i is the number of observations in the ith class. Regarding the covariance operator, the eigenvalue λ ij and the eigenvector q ij are respectively estimated by the jth largest eigenvalue λij and its associated eigenvector function qij of the empirical covariance operator Σi :

Σi (s, t) = 1 n i x ℓ ∈Ci ϕ(x ℓ )(s)ϕ(x ℓ )(t) -μi (s)μ i (t).
Finally, the estimator of λ is:

λ = 1 k i=1 πi (r -d i ) k i=1 πi   trace( Σi ) - di j=1 λij   . ( 3
)
Using the plug-in method, the estimated classification function Di can be written as follows:

Di (ϕ(x)) = di j=1 1 λij - 1 λ < ϕ(x) -μi , qij > 2 L2 + 1 λ ||ϕ(x) -μi || 2 L2 + di j=1 log( λij ) + (d max -d i ) log( λ) -2 log(π i ). (4)
However, as we can see, the estimated classification function Di still depends on the function ϕ and therefore requires computations in the feature space. However, since all these computations involve dot products, it will be shown in Proposition 2 that the estimated classification function can be computed without explicit knowledge of ϕ through a kernel function.

Estimation of the classification function through a kernel

Kernel methods are all based on the so-called "kernel trick" which allows the computation of the classifier in the observation space through a kernel K. Let us therefore introduce the kernel

K : E × E → R defined as K(x, y) =< ϕ(x), ϕ(y) > L2 and ρ i : E × E → R defined as ρ i (x, y) =< ϕ(x) -µ i , ϕ(y) - µ i > L2 .
In the following, it is shown that the classification function D i only involves ρ i which can be computed using K:

ρ i (x, y) = 1 n 2 i x ℓ ,x ℓ ′ ∈Ci < ϕ(x) -ϕ(x ℓ ), ϕ(y) -ϕ(x ℓ ′ ) > L2 (5) = K(x, y) - 1 n i x ℓ ∈Ci (K(x ℓ , y) + K(x, x ℓ )) + 1 n 2 i x ℓ ,x ℓ ′ ∈Ci K(x ℓ , x ℓ ′ ). (6) Kernels K(x, y) r i Linear < x, y > L2 min(n i , p) RBF exp - x-y 2 L 2 2σ 2 n i Polynomial (< x, y > L2 +1) q min n i , p+q p Table 2: Dimension r i for several classical kernels.
For each class C i , let us introduce the n i × n i symmetric matrix M i defined by:

(M i ) ℓ,ℓ ′ = ρ i (x ℓ , x ℓ ′ ) n i .
With these notations, we have the following result (proof is postponed to the appendix).

Proposition 2. For i = 1, . . . , k, the estimated classification function can be computed, in the case of the model M 0 , as follows:

Di (ϕ(x)) = 1 n i di j=1 1 λij 1 λij - 1 λ x ℓ ∈Ci β ijℓ ρ i (x, x ℓ ) 2 + 1 λ ρ i (x, x) + di j=1 log( λij ) + (d max -d i ) log( λ) -2 log(π i ),
where, for j = 1, . . . , d i , β ij is the normed eigenvector associated to the jth largest eigenvalue λij of

M i and λ = 1/ k i=1 πi (r i -d i ) × k i=1 πi trace(M i ) - di j=1 λij .
It thus appears that each new sample point x can be assigned to the class C i with the smallest value of the classification function without knowledge of ϕ. Moreover, this operation only requires the estimation of the d i largest eigenvalues and eigenvectors of M i . In practice, this computation is achieved thanks to the implicitly restarted Arnoldi method (IRAM) [START_REF] Lehoucq | Deflation techniques for an implicitly restarted arnoldi iteration[END_REF] which is proved to be a stable procedure. The methodology based on Proposition 2 is referred to as pgpDA in the sequel. In practice, the value of r i depends on the chosen kernel (see Table 2 for examples).

Intrinsic dimension estimation

The estimation of the intrinsic dimension of a dataset is a difficult problem which occurs frequently in data analysis, such as in principal component analysis. A classical solution in PCA is to look for a break in the eigenvalue scree of the covariance matrix. This strategy relies on the fact that the jth eigenvalue of the covariance matrix corresponds to the fraction of the full variance carried by the jth eigenvector of this matrix. Since, in our case, the class conditional matrix M i shares with the empirical covariance operator of the associated class its largest eigenvalues, we propose to use a similar strategy based on the eigenvalue scree of the matrices M i to estimate d i , i = 1, ..., k. More precisely, we propose to make use of the scree-test of Cattell [START_REF] Cattell | The scree test for the number of factors[END_REF] for estimating the class specific dimension d i , i = 1, ..., k.

For each class, the selected dimension is the one for which the subsequent eigenvalues differences are smaller than a threshold which can be tuned by cross-validation for instance.

Visualization in the feature subspaces

An interesting advantage of the approach is to allow the visualization of the data in subspaces of the feature space. Indeed, even though the chosen mapping function is associated with a space of very high or infinite dimension, the proposed methodology models and classifies the data in low-dimensional subspaces of the feature space. It is therefore possible to visualize the projection of the mapped data on the feature subspaces of each class using Equation ( 10) of the appendix. The projection of ϕ(x) on the jth axis of the class C i is therefore given by:

P ij (ϕ(x)) :=< ϕ(x) -μi , qij >= 1 n i λij x ℓ ∈Ci β ijℓ ρ i (x, x ℓ ).
Thus, even if the observations are non quantitative, it is possible to visualize their projections in the feature subspaces of the classes which are quantitative spaces. Notice that the obtained visualizations of the data are of course dependent on the chosen kernel, such as in kernel PCA for instance.

Links with existing works and extension to clustering

The methodology proposed in Section 4 is made very general by the large choice for the mapping function ϕ(x). We focus in this section on two specific choices for ϕ(x) for which the direct calculation of the classification rule is possible. We also exhibit the links with related methodologies. Finally, an extension to unsupervised classification is considered through the use of an EM algorithm.

The case of the linear kernel for quantitative data

In the case of quantitative data, E = R p and one can choose ϕ(x) = x associated to the standard scalar product which gives rise to the linear kernel K(x, y) = x t y. In such a framework, the estimated classification function can be simplified as follows (proof is provided in the appendix):

Proposition 3. If E = R p and K(x, y) = x t y then, for i = 1, . . . , k, the estimated classification function reduces to Di (x) = di j=1 1 λij - 1 λ qt ij (x -μi ) 2 + 1 λ ||x -μi || 2 R p + di j=1 log( λij ) + (d max -d i ) log( λ) -2 log(π i ).
where μi is the empirical mean of the class C i , qij is the eigenvector of the empirical covariance matrix Σi associated to the jth largest eigenvalue λij and λ is given by ( 3).

It appears that the estimated classification function reduces to the one of the high-dimensional discriminant analysis (HDDA) method (Bouveyron et al., 2007a) with the model [a ij bQ i d] which has constraints similar to M 0 . Therefore, the methodology proposed in this work partially encompasses the method HDDA. Let us however highlight that the proposed methodology is not just a kernelized version of HDDA. One can remark that it is clearly not obvious to obtain Proposition 2 from Proposition 3. Indeed, substituting dot products by kernels in the classification function of the latter proposition does not yield the classification function of Proposition 2.

The case of functional data

Let us consider now functional data observed in E = L 2 ([0, 1]). Let (b j ) j≥1 be a basis of L 2 ([0, 1]) and F = R L where L is a given integer. For all ℓ = 1, . . . , L, the projection of a function x on the jth basis function is computed as

γ j (x) = ˆ1 0 x(t)b j (t)dt
and γ(x) := (γ j (x)) j=1,...,L . Let B the L × L Gram matrix associated to the basis:

B jℓ = ˆ1 0 b j (t)b ℓ (t)dt,
and consider the associated scalar product defined by < u, v >= u t Bv for all u, v ∈ R L . One can then choose ϕ(x) = B -1 γ(x) and K(x, y) = γ(x) t B -1 γ(y) leading to a simple estimated classification function.

Proposition 4. Let E = L 2 ([0, 1]) and K(x, y) = γ(x) t B -1 γ(y). Introduce, for i = 1, . . . , k, the L × L covariance matrix of the γ(x j ) when x j ∈ C i : Σi = 1 n i x ℓ ∈Ci (γ(x ℓ ) -γi )(γ(x ℓ ) -γi ) t where γi = 1 n i xj ∈Ci γ(x j )
Then, for i = 1, . . . , k, the estimated classification function reduces to

Di (ϕ(x)) = di j=1 1 λij - 1 λ qt ij (γ(x) -γi ) 2 + 1 λ (γ(x) -γi ) t B -1 (γ(x) -γi ) + di j=1 log( λij ) + (d max -d i ) log( λ) -2 log(π i ),
where qij and λij are respectively the jth normed eigenvector and eigenvalue of the matrix B -1 Σi and λ is given by (3).

Proof of the above result is provided in the appendix. Remark that B -1

Σi coincides with the matrix of interest in functional PCA (Ramsay and Silverman, 2005, Chap. 8.4) and that, if the basis is orthogonal, then B is the identity matrix. Notice that the proposed method therefore encompasses as well the model proposed in [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF] for the clustering of functional data.

Links with existing works

In addition to the links with existing methods exhibited above, it is also possible to establish some links with other generative techniques proposed in the literature. Kernelized versions of Gaussian mixture models were introduced in the unsupervised [START_REF] Wang | Kernel trick embedded Gaussian mixture model[END_REF] and supervised [START_REF] Xu | A novel kernel-based maximum a posteriori classification method[END_REF] classification contexts. These methods do not estimate the smallest eigenvalues of the kernel matrix from the data but use instead a constant threshold. Furthermore, they numerically regularize the estimated covariance matrices before to invert them. [START_REF] Mika | Fisher discriminant analysis with kernels[END_REF] have also proposed kernel

Fisher discriminant (KFDA) as a non-linear version of FDA which only relies on kernel evaluations.

KFDA shares as well the objectives of pgpDA but assumes a common covariance matrix for all classes and uses a ridge regularization for the inversion of this covariance matrix. Recently, pseudo-inverse and ridge regularization have been proposed to define a kernel quadratic classifier where classes have their own covariance matrices [START_REF] Pekalska | Kernel discriminant analysis for positive definite and indefinite kernels[END_REF]. Similarly, a common covariance assumption has been proposed in [START_REF] Dundar | Toward an optimal supervised classifier for the analysis of hyperspectral data[END_REF] for the construction of a LDA classifier in the kernel feature space. Yet, small eigenvalues are thresholded in the computations. Finally, in the unsupervised context, kernel k-means (Shawe-Taylor and Cristianini, 2004) can be also viewed as using a constraint version of our pgp models. Indeed, it uses model M 8 with d i = d = r and π i = π = 1/K (all mixture parameters are common). To summarize, all kernel methods have to face the inversion of a kernel matrix (denoted by M i in our framework) which is ill-conditioned or singular. The above mentioned techniques to overcome this problem include eigenvalues thresholding, shape constraints on the matrices, shrinkage methods (such as ridge regularization) and arbitrary combinations of these three methods. Here, our approach relies on the two assumptions (A1) and (A2) which benefit from clear interpretations, see Section 3.2. Moreover, no additional regularization has to be made, since the model parameters can be estimated in a stable way, see Section 4.1. Finally, the possible visualization of the data in subspaces of the feature space (Section 4.4) and the ability to deal with both supervised and non-supervised frameworks are additional arguments in favor of our approach.

Extension to unsupervised classification

In model-based classification, the unsupervised and supervised cases mainly differ in the manner to estimate the parameters of the model. The clustering task aims to form k homogeneous groups from a set of n observations {x 1 , ..., x n } without any prior information about their group memberships. Since the labels are not available, it is not possible in this case to directly estimate the model parameters. In such a context, the expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] allows to both estimate the model parameters and predict the class memberships of the observations at hand. In the case of the parsimonious model M 0 introduced above, the EM algorithm takes the following form:

The E step This first step reduces, at iteration s, to the computation of t

(s) ij = E(Z j = i|x j , θ (s-1)
), for j = 1, . . . , n and i = 1, . . . , k, conditionally on the current value of the model parameter θ (s-1) .

We have the following proposition (proof is postponed to the appendix): (s-1) ) are given by:

Proposition 5. The conditional expectations t (s) ij = E(Z j = i|x j , θ
t (s) ij = 1/ k ℓ=1 exp D (s-1) i (ϕ(x j )) -D (s-1) ℓ (ϕ(x j )) , ( 7
)
where

D (s-1) i (ϕ(x)) = 1 n (s-1) i di j=1 1 λ(s-1) ij 1 λ(s-1) ij - 1 λ(s-1) n ℓ=1 β ijℓ √ t iℓ ρ (s-1) i (x, x ℓ ) 2 + 1 λ(q-1) ρ (s-1) i (x, x) + di j=1 log( λ(s-1) ij ) + (d max -d i ) log( λ(s-1) ) -2 log(π (s-1) i
).

is the Gaussian classification function associated with the model parameters estimated in the M step at iteration s -1 and n

(s-1) i = n j=1 t (s-1) ij .
The M step This second step estimates the model parameters conditionally on the posterior probabilities t (q) ij computed in the previous step. In practice, this step reduces to update the estimate of model parameters according to the following formula:

• mixture proportions are estimated by π(q

) i = n (q) i /n where n (q) i = n j=1 t (q) ij ,
• parameters λ ij , λ, β ij and d i are estimated at iteration q using the formula given in Proposition 2 but where the matrix M i is now a n × n matrix, recomputed at each iteration q, and such that, for i = 1, ..., k and ℓ, ℓ ′ = 1, ..., n:

M (q) i ℓ,ℓ ′ = t (q) iℓ t (q) iℓ ′ n (q) i ρ (q) i (x ℓ , x ℓ ′ )
where ρ (q)

i (x ℓ , x ℓ ′ ) can be computed through the kernel K as follows:

ρ (q) i (x ℓ , x ℓ ′ ) =K(x ℓ , x ℓ ′ ) - 1 n (q) i n j=1 t (q) ji (K(x j , x ℓ ) + K(x ℓ ′ , x j )) + 1 (n (q) i ) 2 n j,j ′ =1 t (q) ji t (q) j ′ i K(x j , x j ′ ).
The associated clustering algorithm, denoted by pgpEM in the following, provides at convergence an estimation of the set of parameters maximizing the pseudo-likelihood associated with the parsimonious Gaussian model used. Model selection can be done using likelihood penalized criteria which add to the log-likelihood a penalty depending on the complexity of the model. Classical criteria for model selection include the AIC [START_REF] Akaike | A new look at the statistical model identification[END_REF], BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] and ICL [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]

Numerical experiments

In this section, numerical experiments and comparisons are conducted on simulated and real-world data sets to highlight the main features of pgpDA and pgpEM.

Simulation study: influence of d and n

A two-class non linear classification problem is first considered, see Figure 1, in order to study the influence of both the intrinsic dimension d and the size n of the data on pgpDA. The data have been simulated in R 2 according to:

X |Z=1 = -1 + τ + η, 2 - τ 2 2 + η , X |Z=2 = 1 + τ + η, -2 + τ 2 2 + η ,
where τ ∼ U [-4, 4] and η ∼ N (0, 0.25). The first class is depicted by circles on Figure 1 whereas crosses correspond to observations of the second class. For all the experiments, the radial basis function (RBF) kernel was used and the kernel hyper-parameter was set to σ = 0.5. For models M 0 , M 2 and M 5 , the scree test threshold was fixed to 0.05. We remind that the Cattell scree test selects, for each class, the dimension for which the subsequent eigenvalues differences are smaller than a threshold. The decision boundaries for each pgpDA models are depicted in Figure 1. From panels (a)-(d), it appears that all the pgpDA models perform a non linear classification of the samples. For this toy data set, except M 7 and M 8 , all the models perform similarly and the decision boundaries are very satisfying.

M0 M1 M2 M3 M4 M5 (a) (b) M6 M7 M8 HDDA QDA LDA (c) (d)
As anyone could expect, conventional model-based classification methods (LDA, QDA, HDDA) are not be able to separate correctly those non linear data.

Let us first focus on the influence of the intrinsic dimensions on the pgpDA classifier. For the sake of simplicity, model M 1 is considered in this experiment since it assumes a common intrinsic dimension, let us say d, for all classes. The influence of this intrinsic dimension d on the decision boundaries for model M 1 is illustrated in Figure 2. For this experiment, we vary the intrinsic dimension value between 1 and 80. On the one hand, the parameter d seems, at least for this toy data set, not to have a strong influence since the decision boundaries are similar whatever the value of d. This behavior can be explained by the fact that all dimensions of the feature space are taken into considerations even though only the most informative are modeled precisely. On the other hand, it can be seen that when d is set to a low value, the decision boundary is slightly smoother than when it is set to a high value. Hence, we recommend the use of low values of d in practical situations in order to prevent over-fitting.

Table 3 shows the computing time for learning the pgpDA classifier according to the size of the learning set. Computing times for kernel Fisher discriminant analysis (KFD), support vector machines (SVM, with the SVM-KM toolbox, [START_REF] Canu | SVM and kernel methods matlab toolbox[END_REF]) and Gaussian process classification (GPML toolbox, (Rasmussen and Williams, 2006a)) are also provided for comparison. Notice that all codes used in this experiment are written in Matlab in order to obtain a fair comparison. One can first observe that pgpDA has a similar computing time as SVM. One can also remark that pgpDA has a comparable computing time as KFDA (though the latter is slightly faster), which is unsurprising since both methods are somehow related (see Section 5.3). Let us notice that models M 7 and M 8 of pgpDA are significantly more time consuming than the other models because their parameters (variances inside and outside the subspace, subspace orientation and intrinsic dimensions) are common between groups. As a consequence, the estimation of parameters for those models can not be done for each class independently and requires to work on a unique n × n matrix instead of several n i × n i matrices.

Finally, the GPML implementation of Gaussian process classification turns out to be dramatically more time consuming than pgpDA, KFDA and SVM. This difference in computing times is certainly explained by the need to use the iterative expectation-propagation algorithm for inferring the GP model whereas the inference is direct for KFDA and pgpDA.

Benchmark study on quantitative data

We focus here on the comparison of pgpDA with state-of-the-art methods. To that end, two kernel generative classifiers are considered, kernel Fisher discriminant analysis (KFD) [START_REF] Mika | Fisher discriminant analysis with kernels[END_REF] and kernel LDA (KLDA) [START_REF] Dundar | Toward an optimal supervised classifier for the analysis of hyperspectral data[END_REF], and one kernel discriminative classifier, support vector machine (SVM) [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]. The RBF kernel is used once again in the experiments for all methods, including pgpDA. Since real-world problems are considered, all the hyperparameters of the classifiers have been tuned using 5-fold cross-validation.

Eleven data sets from the UCI Machine Learning Repository (http://archive.ics.uci.edu/ml/) have been selected: glass, ionosphere, iris, sonar, USPS, wine, vowel, letter, protein, satimage and vehicle.

We chose these data sets because they represent a wide range of situations in term of number of observations n, number of variables p and number of groups k. The USPS dataset has been modified to focus on discriminating the three most difficult classes to classify, namely the classes of the digits 3, 5 and 8. This dataset has been called USPS 358. The main feature of the data sets are described in Table 4. Each data set was randomly split into training and testing sets in the hold-out ratio hr given in Table 4.

The data were scaled between -1 and 1 on each variable. The search range for the cross-validation was 5.

Regarding the competitive methods, KFD and SVM provide often better results than KLDA. The model used in KLDA only fits "ionosphere", "iris" and "wine" data, for which classification accuracies are similar to or better than those obtain with KFD and SVM. For the parsimonious pgpDA models, the classification accuracies are globally good. The results of models M 7 and M 8 are not presented here since they perform very poorly and we do not recommend their use. Model M 1 provides the best result in terms of average correct classification rates. In particular, for the "USPS 358", "wine" and "satimage" data sets, it provides the best overall accuracy. Let us remark that pgpDA performs significantly better than SVM (for the RBF kernel) on high-dimensional data (USPS 358).

In conclusion of these experiments, by relying on parsimonious models rather than regularization, pgpDA provides good classification accuracies and it is robust to the situation where few samples are available in regards to the number of variables in the original space. In practice, models M 1 and M 4 should be recommended: intrinsic dimension is common between the classes and the variance inside the class intrinsic subspace is either free or common. Conversely, models M 7 and M 8 must be avoided since they appeared to be too constrained to handle real classification situations.

Clustering of functional data: the Canadian temperatures

We now focus on illustrating the possible range of application of the proposed methodologies to different types of data. We consider here the clustering of functional data with pgpEM for which the mapping function ϕ is explicit (see Section 5.2). The Canadian temperature data used in this study, presented in details in [START_REF] Ramsay | Functional data analysis[END_REF] weather stations across the country. The pgpEM algorithm was applied here with the model M 0 , which is the most general parsimonious Gaussian process model proposed in this work. The mapping function ϕ consists in the projection of the observed curves on a basis of 20 natural cubic splines. Once the pgpEM algorithm has converged, various informations are available and some of them are of particular interest. Group means, intrinsic dimensions of the group-specific subspaces and functional principal components of each group could in particular help the practitioner in understanding the clustering of the dataset at hand. Since pgpEM is a generative clustering technique, it is possible to choose the most appropriate model for the data at hand using model selection criteria, such as BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] or ICL [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. The number of groups has been chosen here with the BIC criterion which picks k = 4 (see Figure 3). Figure 4 presents the clustering of the temperature data set into 4 groups with pgpEM.

It is first interesting to have a look at the name of the weather stations gathered in the different groups formed by pgpEM. It appears that group 1 (black solid curves) is mostly made of continental stations, group 2 (red dashed curves) mostly gathers the stations of the North of Canada, group 3 (green dotted curves) mostly contains the stations of the Atlantic coast whereas the Pacific stations are mostly gathered in group 4 (blue dot-dashed curves). For instance, group 3 contains stations such as Halifax (Nova Scotia) and St Johns (Newfoundland) whereas group 4 has stations such as Vancouver and Victoria (both in British Columbia). Figure 5 provides a map of the weather stations where the colors indicate their group membership. This figure shows that the obtained clustering with pgpEM is very satisfying and rather coherent with the actual geographical positions of the stations (the clustering accuracy is 71% here compared with the geographical classification provided by [START_REF] Ramsay | Functional data analysis[END_REF]). We recall that the geographical positions of the stations have not been used by pgpEM to provide the partition into 4 groups. An important characteristic of the groups, but not necessarily easy to visualize, is the specific functional subspace of each group. A classical way to observe principal component functions is to plot the group mean function μi (t) as well as the functions μi (t) ± 2 λij qij (t) (see [START_REF] Ramsay | Functional data analysis[END_REF] for more details). Figure 6 shows such a plot for the 4 groups of weather stations formed by pgpEM. The center and right panels show the cluster mean as well as the back-projections of the perturbed mean along the two first eigenfunctions (+ and -curves). For instance, the + curve of the center panel of 1st row is the back-projection of the curve μ1 (t) + 2 λ11 q11 (t). It first appears on the first functional principal component of each group that there is more variance between the weather stations in winter than in summer. In particular, the first principal component of group 4 (blue curves, mostly Pacific stations) reveals a specific phenomenon which occurs at the beginning and the end of the winter. Indeed, we can observe a high variance in the temperatures of the Pacific coast stations at these periods of time which can be explained by the presence of mountain stations in this group.

St

The analysis of the second principal components reveals finer phenomena. For instance, the second principal component of group 1 (black curves, mostly continental stations) shows a slight shift between the + andalong the year which indicates a time-shift effect. This may mean that some cities of this group have their seasons shifted, e.g. late entry and exit in the winter. Similarly, the inversion of the + andon the second principal component of the Pacific and Atlantic groups (blue and green curves)

suggests that, for these groups, the coldest cities in winter are also the warmest cities in summer. On the second principal component of group 2 (red curves, mostly Arctic stations), the fact that the + andcurves are almost superimposed shows that the North stations have very similar temperature variations (different temperature means but same amplitude) along the year. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + The first 168 congressmen are republicans whereas the 267 last ones are democrats.

+ + + + + -------------------------------------------------------------------------------------------------------------------------------- 0
+ + + + + -------------------------------------------------------------------------------------------------------------------------------- (a
+ + + + + - ------------------------------------------------------------------------------------------------------------------------------- 0
+ + + + + -------------------------------------------------------------------------------------------------------------------------------- (b)
+ + + + + -------------------------------------------------------------------------------------------------------------------------------- 0
+ + + + + -------------------------------------------------------------------------------------------------------------------------------- (

Clustering of categorical data: the house-vote dataset

We focus now on categorical data which are also very frequent in scientific fields. We consider here the task of clustering (unsupervised classification) and therefore the pgpEM algorithm. To evaluate the ability of pgpEM to classify categorical data, we used the U.S. House Votes data set from the UCI repository and compare its performance with kernel k-means [START_REF] Girolami | Mercer kernel-based clustering in feature space[END_REF][START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. This data set is a record of the votes (yea, nay or unknown) for each of the U.S. House of Representatives congressmen on 16 key votes in 1984. These data were recorded during the third and fourth years of Ronald Reagan's Presidency. At this time, the republicans controlled the Senate, while the democrats controlled the House of Representatives. Figure 7 shows the database where yeas are in indicated in white, nays in gray and missing values in black. The first 168 congressmen are republicans whereas the 267 last ones are democrats. As we can see, the considered votes are very discriminative since republicans and democrats vote differently in almost all cases while most of the congressmen follow the majority vote in their group. We can however notice that a significant part (around 50 congressmen) of the democrats tend to vote differently from the other democrats.

To cluster this dataset, we first build a kernel from the categorical observations (16 qualitative variables with 3 possible values: yea, nay or unknown). We chose a kernel, proposed in [START_REF] Couto | Kernel k-means for categorical data[END_REF], based on the Hamming distance which measures the minimum number of substitutions required to change one observation into another one. Naturally, pgpEM and kernel k-means worked on the same kernel to have a fair comparison and with a number of group equals to 2. The pgpEM algorithm was used with the model M 0 and the Cattell's threshold set to 0.2 (cf. Section 4.3). Figure 8 presents the clustering result obtained with kernel k-means (left) and pgpEM (right). The clustering results are presented through a binary matrix where a black pixel indicates a common membership between two congressmen and a white pixel means different memberships for the two congressmen. The clustering accuracy between the obtained partition of the data and the democrat/republican partition was 88.97% for pgpEM and 87.59% for kernel k-means on this example. As one can observe, both algorithms globally succeed in recovering the partition of the House of Representatives. It is also interesting to notice that most of the congressmen which are not correctly classified are those who tend to vote differently from the majority vote in their group. Nevertheless, pgpEM correctly classifies 6 congressmen more than kernel k-means on this run. To confirm the apparent good behavior of pgpEM, the experiment was replicated. Figure 9 presents boxplots of the clustering accuracies for kernel k-means and pgpEM over 25 replications. It turns out that pgpEM is significantly better than kernel k-means to cluster this kind of data.

Classification of mixed data: the Thyroid dataset

In this final experiment, we consider the supervised classification of mixed data which is more and more a frequent case. Indeed, it is usual to collect for the same individuals both quantitative and categorical data. For instance, in Medicine, several quantitative features can be measured for a patient (blood test results, blood pressure, morphological characteristics, ...) and these data can be completed by answers of the patient on its general health conditions (pregnancy, surgery, tabacco, ...). The Thyroid dataset considered here is from the UCI repository and contains thyroid disease records supplied by the Garavan Institute, Sydney, Australia. The dataset contains 665 records on male patients for which the answers (true of false) on 14 questions have been collected as well as 6 blood test results (quantitative measures). Among the 665 patients of the study, 61 suffer from a thyroid disease.

To make pgpDA able to deal with such data, we built a combined kernel by mixing a kernel based on the Hamming distance [START_REF] Couto | Kernel k-means for categorical data[END_REF] for the categorical features and a RBF kernel for the quantitative data. We chose to combine both kernels simply as follows:

K(x j , x ℓ ) = αK 1 (x j , x ℓ ) + (1 -α)K 2 (x j , x ℓ ),
where K 1 and K 2 are the kernels computed respectively on the categorical and quantitative features.

Another solution would be to multiply both kernels. We refer to [START_REF] Gönen | Multiple kernel learning algorithms[END_REF] for further details on multiple kernel learning.

The model for pgpDA was the model M 0 with the Cattell's threshold set to 0.2 (cf. Section 4.3). On each replication of the experiment, the data set was randomly split into a learning set of 599 observations and a test set of 66 observations. We selected the optimal set of kernel parameters by cross-validation on the learning part of the data. The average chosen parameter for the Hamming kernel was ξ = 0.5 (K 1 = exp(-D 2 )/ξ 2 where D is the Hamming distance matrix), the RBF kernel parameter was σ = 15 and the retained value for the mixing parameter α was 0.3. Figure 10 shows the choice of α and the kernels associated to the chosen values are presented in Figure 11. The rows and columns of the matrices are sorted according to the class memberships (healthy or sick) and the sick patients are the last ones. We then compared the performance of pgpDA with the combined kernel to pgpDA with, on the one hand, a simple RBF kernel built only on the quantitative variables of the dataset and, on the other hand, a Hamming kernel built only on the categorical variables. Table 6 presents both the true positive (TP) and false positive (FP) rates obtained on 25 replications of the classification experiment for pgpDA on quantitative data, on categorical data and on the mixed data.

It turns out that quantitative data contains most of the important information to discriminate the patients with thyroid diseases and that categorical data, when considered alone, are not enough to build an efficient classifier. However, it appears that the use of the categorical features in combination with the quantitative data allows to slightly improve the prediction of thyroid diseases (increases the TP rate and decreases the FP rate).

Conclusion

This work has introduced a family of parsimonious Gaussian process models for the supervised and unsupervised classification of quantitative and non-quantitative data. The proposed parsimonious models are obtained by constraining the eigen-decomposition of the Gaussian processes modeling each class. They allow in particular to use non-linear mapping functions which project the observations into an infinite dimensional space and to build, from a finite sample, a model-based classifier in this space. It has been also demonstrated that the building of the classifier can be directly done from the observation space through a kernel, avoiding the explicit knowledge of the mapping function. It has been possible to classify data of various nature including categorical data, functional data, networks and even mixed data. The methodology is as well extended to the unsupervised classification case. Numerical experiments on benchmark data sets have shown that pgpDA performs similarly or better compared to the best kernel methods of the state of the art. The possibility to examine the model parameters and to visualize the data into the class-specific feature subspaces permits a finer interpretation of the results than with conventional discriminative kernel methods. Among the possible extensions of this work, it would be interesting to extend the methodology to the semi-supervised classification case or to the classification with uncertain labels [START_REF] Bouveyron | Robust supervised classification with mixture models : Learning from data with uncertain labels[END_REF].

Appendix: Proofs

Proof of Proposition 1 Recalling that d max = max(d 1 , ..., d k ), the classification function can be rewritten as:

D i (ϕ(x)) = r j=1 1 λ ij < ϕ(x) -µ i , q ij > 2 L2 + di j=1 log(λ ij ) + dmax j=di+1 log(λ) -2 log(π i ) + γ,
where γ = (rd max ) log(λ) is a constant term which does not depend on the index i of the class. In view of the assumptions, D i (ϕ(x)) can be also rewritten as:

D i (ϕ(x)) = di j=1 1 λ ij < ϕ(x) -µ i , q ij > 2 L2 + 1 λ r j=di+1 < ϕ(x) -µ i , q ij > 2 L2 + di j=1 log(λ ij ) + (d max -d i ) log(λ) -2 log(π i ) + γ.
Introducing the norm ||.|| L2 associated with the scalar product < ., . > L2 and in view of Proposition 1 of (Shorack and Wellner, 1986, p. 208), we finally obtain:

D i (ϕ(x)) = di j=1 1 λ ij - 1 λ < ϕ(x) -µ i , q ij > 2 L2 + 1 λ ||ϕ(x) -µ i || 2 L2 + di j=1 log(λ ij ) + (d max -d i ) log(λ) -2 log(π i ) + γ,
which is the desired result.

Proof of Proposition 2

The proof involves three steps. i) Computation of the projection < ϕ(x)-μi , qij > L2 : Since ( λij , qij ) is solution of the Fredholm-type equation, it follows that, for all t ∈ J, λij qij (t) = ˆJ Σi (s, t)q ij (s)ds = 1 n i x ℓ ∈Ci < ϕ(x ℓ )μi , qij > L2 (ϕ(x ℓ )(t)μi (t)).

This implies that qij lies in the linear subspace spanned by the (ϕ(x ℓ )-μi ), x ℓ ∈ C i . As a consequence, the rank of the operator Σi is finite and is at most r i = min(n i , r). It therefore exists β ijℓ ∈ R such that: qij = 1

n i λij x ℓ ∈Ci β ijℓ (ϕ(x ℓ ) -μi ) (9)
leading to:

< ϕ(x) -μi , qij > L2 = 1 n i λij x ℓ ∈Ci β ijℓ ρ i (x, x ℓ ), (10)
for all j = 1, . . . , r i . The estimated classification function has therefore the following form:

Di (ϕ(x)) = 1 n i di j=1 1 λij 1 λij - 1 λ x ℓ ∈Ci β ijℓ ρ i (x, x ℓ ) 2 + 1 λ ρ i (x, x) + di j=1 log( λij ) + (d max -d i ) log( λ) -2 log(π i ),
for all i = 1, . . . , k.

ii) Computation of the β ijℓ and λij : Replacing (9) in the Fredholm-type equation ( 8) it follows that

1 n i x ℓ ,x ℓ ′ ∈Ci β ijℓ ′ (ϕ(x ℓ ) -μi )ρ i (x ℓ , x ℓ ′ ) = λij x ℓ ′ ∈Ci β ijℓ ′ (ϕ(x ℓ ′ ) -μi ).
Finally, projecting this equation on ϕ(x m )μi for x m ∈ C i yields

1 n i x ℓ ,x ℓ ′ ∈Ci β ijℓ ′ ρ i (x ℓ , x m )ρ i (x ℓ , x ℓ ′ ) = λij x ℓ ′ ∈Ci β ijℓ ′ ρ i (x m , x ℓ ′ ).
Recalling that M i is the matrix n i × n i defined by (M i ) ℓ,ℓ ′ = ρ i (x ℓ , x ℓ ′ )/n i and introducing β ij the vector of R ni defined by (β ij ) ℓ = β ijℓ , the above equation can be rewritten as M 2 i β ij = λij M i β ij or, after simplification M i β ij = λij β ij . As a consequence, λij is the jth largest eigenvalue of M i and β ij is the associated eigenvector for all 1 ≤ j ≤ d i . Let us note that the constraint qij = 1 can be rewritten as β t ij β ij = 1. iii) Computation of λ: Remarking that trace( Σi ) = trace(M i ) +

n i λij 1 n i x ℓ ′ ∈Ci (x ℓ ′ -μi )(x ℓ ′ -μi ) t x ℓ ∈Ci β ijℓ (x ℓ -μi ) = 1 n i λij 1 n i x ℓ ′ ,x ℓ ∈Ci β ijℓ (x ℓ ′ -μi )(x ℓ ′ -μi ) t (x ℓ -μi ) = 1 n i λij 1 n i x ℓ ′ ,x ℓ ∈Ci β ijℓ (x ℓ ′ -μi )ρ i (x ℓ , x ℓ ′ ) = 1 n i λij x ℓ ′ ,x ℓ ∈Ci (M i ) ℓ,ℓ ′ β ijℓ (x ℓ ′ -μi ) = 1 n i λij B -1 x ℓ ′ ∈Ci (M i β ij ) ℓ ′ (x ℓ ′ -μi ),
and remarking that β ij is eigenvector of M i , it follows:

Σi qij = λij 1

n i λij B -1 x ℓ ′ ∈Ci β ijℓ ′ (x ℓ ′ -μi ) = λij qij .
Second, straightforward algebra shows that

||q ij || 2 = 1 n i λij x ℓ ∈Ci β ijℓ (x ℓ ) -μi ) t x ℓ ′ ∈Ci β ijℓ ′ (x ℓ ′ -μi ) = 1 n i λij x ℓ ′ ,x ℓ ∈Ci β ijℓ β ijℓ ′ (x ℓ -μi ) t (x ℓ ′ -μi ) = 1 λij x ℓ ′ ,x ℓ ∈Ci (M i ) ℓ,ℓ ′ β ijℓ β ijℓ ′ = 1 λij qt ij M i qij = 1,
and the result is proved. As a consequence, ϕ(x) = B -1 γ(x) and K(x, y) = γ(x) t B -1 γ(y). Introducing γi = 1 n i xj ∈Ci γ(x j ), it follows that ρ i (x, y) = (γ(x)γi ) t B -1 (γ(y)γi ). Let us first show that qij is eigenvector of

Proof of

B -1 Σi . Recalling that qij = 1 n i λij B -1 x ℓ ∈Ci β ijℓ (γ(x ℓ ) -γi ),
we have

B -1 Σi qij = 1 n i λij B -1 1 n i x ℓ ′ ∈Ci (γ(x ℓ ′ ) -γi )(γ(x ℓ ′ ) -γi ) t B -1 x ℓ ∈Ci β ijℓ (γ(x ℓ ) -γi ) = 1 n i λij B -1 1 n i x ℓ ′ ,x ℓ ∈Ci β ijℓ (γ(x ℓ ′ ) -γi )(γ(x ℓ ′ ) -γi ) t B -1 (γ(x ℓ ) -γi ) = 1 n i λij B -1 1 n i x ℓ ′ ,x ℓ ∈Ci β ijℓ (γ(x ℓ ′ ) -γi )ρ i (x ℓ , x ℓ ′ ) = 1 n i λij B -1 x ℓ ′ ,x ℓ ∈Ci (M i ) ℓ,ℓ ′ β ijℓ (γ(x ℓ ′ ) -γi ) = 1 n i λij B -1 x ℓ ′ ∈Ci (M i β ij ) ℓ ′ (γ(x ℓ ′ ) -γi ).
Remarking that β ij is eigenvector of M i , it follows:

B -1 Σi qij = λij 1 n i λij B -1 x ℓ ′ ∈Ci β ijℓ ′ (γ(x ℓ ′ ) -γi ) = λij qij .
Let us finally compute the norm of qij :

||q ij || 2 = 1 n i λij x ℓ ∈Ci β ijℓ (γ(x ℓ ) -γi ) t B -1 x ℓ ′ ∈Ci β ijℓ ′ (γ(x ℓ ′ ) -γi ) = 1 n i λij x ℓ ′ ,x ℓ ∈Ci β ijℓ β ijℓ ′ (γ(x ℓ ) -γi ) t B -1 (γ(x ℓ ′ ) -γi ) = 1 λij x ℓ ′ ,x ℓ ∈Ci (M i ) ℓ,ℓ ′ β ijℓ β ijℓ ′ = 1 λij qt ij M i qij = 1,
and the result is proved.

Proof of Proposition 5

Let us first set the estimate of µ i at iteration s to its empirical counterpart conditionally on the current value of the model parameter θ (s-1) :

μ(s) i (t

  criteria. The Bayesian Information Criterion (BIC) is certainly the most popular and consists in selecting the model which penalizes the likelihood by γ(M) 2 log(n) where γ(M) is the number of parameters in model M and n is the number of observations. On the other hand, the AIC criterion penalizes the log-likelihood by γ(M) whereas the ICL criterion add the penalty n i=1 K k=1 t ik log(t ik ) to the one of the BIC criterion in order to favor well separated models. The value of γ(M) is of course specific to the model selected by the practitioner.

Figure 1 :

 1 Figure 1: Decision boundaries for the simulated classification problem are depicted in color: (a) models M 0 , M 1 and M 2 , (b) models M 3 , M 4 and M 5 , (c) models M 6 , M 7 and M 8 , (d) QDA, LDA and HDDA.

Figure 2 :

 2 Figure 2: Simulated classification problem: Influence of the intrinsic dimension for the model M 1 . In blue d = 1, in red d = 5, in black d = 10, in magenta d = 20, in yellow d = 40 and in green d = 80.

Figure 3 :

 3 Figure 3: Choice of the number k of groups with the BIC criterion for the Canadian temperature data set.

  + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -----------------------------------------------------------------------------------------------------------------------------++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + --------------------------------------------------------------------------------------------------------------------------------(c) Group 3 (mostly Atlantic stations)

Figure 6 :

 6 Figure 6: The group means of the Canadian temperature data obtained with pgpEM (left panel) and the effect of adding (+) or subtracting (-) twice the square root of the variance to the cluster means along the two first eigenfunctions (center and right panels, see text for details).

Figure 7 :

 7 Figure 7: Votes (yea, nay or unknown) for each of the U.S. House of Representatives congressmen on 16 key votes in 1984. Yeas are in indicated in white, nays in gray and missing values in black. The first 168 congressmen are republicans whereas the 267 last ones are democrats.

Figure 8 :

 8 Figure 8: Clustering results of the house-vote dataset obtained with kernel k-means (left) and pgpEM (right). The clustering results are presented through a binary matrix where a black pixel indicates a common membership between two congressmen. Again, the first 168 congressmen are republicans whereas the 267 last ones are democrats. On this run, pgpEM correctly classifies 6 congressmen more than kernel k-means.

Figure 9 :Figure 10 :

 910 Figure 9: Boxplots of the clustering accuracies for kernel k-means (left) and pgpEM (right) over 25 replications on the house-vote dataset.

  It is sufficient to prove that qij and λij are respectively the jth normed eigenvector and eigenvalue of Σi . First, Σi qij = 1

B

  Proposition 4 For all ℓ = 1, . . . , L, the ℓth coordinate of the mapping function ϕ(x) is defined as the ℓth coordinate of the function x expressed in the truncated basis {b 1 , . . . , b L }. x)b ℓ (t), for all t ∈ [0, 1] and thus, for all j = 1, . . . , L, we have γ j (x) = ˆ1 0 x(t)b j (t)dt = jℓ ϕ ℓ (x).

  It indicates the memberships of the learning data to the k classes denoted by C 1 , . . . , C k , i.e., z ℓ = i indicates that x ℓ belongs to C i .

The class labels {z 1 , ..., z n } are assumed to be realizations of a discrete random variable Z ∈ {1, ..., k}.

Table 1 :

 1 Submodels of the parsimonious Gaussian process model (referred to as M 0 here). Letting d max = max(d 1 , ..., d k ), the classification function D i can be written as follows in the case of a parsimonious Gaussian process model pgpM 0 :

	Model	Variance inside the subspace Fi	Variance outside Fi	Subspace orientation Qi	Intrinsic dimension di
	M0	Free	Common	Free	Free
	M1	Free	Common	Free	Common
	M2	Common within groups	Common	Free	Free
	M3	Common within groups	Common	Free	Common
	M4	Common between groups	Common	Free	Common
	M5	Common within and between groups	Common	Free	Free
	M6	Common within and between groups	Common	Free	Common
	M7	Common between groups	Common	Common	Common
	M8	Common within and between groups	Common	Common	Common
	pgpM 0 (or M 0 for short) hereafter. With these assumptions, we have the following result (proof is provided in the appendix).
	Proposition 1.

Table 3 :

 3 Average computing time (in seconds, ± standard deviation) over 30 replications for learning the classifiers KFDA, SVM, GPML and pgpDA according to the size of the learning set.

	Dataset	n	p	n/p	k	hr
	Iris	150	4	37.5	3	0.5
	Glass	214	9	23.7	6 0.25
	Wine	178	13	13.7	3	0.5
	Ionosphere	351	34	10.3	2	0.5
	Sonar	208	60	3.5	2	0.5
	USPS 358	2248 256	8.8	3	0.5
	Vowel	990	10	99	10 0.5
	Letter	15000 16 937.5 26 0.1
	Protein	17766 357 49.8	3 0.05
	Satimage	4435	36 123.2 6 0.25
	Vehicle	846	18	47	4	0.5

Table 4 :

 4 Data used in the experiments. n is the number of samples, p is the number of features, k is the number of classes and hr is the hold-out ratio used in the experiments.

  , consist in the daily measured temperatures at 35 Canadian

	Method	Iris	Glass	Wine	Ionosphere	Sonar	USPS 358	Vowel	Letter	Protein	Satimage	Vehicle	Mean (rank)
	pgpDA M 0	95.9± 2.1	64.9 ± 6.3	96.8 ± 1.7	90.5 ± 2.3	77.9 ± .9	92.2 ± 1.0	95.2 ± 1.6	86.5 ± 0.6	54.8 ± 1.3	89.4 ± 0.5	78.3 ± 2.7	83.85 (4)
	pgpDA M 1	95.2± 2.1	62.6 ± 12.5	96.7 ± 2.3	93.7 ± 1.6	81.8 ± 4.9	96.6 ± 0.4	95.2 ± 1.6	85.4 ± 0.8	57.6 ± 0.9	89.4 ± 0.5	79.9 ± 1.9	84.92 (1)
	pgpDA M 2	94.4± 6.2	64.4 ± 6.7	96.8 ± 1.8	91.0 ± 2.8	71.6 ± 13.4	95.4 ± 0.8	91.2 ± 1.6	79.0 ± 0.7	54.6 ± 1.0	66.3 ± 2.5	66.8 ± 2.7	79.23 (8)
	pgpDA M 3	95.8± 2.3	64.3 ± 6.8	96.9 ± 2.0	93.2 ± 2.1	79.3 ± 4.9	96.2 ± 0.5	92.2 ± 1.5	78.7 ± 1.1	57.2 ± 1.0	77.6 ± 2.0	68.4 ± 2.3	81.80 (5)
	pgpDA M 4	94.4± 2.2	65.3 ± 6.4	97.2 ± 1.8	93.4 ± 2.0	81.6 ± 4.5	96.3 ± 0.7	95.0 ± 1.5	84.8 ± 0.7	58.2 ± 1.3	87.8 ± 0.9	72.8 ± 1.9	84.25 (3)
	pgpDA M 5	94.2± 7.1	59.8 ± 10.9	96.4 ± 2.0	92.0 ± 1.8	72.5 ± 12.6	96.0 ± 0.5	92.0 ± 1.5	79.0 ± 0.9	48.5 ± 1.9	63.4 ± 3.1	66.0 ± 2.9	78.16 (10)
	pgpDA M 6	94.8± 2.1	65.2 ± 5.6	97.2 ± 1.8	92.5 ± 2.1	79.8 ± 4.9	96.1 ± 0.5	92.4 ± 1.8	79.2 ± 1.0	56.3 ± 1.3	75.6 ± 2.3	66.9 ± 2.5	81.45 (6)
	KFD	93.4± 3.7	47.3 ± 10.1	95.9 ± 2.3	94.1 ± 1.7	82.9 ± 3.1	93.6 ± 0.5	95.1 ± 1.3	83.6 ± 0.9	53.5 ± 3.8	86.1 ± 0.7	66.4 ± 12	81.08 (7)
	KLDA	96.6± 2.3	64.5 ± 6.3	96.6 ± 1.7	88.1 ± 2.3	68.9±18.1	64.7 ± 37.5	92.4 ± 1.4	85.2 ± 0.9	41.5 ± 10.4	86.1 ± 0.7	82.1 ± 1.4	78.79 (9)
	SVM	95.7± 2.0	69.1 ± 5.5	96.8 ± 1.4	92.8 ± 1.8	84.8 ± 4.0	77.6 ± 5.4	96.2 ± 1.5	87.0 ± 0.5	60.7 ± 0.9	88.5 ± 0.5	82.7 ± 1.5	84.72 (2)

Table 5 :

 5 Classification results on real-world datasets: reported values are average correct classification rates and standard deviation computed on validation sets. Clustering of the 35 times series of the Canadian temperature data set into 4 groups with pgpEM. The colors indicate the group memberships: group 1 in black, group 2 in red, group 3 in green and group 4 in blue.

	40
	20
	0
	-20
	-40
	0	100	200	300
	Figure 4:

  Geographical positions of the weather stations according to their group belonging. The colors indicate the group memberships: group 1 in black, group 2 in red, group 3 in green and group 4 in blue.

			Resolute
	Inuvik
	Dawson	Yellowknife	Iqaluit
	Whitehorse		Churchill Uranium Cty
	Pr. George Pr. Rupert	The Pas Pr. Albert Edmonton	Scheffervll
	Yarmouth Charlottvl London Toronto Halifax Sydney Fredericton Sherbrooke Montreal Ottawa Quebec . Johns Arvida Bagottville Thunderbay Winnipeg Regina Calgary Kamloops Vancouver Victoria
	Figure 5:

Table 6 :

 6 Classification results on test sets for the Thyroid dataset (mixed data). Results are averaged on 25 replications of the experiment.

	Method	pgpDA on quantitative data	pgpDA on categorical data	pgpDA with the combined kernel
	TP rate	74.86	96.00	75.88
	FP rate	22.16	95.53	21.97

Acknowledgments

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

According to Bayes' rule and Equation ( 1), the computation of t (s) ij = E(Z j = i|x j , θ (s-1) ) conditionally on the current value of the model parameter θ (s-1) can be finally written as:

, for j = 1, . . . , n and i = 1, . . . , k. The result follows.