Multiclass Probabilistic Kernel Discriminant Analysis

1 Introduction

Learning in high-dimensional spaces is challenging due to the curse of dimensionality [Hastie et al., 2001]. Linear and kernel discriminant analysis provide effective means to reduce dimensionality [McLachlan, 1992]. As an extension of the linear discriminant analysis (LDA), kernel discriminant analysis (KDA) extends LDA to the kernel-induced feature space [Scholk¨opf and Smola, 2002], so that the non-linear structures in the data can be handled effectively [Mika et al., 1999]. KDA has been shown to be effective in many applications, such as image processing [Belhumeur et al., 1997] and text information retrieval [Howland et al., 2003], where KDA is used to generate low dimensional representations of the original data for subsequent analysis. KDA is known to be prone to overfitting, which causes the model to be sensitive to noise. One effective way to address the problem is to apply regularization for variance reduction [Friedman, 1989]. Various model selection approaches have been proposed to determine a good regularization parameter value from a finite set of candidate values based on cross-validation, which is, however, time consuming. In addition, the performance of cross-validation could also be influenced by the quality of the candidate set. Therefore it is important to study how to automatically determine the model parameters for KDA to improve its robustness.

Probabilistic models have been proposed for LDA and KDA. For example in [Centeno and Lawrence, 2006], by relating Rayleighs coefficient to a noise model, the authors showed that their model is equivalent to KDA. The authors also demonstrated that through a Bayesian model selection approach [Gelman et al., 1995], parameters in their model can be efficiently tuned to avoid the cost of cross-validation. However, due to the coding scheme used in designing, existing models are limited to binary class problems, which significantly confines their application in real world problems. To address the limitation, in this paper, we develop a nonparametric probabilistic model for KDA based on Gaussian Process (GP) to handle multiclass problems. We show that under a mild assumption, which holds in most real applications, the proposed model is equivalent to KDA, when the noise terms approach zero. The proposed model provides a novel Bayesian interpretation for KDA, which allows its parameters to be automatically tuned through the optimization of the marginal log-likelihood of the data. We show that the optimization problem can be formulated as a DC-programming problem [Horst and Thoai, 1999], which can be solved efficiently. As a probabilistic model based on GP, many well studied techniques for GP can be directly applied to the proposed model and enable it, for example, to handle problems with large scale data [Snelson, 2007], and learn from multiple heterogenous data sources [Lanckriet et al., 2004; Rasmussen and Ghahramani, 2002]. The proposed model substantially extends the capability of KDA. To evaluate its performance, we conduct experiments on benchmark data sets. Experimental results demonstrate its efficacy.

2 Background

We first define the notations used in the rest of this paper. We use φ(X) = {φ(x1) , . . . , φ(xn)} to denote a data set of n samples, where φ(xi) is the mapping of the ith sample xi in the kernel-induced feature space. Assume the data has c classes, and y = {y1, . . . , yn} denotes the class label, with yi ∈ {1, . . . , c} being the label of the ith sample. Let ki,j = φ(xi)T φ(xj), and let K denote the kernel matrix with ki,j as its i-j th element. Let P = I - n-111T be the centering matrix. Kc = PKP denotes the centered kernel matrix, such that samples in the feature space induced from Kc are centered: φc(x) = φ(xi)- ¯ φ(x), ¯ φ(x) = 1n ∑ φ(xi).

In this paper, we use uppercase characters to denote matrices, boldface lowercase characters to denote vectors, and standard lowercase characters to denote scalars. Also, we use I to denote the identity matrix and 1 the vector of all ones. Below we give a brief introduction to KDA and GP.

Kernel Discriminant Analysis (KDA): Given the kernel matrix K and the class label y, KDA determines a transformation matrix B to project samples, such that they can be best separated. Analogous to LDA, given Kc, the centered K, we denote SK t = KcKc as the total scatter matrix, SK b = n-1KcY Y T Kc as the between-class scatter matrix, and SK w = SK t - SK b as the within-class scatter matrix. In the definition of SK b , Y ∈ Rn×c is the coding matrix derived from y as [Ye, 2007]:

KDA maximizes the separability of the samples in the dimensionality-reduced space by simultaneously minimizing trace(BT SK w B) and maximizing trace(BT SK b B), which correspond to the within-class distance and the between-class distance, respectively. KDA solves the following optimization problem:

However, the above formulation is prone to overfitting. Regularization is commonly applied to alleviate the problem:

And its optimal solution can be obtained by computing the principal eigenvectors of the following eigenvalue problem:

Here (·)+ denotes the matrix pseudo inverse [Golub and Van Loan, 1996]. Choosing a proper λ value is crucial to KDA.

Gaussian Process (GP) Regression: A Gaussian Process (GP) is a stochastic process defining a nonparametric prior over functions [Rasmussen and Williams, 2006]. A real function f : Rd → R follows a GP, denoted as GP(h; κ), if for any finite number of data points x1, . . . , xn, f = {f (xi)}ni=1 follows a multivariate Gaussian distribution N(h,K), with a mean function h = {h(xi)}ni=1 and a covariance matrix K = {κ(xi, xj)}ni,j=1 . It is common to set h ≡ 0. Given the training data and a test point x∗, GP regression assumes p(y∗|x∗) = ∫ p (y∗|f∗) p (f∗) df∗, where f∗ = f (x∗) and p (y∗|f∗) is an isotropic Gaussian with the variance λ specifying the system noise. By marginalizing f∗, we can obtain the joint distribution of y and y∗, which takes the form:

where σ corresponds to the system noise. Using Equation (5), we can obtain the conditional distribution of p (f∗ | y ), which is also a Gaussian with its mean and variance given by:

3 PKDA: A Probabilistic Model for Kernel Discriminant Analysis

We propose a probabilistic model for KDA based on Gaussian Process. We call the model, "Probabilistic Kernel Discriminant Analysis" or PKDA. Existing probabilistic models for discriminant analysis, such as the ones proposed in [Ioffe, 2006] and [Centeno and Lawrence, 2006] rely on regressing the samples to [-1, +1] or [n/n1,-n/n2], where n1 and n2 are the numbers of samples in positive and negative classes, respectively. The coding scheme only works for binary class problems. In PKDA, we propose to regress samples to Y defined in Equation (1). We can show that assuming the samples are linearly independent in the kernel-induced feature space, PKDA is equivalent to KDA, when the noise terms approach zero. We first present the PKDA model as follows:

Definition 1 Given a centered kernel matrix Kc, the Probabilistic Kernel Discriminant Analysis, or the PKDA, consists of a set of c Gaussian Processes,

where c is the number of classes, and the lth GP, GPl, uses yl ∈ Rn as its observations with yl being the lth column of Y defined in Equation (1).

Figure 1 shows the graphical model of PKDA. The model contains c GPs making predictions with noisy observations. In the model, the c GPs share the same covariance matrix K. For the lth GP, λl = σ2 l is its noise term, which is analogous to the regularization parameter in KDA, and hl is its mean vector. Following the convention, we fix hi ≡ 0. yl is its observation, which is the lth column of Y , where Y = [y1, . . . , yc] is defined as in Equation (1). Given a test data point φ (x∗), the PKDA model projects φ (x∗) to a vector of c dimensions: y∗ = (y∗,1, . . . , y∗,c) with the ith element following the predictive distribution specified as:

and mi (x∗) and σ2 i (x∗) are given by Equations (7) and (8). With the computed mi (x∗) and σ2 i (x∗), the projections of the test point can be either sampled from the normal distributions specified in Equation (10), or obtained directly using m(x∗) = (m1(x∗), . . . ,mc(x∗)), which maximizes the likelihood of the observations. In theorem 1 below, we show the equivalent relationship between KDA and PKDA under certain mild conditions.

Theorem 1 Assume that samples are linearly independent in the kernel-induced feature space. When the noise terms λl → 0, l = 1, . . . , c, the projection determined by the expectation of the predictive distributions of PKDA is equivalent to that generated by KDA. More specifically, for m(x∗) = (m1 (x∗) , . . . ,mc (x∗))T and k∗ = KX,x∗ , we have m(x∗) = BT k∗ up to an expansion with a dummy variable and an orthogonal transformation.

Here the equivalence means that the distance among samples under different projections are the same 1 . Various projections can be obtained by applying orthogonal transformations on an existing projection or by increasing the dimensionality of the inputs by adding dummy variables which has 0 as their only value. We first present two lemmas, which pave the way for the proof. The first lemma tells that the assumption used in the theorem is indeed very mild.

Lemma 1 When the RBF kernel function is used, as long as x1, . . . , xn are all distinct, the kernel matrix K is of full rank.

Proof of the lemma can be found in [Micchelli, 1984]. The Lemma tells that when RBF kernel function is used, as long as the samples are different, K will be of full rank, which means that the samples will be linearly independent in the feature space induced by K. In real applications, it is usually sensible to assume that the given samples are all distinct, therefore when RBF kernel function is used, the linear independent assumption will always hold. Let the compact SVD [Golub and Van Loan, 1996] of Kc be Kc = U1ΣtU T 1 , and the full SVD of U T 1 Y be U T 1 Y = PΣbQ, we have:

Lemma 2 When samples are linearly independent in the kernel induced feature space, we have Σ2b = diag(1, . . . , 1, 0).

Proof: Let SK w = SK t - SK b = Kc(I - 1n Y Y T )Kc. I - n-1Y Y T is positive semidefinite (psd) since:

This means SK w is also psd. Let G be a matrix defined as:

where U consists of the whole set of singular vectors of Kc, and U2 contains the singular vectors corresponding to the zero singular values. Under the assumption, we know that rank(K)=n and rank(Kc)=n1. Therefore U2 contains only one column: U2 = [(1/√n)1]. With the definition of G, it can be shown that, GT SK t G = diag (I, 0) , GT SK b G = diag (Σ2b , 0) and GT SK w G = diag(Σw, 0), therefore we have: diag (I, 0) = diag(Σw, 0) + diag(Σb, 0). Since SK w and SK b are all positive semidefinite, the diagonal elements of Σw and Σb must be nonnegative. It is easy to verify that rank(SK b )≤ c-1 and rank(Sw) ≤ n - c. Therefore we have: rank(diag(I, 0)) ≤ rank(diag(Σw, 0)) + rank(diag(Σb, 0))=n - 1. Thus, we can conclude that all nonzero diagonal entries in Σb are 1. Note that rank(U T 1 Y ) = rank(Y ) = c - 1 since U T 1 contains orthonormal columns, thus the number of nonzero diagonal entries in Σb is c1. This completes the proof. ■

1 This definition for equivalence is sensible , since in many classification and clustering approaches, such as the SVM and k-means, only the distance among samples are used to fit model.

We are now ready to prove the theorem, which establishes the equivalence between KDA and PKDA.

Proof of Theorem 1: We can show that for any input φ(x∗), ml(x∗) = kT∗ (Kc +λlI)-1 yl. Therefore, when λl → 0, the expectation of the predictive distributions in PKDA actually projects data with a transformation matrix defined as ˆB = (Kc)+ Y . The equivalence between KDA and PKDA can be established by studying the relationship between ˆB and B. We first study the structure of matrix ˆB. ˆB = K+ c Y = U1Σ-1 t U T 1 Y . Recall that U T 1 Y = PΣbQ. Let P = [p1, . . . , pc]. We have:

Next we see the structure of B can be expressed as:

Since P T ΣtU T 1 U1Σ-1 t P = I, it can be verified that the top c - 1 eigenvectors of (SK t )+ SK b are given by the first c - 1 columns of U1Σ-1 t P. Therefore we have:

As Q is orthogonal, and Σ2b =diag(1, . . . , 1, 0), the two projections are essentially equivalent. ■

3.1 Efficient Model Selection for PKDA

Theorem 1 establishes the connection between PKDA and KDA. Given a kernel K and the class label vector y, PKDA projects test points onto the reduced space where the data can be well separated. Below we show how to automatically determine the proper value for the regularization parameters (or the noise terms) by minimizing the negative log marginal likelihood. Denote ˜K = Kc + λI, where λ ≥ 0. For each GP in PKDAmodel, its negative log marginal likelihood is:

-

Minimizing the negative log marginal likelihood specifies the following nonlinear constrained optimization problem:

Let Kc = UΣU T , U = [u1, . . . , uN], yi = uTi y and Σ =diag(α1, . . . , αN). We have:

-

In the equation, the constant term n2 log 2π is ignored. The equation shows that the optimization problem specified in Equation (15) is a DC (difference of convex functions) Programming [Horst and Thoai, 1999] problem, and efficient techniques such as the DC algorithm (DCA) [An and Tao, 2005] can be applied. To minimize Equation (16), we define:

where

Let g∗ (z) = supλ {zλ - g (λ)} be the conjugate function of g(λ) and ∂ be the subdifferential operator, it is easy to see that the following equations for ∂h(λ) and ∂g∗ (z) hold:

Equation (21) can be solved by minimizing the convex function g (λ) - z1λ, subject to the nonnegative constrains on λ. The above equations, can be used in DCA proposed in [An and Tao, 2005] to compute solutions for the optimization problem specified in Equation (16). In our experiments we found that DCA usually returns globally optimal solutions.

The model parameters of PKDA can also be tuned by optimizing the leave-one-out log predictive probability [Rasmussen and Williams, 2006] in a similar way.

Time Complexity of PKDA: It turns out that the model selection step for PKDA does not increase the computation cost significantly, since the SVD of Kc generated in the model selection step is also used in the subsequent prediction step. We provide the complexity analysis for PKDA. Assume there are n training and m test samples. Calculating the SVD for Kc has a cost of O(n3 ). The results from SVD will be used for model selection and the prediction of the test points. DC-Programming can usually be solved in O(n3 ) operations. Given the SVD of Kc, computing its inverse costs O(n3 ) operations and the projection step for the m test samples costs O(cmn2 ) operations. Therefore, the total time complexity of PKDA is O(max(n, cm)n2 ). The complexities of PKDA and KDA are of the same order.

4 Empirical Study

In this section, we empirically evaluate the performance of PKDA. Nine data sets are used in the experiments. They are: two image data sets, AR10P 2 and ORL10P 3 ; five UCI data sets [Murphy and Aha, 1994]: DNA, SPLICE, SOLAR, SOYBEAN and ADVERTISE; and two text data sets from the 20-news-group data: TXT4C and TXT2C 4 . The details of the benchmark data sets are summarized in Table 1. The last column of Table 1 lists the rank of the kernels obtained on each data set using the RBF kernel. It shows that all kernels are of full rank, that is, samples are linearly independent in the kernel-induced space. The results suggest that the condition assumed in Theorem 1 is indeed very mild.

In the experiment, the parameters in PKDA are obtained in two ways: using the proposed automatic tuning mechanism (PKDA-OP) and using cross-validation (PKDA-CV). We compare PKDA with KDA, whose parameters are also obtained in two ways: using a default value (λ=0, KDA) and using cross-validation (KDA-CV). To evaluate the qualities of the dimensionality-reduced spaces found by different approaches, we train a 1-nearest-neighbor (1-nn) classifier in each space, and use its accuracy to determine how well the samples from different classes are separated in that space. A higher accuracy suggests a higher quality of the space in terms of sample separability. We also train a support vector machine (SVM) [Vapnik, 1995] using the given kernel and record its accuracy. Note that SVM is a classifier rather than a dimension reduction approach. We include SVM in the experiment to verify whether the accuracy of the 1-nn obtained in the dimensionality-reduced space is reasonably good. The parameter of the SVM is obtained in two ways: using a default parameter (C=1, SVM) and using cross-validation (SVM-CV). We implemented PKDA and KDA in Matlab and LibSVM [Chang and Lin, 2001] is used for SVM. We use the RBF kernel function to construct the kernels. Each algorithm is tested for 25 times on each data set by randomly sampling (at most) 200 instances from each class and we split the data into training and test sets of a ratio 4:1. The obtained averaged accuracy rates are presented in the paper. For PKDACV, KDA-CV and SVM-CV we apply 5-fold cross-validation on the training data. For CV based model selection, choosing the candidate parameter values is a difficult problem. In the experiments, to improve the probability that there is at least one good regularization parameter value for each CV based algorithm on each data set, we picked 20 candidates from a wide range between 0 and 104 .

2 http://rvl1.ecn.purdue.edu/ aleix/aleix face DB.html. Data set is subsampled down to the size of 60×40 = 2400

3 http://www.uk.research.att.com/facedatabase.html. Data set is subsampled down to the size of 100×100 = 10000

4 TXT4C: Baseball, Hockey, PC and MAC; TXT2C: Religion and Atheism. http://people.csail.mit.edu/jrennie/20Newsgroups/

4.1 Experimental Results

Below we present the experimental results comparing different approaches on the quality of the dimensionality-reduced spaces they generated, and their efficiency.

Accuracy Comparison

Table 2 presents the accuracy of different approaches on the nine benchmark data sets. Based on the accuracy results, we summarize the following observations.

First, compared to KDA using the default parameter value λ = 0, PKDA-OP performs much better. On 7 out of 9 data sets, PKDA-OP performs significantly better than KDA. We observed that on several data sets, such as the DNA and the SOYBEAN data, the performance of KDA is poor. This is likely due to overfitting. This suggests that it is necessary to apply regularization to make KDA more robust.

Second, when comparing PKDA-OP with PKDA-CV, we observed that the two approaches perform equally well. Compared to cross-validation, the automatic tuning mechanism does not need to run PKDA multiple times on different splits of the training data, therefore it is more efficient. In our experiments, we observed that the computational cost of PKDAOP is usually comparable to that of KDA and is significantly lower than that of KDA-CV and PKDA-CV. It is known that picking the right candidate parameter values in crossvalidation may not be easy. The proposed tuning process can also automatically find good values for the parameters.

Third, in comparison with KDA-CV, we found that PKDA is significantly better on 4 out of 9 data sets. We also noticed that comparing KDA with KDA-CV, KDA-CV performs significantly better, which verifies the effectiveness of crossvalidation as well as the regularization.

Finally, comparing PKDA-OP with SVM-CV, we observed that by using the dimensionality-reduced spaces generated by PKDA-OP, the 1-nn classifier is able to achieve accuracy rates comparable to those of the SVM with cross-validation. This result clearly shows the high quality of the dimensionalityreduced spaces generated by PKDA-OP.

Overall, on the 9 benchmark data sets, PKDA-OP+1nn achieved the best average accuracy of 0.88 which is followed by SVM-CV (0.87), PKDA-CV+1nn (0.87), SVM (0.84), KDA-CV+1nn (0.83) and KDA+1nn (0.67).

Efficiency Comparison

On the nine benchmark data sets, the averaged running time of PKDA-OP, KDA PKDA-CV and KDA-CV are: 1.58s, 1.27s, 114.63s and 114.04s, respectively. PKDA-OP has almost the same running time as KDA, while having the advantage of being able to automatically tune its model parameters. Figure 2 plots the running time of PKDA-OP, KDA, PKDACV, KDA-CV on five benchmark data sets, on which algorithms have relatively longer running time than on the other four benchmark data sets. The plots in the figure show that compared with cross-validation, the automatic tuning mechanism is significantly faster. We observed similar trends on the other four data sets. The results on efficiency comparison is consistent with our complexity analysis for PKDA.

The results demonstrated that the proposed multiclass probabilistic model for KDA is robust and efficient.

5 Conclusion

In this paper, we proposed a probabilistic model for KDA, which is able to handle multiclass problems. The proposed model is based on Gaussian Process and its model parameters can be automatically tuned in an efficient way. Experimental results demonstrated that the proposed model offers good performance, and is very efficient. Based on Gaussian Process, the proposed model can handle large scale data via approximation methods, such as BCM [Tresp, 2000] and SPGP [Snelson and Ghahramani, 2006]. Also as a probabilistic model, PKDA allows the integration of multiple heterogeneous data sources via Bayesian mixture models [Svensen and Bishop, 2005], which leads to an interesting nonlinear kernel combination formulation: Mixture of Discriminant Gaussian Process (MPKDA). Our preliminary results show that MPKDA outperforms the existing linear kernel combination approaches [Ye et al., 2008]. This forms one line of our ongoing research work.

Acknowledgment

This work was supported by NSF IIS-0612069, IIS-0812551, CCF-0811790, NIH R01-HG002516, and NGA HM1582-081-0016.