Journal of Machine Learning Research 8 (2007) 1027-1061

Submitted 3/06; Revised 12/06; Published 5/07

Dimensionality Reduction of Multimodal Labeled Data
by Local Fisher Discriminant Analysis∗
Masashi Sugiyama

SUGI @ CS . TITECH . AC . JP

Department of Computer Science
Tokyo Institute of Technology
2-12-1, O-okayama, Meguro-ku, Tokyo, 152-8552, Japan

Editor: Sam Roweis

Abstract
Reducing the dimensionality of data without losing intrinsic information is an important preprocessing step in high-dimensional data analysis. Fisher discriminant analysis (FDA) is a traditional
technique for supervised dimensionality reduction, but it tends to give undesired results if samples in a class are multimodal. An unsupervised dimensionality reduction method called localitypreserving projection (LPP) can work well with multimodal data due to its locality preserving
property. However, since LPP does not take the label information into account, it is not necessarily
useful in supervised learning scenarios. In this paper, we propose a new linear supervised dimensionality reduction method called local Fisher discriminant analysis (LFDA), which effectively
combines the ideas of FDA and LPP. LFDA has an analytic form of the embedding transformationand the solution can be easily computed just by solving a generalized eigenvalue problem. We
demonstrate the practical usefulness and high scalability of the LFDA method in data visualization and classification tasks through extensive simulation studies. We also show that LFDA can beextended to non-linear dimensionality reduction scenarios by applying the kernel trick.
Keywords: dimensionality reduction, supervised learning, Fisher discriminant analysis, locality
preserving projection, affinity matrix

1. Introduction
The goal of dimensionality reduction is to embed high-dimensional data samples in a low-dimensional
space so that most of ‘intrinsic information’ contained in the data is preserved (e.g., Roweis andSaul, 2000; Tenenbaum et al., 2000; Hinton and Salakhutdinov, 2006). Once dimensionality reduction is carried out appropriately, the compact representation of the data can be used for various
succeeding tasks such as visualization, classification, etc. In this paper, we consider the supervised
dimensionality reduction problem, that is, samples are accompanied with class labels.
Fisher discriminant analysis (FDA) (Fisher, 1936; Fukunaga, 1990) is a popular method forlinear supervised dimensionality reduction.1 FDA seeks for an embedding transformation such∗. An efficient MATLAB implementation of local Fisher discriminant analysis is available from the author’s website:‘http://sugiyama-www.cs.titech.ac.jp/˜sugi/software/LFDA/’.
1. FDA may refer to the classification method which first projects data samples onto a one-dimensional subspace andthen classifies the samples by thresholding (Fisher, 1936; Duda et al., 2001). The one-dimensional embedding spaceused here is obtained as the maximizer of the so-called Fisher criterion. This Fisher criterion can be used for
dimensionality reduction onto a space with dimension more than one in multi-class problems (Fukunaga, 1990). With
some abuse, we refer to the dimensionality reduction method based on the Fisher criterion as FDA (see Section 2.2
for detail).
c 2007 Masashi Sugiyama.

S UGIYAMA

that the between-class scatter is maximized and the within-class scatter is minimized. FDA is a
traditional but useful method for dimensionality reduction. However, it tends to give undesired
results if samples in a class form several separate clusters (i.e., multimodal) (see, e.g., Fukunaga,
1990).
Within-class multimodality can be observed in many practical applications. For example, indisease diagnosis, the distribution of medial checkup samples of sick patients could be multimodal
since there may be several different causes even for a single disease. In a traditional task of handwritten digit recognition, within-class multimodality appears if digits are classified into, for example, even and odd numbers. More generally, solving multi-class classification problems by a set oftwo-class ‘one-versus-rest’ problems naturally induces within-class multimodality. For this reason,
there is a universal need for reducing the dimensionality of multimodal data.
In order to reduce the dimensionality of multimodal data appropriately, it is important to preserve the local structure of the data. Locality-preserving projection (LPP) (He and Niyogi, 2004)
meets this requirement; LPP seeks for an embedding transformation such that nearby data pairs
in the original space close in the embedding space. Thus LPP can reduce the dimensionality ofmultimodal data without losing the local structure. However, LPP is an unsupervised dimensionality reduction method and does not take the label information into account. Therefore, it does notnecessarily work appropriately in supervised dimensionality reduction scenarios.
In this paper, we propose a new dimensionality reduction method called local Fisher discriminant analysis (LFDA). LFDA effectively combines the ideas of FDA and LPP, that is, LFDA maximizes between-class separability and preserves within-class local structure at the same time. Thus
LFDA is useful for dimensionality reduction of multimodal labeled data.
The original FDA provides a meaningful result only when the dimensionality of the embedding
space is smaller than the number of classes because of the rank deficiency of the between-classscatter matrix (Fukunaga, 1990). This is an essential limitation of FDA in dimensionality reduction.
On the other hand, the proposed LFDA does not generally suffer from this problem and can be
employed for dimensionality reduction into an arbitrary dimensional space. Furthermore, LFDA
inherits an excellent property from FDA—it has an analytic form of the embedding matrix andthe solution can be easily computed just by solving a generalized eigenvalue problem. This is anadvantage over recently proposed supervised dimensionality reduction methods (e.g., Goldberger
et al., 2005; Globerson and Roweis, 2006). Furthermore, LFDA can be naturally extended to nonlinear dimensionality reduction scenarios by applying the kernel trick (Sch ölkopf and Smola, 2002).
The rest of this paper is organized as follows. In Section 2, we formulate the linear dimensionality reduction problem, briefly review FDA and LPP, and illustrate how they typically behave.
In Section 3, we define LFDA and show its fundamental properties. In Section 4, we discuss therelation between LFDA and other methods. In Section 5, we numerically evaluate the performanceof LFDA and existing methods in visualization and classification tasks using benchmark data sets.
Finally, we give concluding remarks and future prospects in Section 6.

2. Linear Dimensionality Reduction
In this section, we formulate the problem of linear dimensionality reduction and review existing
methods.
1028

L OCAL F ISHER D ISCRIMINANT A NALYSIS

2.1 Formulation
Let xi ∈ Rd (i = 1, 2, . . . , n) be d-dimensional samples and yi ∈ {1, 2, . . . , c} be associated classlabels, where n is the number of samples and c is the number of classes. Let n ` be the number ofsamples in class `:c∑ n` = n.

`=1

Let X be the matrix of all samples:X ≡ (x1 |x2 | · · · |xn ).
Let zi ∈ Rr (1 ≤ r ≤ d) be low-dimensional representations of x i , where r is the reduced dimension(i.e., the dimension of the embedding space). Effectively we consider d to be large and r to be small,but not limited to such cases.
For the moment, we focus on linear dimensionality reduction, that is, using a d × r transformation matrix T , the embedded samples zi are given byzi = T > x i ,where > denotes the transpose of a matrix or vector. In Section 3.4, we extend our discussion to thenon-linear dimensionality reduction scenarios where the mapping from x i to zi is non-linear.
2.2 Fisher Discriminant Analysis for Dimensionality Reduction
One of the most popular dimensionality reduction techniques is Fisher discriminant analysis (FDA)(Fisher, 1936; Fukunaga, 1990; Duda et al., 2001). Here we briefly describe the definition of FDA.
Let S(w) and S(b) be the within-class scatter matrix and the between-class scatter matrix:cS(w) ≡ ∑ ∑ (xi − µ` )(xi − µ` )> ,(1)

S(b) ≡ ∑ n` (µ` − µ)(µ` − µ)> ,(2)`=1 i:yi =`c`=1where ∑i:yi =` denotes the summation over i such that yi = `, µ` is the mean of the samples in class`, and µ is the mean of all samples:
µ` ≡

1
xi ,
n` i:y∑
i =`

µ≡

1 n
1 c
x i = ∑ n` µ` .
∑
n i=1
n `=1

We assume that S(w) has full rank. The FDA transformation matrix T FDA is defined as follows:2
h 
i
> (w)
−1 > (b)T FDA ≡ argmax tr (T S T ) T S T .
(3)
T ∈Rd×r

2. The following definition is also used in the literature (e.g., Fukunaga, 1990) and yields the same solution.

det T > S(b) T,
T FDA = argmax  
T ∈Rd×r det T > S(w) T
where det (·) denotes the determinant of a matrix.

1029

S UGIYAMA

That is, FDA seeks a transformation matrix T such that the between-class scatter is ‘maximized’
while the within-class scatter is ‘minimized’. In the above formulation, we implicitly assumed that
T > S(w) T is invertible. This implies that the above optimization is subject torank (T ) = r.
Let {ϕk }dk=1 be the generalized eigenvectors associated with the generalized eigenvalues λ 1 ≥λ2 ≥ · · · ≥ λd of the following generalized eigenvalue problem:
S(b) ϕ = λS(w) ϕ.
Then a solution T FDA of the above maximization problem is analytically given by
T FDA = (ϕ1 |ϕ2 | · · · |ϕr ).
Note that the solution is not unique and the following simple constraint is sometimes imposed
additionally (Fukunaga, 1990).
(w)
T>
T FDA = I r ,FDA S
where I r is the identity matrix on Rr . This constraint makes the within-class scatter in the embedding
space sphered.
The between-class scatter matrix S(b) has at most rank c − 1 (Fukunaga, 1990). This implies
that the multiplicity of λ = 0 is at least d − c + 1. Therefore, FDA can find at most c − 1 meaningful
features; the remaining features found by FDA are arbitrary. This is an essential limitation of FDA
for dimensionality reduction and is very restrictive in practice.
2.3 Locality-Preserving Projection
Another dimensionality reduction technique that is relevant to the current setting is locality-preserving
projection (LPP) (He and Niyogi, 2004). Here we review LPP.
Let A be an affinity matrix, that is, the n-dimensional matrix with the (i, j)-th element A i, j being
the affinity between xi and x j . We assume that Ai, j ∈ [0, 1]; Ai, j is large if xi and x j are ‘close’and Ai, j is small if xi and x j are ‘far apart’. There are several different manners of defining A. We
briefly describe typical definitions in Appendix D. The LPP transformation matrix T LPP is defined
as follows:3!
1 n
T LPP ≡ argmin∑ Ai, j kT > xi − T > x j k2
2
d×r
T ∈R
i, j=1subject to T > XDX > T = I r ,(4)where D is the n-dimensional diagonal matrix with i-th diagonal element beingnDi,i ≡ ∑ Ai, j .
j=1

3. The matrix D in the constraint (4) is motivated by a geometric argument (Belkin and Niyogi, 2003). However, it issometimes dropped for the sake of simplicity (Ham et al., 2004).

1030

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Eq. (4) implies that LPP looks for a transformation matrix T such that nearby data pairs in theoriginal space Rd are kept close in the embedding space. The constraint (4) is imposed for avoiding
degeneracy.
Let {ψk }dk=1 be the generalized eigenvectors associated with the generalized eigenvalues γ 1 ≥γ2 ≥ · · · ≥ γd of the following generalized eigenvalue problem:
XLX > ψ = γXDX > ψ,
whereL ≡ D − A.
L is called the graph-Laplacian matrix in the spectral graph theory (Chung, 1997), where A is seenas the adjacency matrix of a graph. He and Niyogi (2004) showed that a solution of Eq. (4) is givenby
T LPP = (ψd |ψd−1 | · · · |ψd−r+1 ).
2.4 Typical Behavior of FDA and LPP
Dimensionality reduction results obtained by FDA and LPP are illustrated in Figure 1 (LFDA will
be defined and explained in Section 3)—two-dimensional two-class data samples are embedded intoa one-dimensional space. In LPP, the affinity matrix A is determined by the local scaling method(Zelnik-Manor and Perona, 2005, see also Appendix D.4).
For the simplest data set depicted in Figure 1(a), both FDA and LPP nicely separate the samples
in different classes (‘◦’ and ‘×’) from each other. For the data set depicted in Figure 1(b), FDA
still works well, but LPP mixes samples in different classes into a single cluster. This is caused by
the unsupervised nature of LPP. On the other hand, for the data set depicted in Figure 1(c), LPP
works well but FDA collapses the samples in different classes into a single cluster. The reason for
the failure of FDA is that the ‘levels’ of the between-class scatter and the within-class scatter are
not evaluated in an intuitively natural way because of the two separate clusters in ‘◦’-class (see also
Fukunaga, 1990).

3. Local Fisher Discriminant Analysis
As illustrated in Figure 1, FDA can perform poorly if samples in a class form several separate
clusters (i.e., multimodal). In other words, the undesired behavior of FDA is caused by the globality
when evaluating the within-class scatter and the between-class scatter (e.g., Figure 1(c)). On theother hand, because of the unsupervised nature of LPP, it can overlap samples in different classes ifthey are close in the original high-dimensional space Rd (e.g., Figure 1(b)).
To overcome these problems, we propose combining the ideas of FDA and LPP; more specifically, we evaluate the levels of the between-class scatter and the within-class scatter in a local
manner. This allows us to attain between-class separation and within-class local structure preservation at the same time. We call our new method local Fisher discriminant analysis (LFDA).
3.1 Reformulating FDA
In order to introduce LFDA, let us first reformulate FDA in a pairwise manner.
1031

S UGIYAMA

FDA
LPP
LFDA

10

5

0

−5

−10

−10

−5

0

5

10

(a) Toy data set 1

FDA
LPP
LFDA

10

5

5

0

0

−5

−5

−10

−10

−5

0

5

10

FDA
LPP
LFDA

10

−10

(b) Toy data set 2

−10

−5

0

5

10

(c) Toy data set 3

Figure 1: Examples of dimensionality reduction by FDA, LPP and LFDA. Two-dimensional twoclass samples are embedded into a one-dimensional space. The line in the figure denotes
the one-dimensional embedding space (which the data samples are projected on) obtainedby each method.

1032

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Lemma 1 S(w) and S(b) defined by Eqs. (1) and (2) can be expressed as
S(w) =

1 n(w)
Wi, j (xi − x j )(xi − x j )> ,
2 i,∑
j=1

(5)

S(b) =

1 n(b)
Wi, j (xi − x j )(xi − x j )> ,
2 i,∑
j=1

(6)

(

(7)

where
(w)
Wi, j ≡(b)

Wi, j ≡(

1/n` if yi = y j = `,0

if yi 6= y j ,1/n − 1/n` if yi = y j = `,1/nif yi 6= y j .

(8)

A proof of Lemma 1 is given in Appendix A. Note that 1/n − 1/n ` in Eq. (8) is negative while
1/n` and 1/n in Eqs. (7) and (8) are positive. This implies that if the data pairs in the same class are
made close, the within-class scatter matrix S (w) gets ‘small’ and the between-class scatter matrix S (b)gets ‘large’. On the other hand, if the data pairs in different classes are separated from each other,the between-class scatter matrix S(b) gets ‘large’. Therefore, we may interpret FDA as keeping thesample pairs in the same class close and the sample pairs in different classes apart. A more formal
discussion on the above interpretation is given in Appendix B.
3.2 Definition and Typical Behavior of LFDA(w)
Based on the above pairwise expression, let us define the local within-class scatter matrix e
S and
(b)
the local between-class scatter matrix e
S as follows.
n

(w)
1
e
e (w) (xi − x j )(xi − x j )> ,S ≡ ∑W
2 i, j=1 i, j

(9)

n

where

(b)
1
e
e (b) (xi − x j )(xi − x j )> ,S ≡ ∑W
2 i, j=1 i, je (w) ≡
Wi, je (b) ≡
W
i, j

(

(

Ai, j /n` if yi = y j = `,0

if yi 6= y j ,Ai, j (1/n − 1/n` ) if yi = y j = `,1/nif yi 6= y j .

(10)
(11)

Namely, according to the affinity Ai, j , we weight the values for the sample pairs in the same class.
(w)(b)
This means that far apart sample pairs in the same class have less influence on e
S and eS . Notethat we do not weight the values for the sample pairs in different classes since we want to separate
them from each other irrespective of the affinity in the original space. From here on, we denote thelocal counterparts of matrices by symbols with tilde.
1033

S UGIYAMAWe define the LFDA transformation matrix T LFDA as
h 
i
(w)
(b)
T LFDA ≡ argmax tr (T >eS T )−1 T >e
S T .

(12)

T ∈Rd×r

That is, we look for a transformation matrix T such that nearby data pairs in the same class are
made close and the data pairs in different classes are separated from each other; far apart data pairs
in the same class are not imposed to be close.
Eq. (12) is of the same form as Eq. (3). Therefore, we can similarly compute an analytic form(b)(w)
of T LFDA by solving a generalized eigenvalue problem of e
S and eS . An efficient implementationof LFDA is summarized as a pseudo code in Figure 2 (see Appendix C for detail).
Toy examples of dimensionality reduction by LFDA are illustrated in Figure 1. We used the
local scaling method for computing the affinity matrix A (see Appendix D.4). Note that we performthe nearest neighbor search in the local scaling method in a classwise manner since we do notneed the affinity values for the sample pairs in different classes (see Eqs. 10 and 11). This highly
contributes to reducing the computational cost (see Appendix C). Figure 1 shows that LFDA gives
desirable results for all three data sets, that is, LFDA can compensate for the drawbacks of FDA and
LPP by effectively combining the ideas of FDA and LPP.
If the affinity value Ai, j is set to 1 for all sample pairs (i.e., all pairs are ‘equally close’ to each(w)(b)
other), e
S and eS agree with S(w) and S(b) , respectively, and LFDA is reduced to the original
FDA. Therefore, LFDA may be regarded as a natural localized variant of FDA.
3.3 Properties of LFDA
Here we discuss fundamental properties of LFDA.
(w)

First, we give an interpretation of LFDA in terms of the ‘pointwise scatter’. eSpressed as(w)
1 n 1 e(w)e,
S = ∑ P2 i=1 nyi i

can be ex-e(w)where nyi is the number of samples in the class to which the sample x i belongs and P
is theipointwise local within-class scatter matrix around xi :(w)

Therefore, ‘minimizing’ eSei(w) ≡P

∑ Ai, j (x j − xi )(x j − xi )> .

j:y j =yicorresponds to minimizing the weighted sum of the pointwise local(b)
within-class scatter matrices over all samples. e
S can also be expressed in a similar way as(b)where Pi



(b)
1 n 1
1 e(w) 1 n (b)
e
S = ∑
Pi +
−
∑ Pi ,
2 i=1 n nyi
2n i=1is the pointwise between-class scatter matrix around x i :
(b)

Pi ≡

∑ (x j − xi )(x j − xi )> .

j:y j 6=yi

1034

(13)

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Input:
Output:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:

Labeled samples {(xi , yi ) | xi ∈ Rd , yi ∈ {1, 2, . . . , c}}ni=1
Dimensionality of embedding space r (1 ≤ r ≤ d)d × r transformation matrix T LFDA(b)e
S ←− 0d×d ;(w)eS ←− 0d×d ;for ` = 1, 2, . . . , c% Compute scatter matrices in a classwise mannern`{xi }i=1 ←− {x j } j:y j =` ;for i = 1, 2, . . . , n`
% Determine local scaling(7)`xi ←− 7th nearest neighbor of xi among {x j }nj=1;(7)

σi ←− kxi − xi k;
end
for i, j = 1, 2, . . . , n`
% Define affinity matrix
Ai, j ←− exp(−kxi − x j k2 /(σi σ j ));end
X ←− (x1 |x2 | · · · |xn` );G ←− Xdiag(A1n` )X > − X A X > ;
(b)
(b)
e
S ←− eS + G/n + (1 − n` /n)X X > + X1n (X1n )> /n;
`

(w)
(w)
e
S ←− eS + G/n` ;`15:
16: end
(b)
(b)
(w)
17: e
S ←− eS − X1n (X1n )> /n − eS ;e }r ←− generalized eigenvalues and normalized eigenvectors of
18: {eλk , ϕ
k k=1

(w)
(b)
e
e;
e = eλe% eλ1 ≥ eλ2 ≥ · · · ≥ eλd
S ϕ
S qϕ
q
q
e1 | eλ2 ϕ
e2 | · · · | eλr ϕer );19: T LFDA = ( eλ1 ϕFigure 2: Efficient implementation of LFDA (see Appendix C for detail). The affinity matrix is
computed by the local scaling method (see Appendix D.4). Matrices and vectors denotedwith underline are classwise counterparts of the original ones. 0 d×d denotes the d × dmatrix with zeros, 1n` denotes the n` -dimensional vector with ones, and diag(A1n` ) denotes the diagonal matrix with diagonal elements A1n` . The generalized eigenvectors inline 18 are normalized by Eq. (14), which is often automatically carried out by an eigensolver. The weighting scheme of the eigenvectors in line 19 is explained in Section 3.3.
A possible bottleneck of the above implementation is the nearest neighbor search in line
6. This could be alleviated by incorporating the prior knowledge of the data structure orby approximation (see Saul and Roweis, 2003, and references therein). Another possible
bottleneck is the computation of X A X > in line 13, which could be eased by sparsely
defining the affinity matrix (see Appendix D). A MATLAB implementation is available
from ‘http://sugiyama-www.cs.titech.ac.jp/˜sugi/software/LFDA/’.

1035

S UGIYAMA(b)(b)
Note that Pi does not include the localization factor Ai, j . Eq. (13) implies that ‘maximizing’ eScorresponds to minimizing the weighted sum of the pointwise local within-class scatter matrices
and maximizing the sum of the pointwise between-class scater matrices.
Next, we discuss the issue of eigenvalue multiplicity in LFDA. The original FDA allows us to
extract at most c − 1 meaningful features since the between-class scatter matrix S (b) has rank at most(b)c − 1 (Fukunaga, 1990). On the other hand, the local between-class scatter matrix e
S generally has
a much higher rank with less eigenvalue multiplicity, thanks to the localization factor A i, j included(b)e (b) (see Eq. 11). In the simulation shown in Section 5, e
in W
S is always full rank for various data
sets. Therefore, the proposed LFDA can be practically employed for dimensionality reduction into
any dimensional spaces. This is a very important and significant improvement over the original
FDA.
Finally, we discuss the invariance property of LFDA. The value of the LFDA criterion (12) isinvariant under linear transformations, that is, for any r-dimensional invertible matrix H, T LFDA H isalso a solution of Eq. (12). Therefore, the solution T LFDA is not unique—the range of the transformation H > T >LFDA is uniquely determined, but the distance metric (Goldberger et al., 2005; Globerson and Roweis, 2006; Weinberger et al., 2006) in the embedding space can be arbitrary becauseof the arbitrariness of the matrix H. In practice, we propose determining the LFDA transformationek }dk=1 so thatmatrix T LFDA as follows. First, we rescale the generalized eigenvectors { ϕ
(w)

ek e
ϕ
S

e k0 =
ϕ(

1 if k = k0 ,0 if k 6= k0 .

(14)

Note that this rescaling is often automatically carried out by an eigensolver. Then we weight eachgeneralized eigenvector by the square root of its associated generalized eigenvalue, that is,
q
q
q
e
e
e
e
er ),
(15)
T LFDA = ( λ1 ϕ1 | λ2 ϕ2 | · · · | eλr ϕwhere eλ1 ≥ eλ2 ≥ · · · ≥ eλd . This weighting scheme weakens the influence of minor eigenvectors andis shown to work well in experiments (see Section 5).
3.4 Kernel LFDA for Non-Linear Dimensionality Reduction
Here we show how LFDA can be extended to non-linear dimensionality reduction scenarios.
As detailed in Appendix C, the generalized eigenvalue problem that needs to be solved in LFDA
can be expressed ase(b) X > ϕe(w) X > ϕe = eλX L
e,
XL
(16)

e(b) = Le(m) − Le(w) and Le(m) and Le(w) are defined by Eqs. (33) and (35), respectively. Since
where Le in Eq. (16) belongs to the range of X > , it can be expressed by using some vector αe ∈ Rn as
X >ϕe = X >X α
e = Kαe,
X >ϕwhere K is the n-dimensional matrix with the (i, j)-th element being
Ki, j ≡ x>
i x j.
1036

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Then multiplying Eq. (16) by X > from the left-hand side yields
e(b) K αe(w) K αe = eλK L
e.
KL

(17)

This implies that {xi }ni=1 appear only in terms of their inner products. Therefore, we can obtaina non-linear variant of LFDA by the kernel trick (Vapnik, 1998; Sch ölkopf et al., 1998), which is
explained below.
Let us consider a non-linear mapping φ(x) from Rd to a reproducing kernel Hilbert space H(Aronszajn, 1950). Let K(x, x0 ) be the reproducing kernel of H . A typical choice of the kernel
function would be the Gaussian kernel:kx − x0 k2,
K(x, x0 ) = exp −2σ2with σ > 0. For other choices, see, for example, Wahba (1990), Vapnik (1998), and Sch ölkopf andSmola (2002). Because of the reproducing property of K(x, x 0 ), K is now the kernel matrix, that is,the (i, j)-th element is given by
Ki, j = hφ(xi ), φ(x j )i = K(xi , x j ),where h·, ·i denotes the inner product in H .
e(w) is always degenerated (since Le(w) (1, 1, . . . , 1)> always vanishes;
It can be confirmed that Le(w) K is always degenerated and we cannot directly solve thesee Eq. 35 for detail). Therefore, K Le(w) K
generalized eigenvalue problem (17). To cope with this problem, we propose regularizing K L
and solving the following generalized eigenvalue problem instead (cf. Friedman, 1989).
e(b) K αe(w) K + εI n )αe = eλ(K L
e,
KL

(18)

e k }nk=1 be the generalized eigenvectors associated with the generwhere ε is a small constant. Let {αalized eigenvalues eλ1 ≥ eλ2 ≥ · · · ≥ eλn of Eq. (18). Then the embedded image of φ(x0 ) in H is givenby
K(x1 , x0 )qqq
K(x2 , x0 )

>
e
e
e
e
e
e
( λ 1 α1 | λ 2 α2 | · · · | λ r αr ) 
.
..


.
K(xn , x0 )We call this kernelized variant of LFDA kernel LFDA (KLFDA).
Recently, kernel functions for non-vectorial structured data such as strings, trees, and graphs
have been proposed (see, e.g., Lodhi et al., 2002; Duffy and Collins, 2002; Kashima and Koyanagi,
2002; Kondor and Lafferty, 2002; Kashima et al., 2003; Gärtner et al., 2003; Gärtner, 2003). Since
KLFDA uses the samples only via the kernel function K(x, x 0 ), it allows us to reduce the dimensionality of such non-vectorial data.

4. Comparison with Related Methods
In this section, we discuss the relation between the proposed LFDA and other methods.
1037

S UGIYAMA

4.1 Dimensionality Reduction Using Local Discriminant Information
A discriminant adaptive nearest neighbor (DANN) classifier (Hastie and Tibshirani, 1996a) employs an adapted distance metric at each test point for classification. Based on a similar idea, theyalso proposed a global supervised dimensionality reduction method using local discriminant information (LDI) in the same paper. We refer to this supervised dimensionality reduction method as
LDI. The main idea of LDI is to localize FDA—which is very similar to the proposed LFDA. Here
we discuss the relation between LDI and LFDA.
In LDI, the data samples {xi }ni=1 are first sphered according to the within-class scatter matrix(w)
S , that is, for i = 1, 2, . . . , n,1xi ≡ (S(w) )− 2 xi .
Let Ai, j be the weight of sample x j around xi defined by " 3 #3kxi −x j k(K)
if kxi − x j k < kxi − xi k,
1−(K)
Ai, j ≡
kxi −xi k0otherwise.
(K)

where xi is the K-th nearest neighbor of xi in the sphered space. Note that 0 ≤ Ai, j ≤ 1 and Ai, j
is non-increasing as kxi − x j k increases. Thus it has the same meaning as our affinity matrix. K issuggested to be determined by
K = max(n/5, 50).
[i]

Let µ` be the local weighted mean of the sphered samples in class ` around x i , and let µ[i] be thelocal weighted mean of the sphered samples around x i :
[i]

µ` ≡
µ[i] ≡

1[i]∑ Ai, j x j ,n` j:y j =`1 n
1 c [i] [i]Ax=∑ i, j j n[i] ∑ n` µ` ,n[i] j=1
`=1

where

∑ Ai, j ,

[i]

n` ≡

j:y j =`
nn[i] ≡ ∑ Ai, j .
j=1

Let S

(b)

be the average between sum-of-squares matrix defined as
S

(b)

n

1 c [i] [i][i]n` (µ` − µ[i] )(µ` − µ[i] )> .
∑[i]
i=1 n `=1≡∑

The LDI transformation matrix T LDI is defined as

h
i
> (b)
T LDI ≡ argmax T S T
T ∈Rd×r>

subject to T T = I r .
1038

L OCAL F ISHER D ISCRIMINANT A NALYSIST LDI is a transformation matrix for sphered samples; the LDI transformation matrix T LDI for nonsphered samples is given by
1
T LDI = (S(w) )− 2 T LDI .
Similar to FDA (and LFDA), T LDI can be efficiently computed by solving a generalized eigenvalue
problem.
(b)
The average between sum-of-squares matrix S is conceptually very similar to the local between(b)(b)
class scatter matrix e
S in LFDA. Indeed, as proved in Appendix E, we can express S in a pairwisemanner as
1 n
(b)
(b)
S = ∑ W i, j (xi − x j )(xi − x j )> ,(19)
2 i, j=1
where

!

n

1
1
1


−
Ai,k A j,k if yi = y j = `, ∑ [k](b)n[k] n`[k]
k=1 n
W i, j ≡n1
if yi 6= y j .
∑ (n[k] )2 Ai,k A j,k
k=1

(20)

However, there exist critical differences between LDI and LFDA. A significant difference is that thevalues for the sample pairs in different classes are also localized in LDI (see Eq. 20), while they are
kept unlocalized in LFDA (see Eq. 11). This implies that far apart sample pairs in different classes
could be made close in LDI, which is not desirable in supervised dimensionality reduction. Further(b)(b)(b)
more, the computation of S is slightly less efficient than e
S since W includes the summation
over k.
Another important difference between LDI and LFDA is that the within-class scatter matrix S (w)is not localized in LDI. However, as we showed in Section 3.1, the within-class scatter matrix S (w)also accounts for collapsing the within-class multimodal structure (i.e., far apart sample pairs in thesame class are made close). This phenomenon is experimentally confirmed in Section 5.2.
4.2 Mixture Discriminant Analysis
FDA can be interpreted as maximum likelihood estimation of Gaussian distributions with common
covariance and different means for each class. Based on this view, Hastie and Tibshirani (1996b)
proposed mixture discriminant analysis (MDA), which extends FDA to maximum likelihood estimation of Gaussian mixture distributions.
A maximum likelihood solution is obtained by an EM-type algorithm (cf. Dempster et al., 1977).
However, this is an iterative algorithm and gives only a local optimal solution. Therefore, the computation of MDA is rather slow and there is no guarantee that the global solution can be obtained.
Furthermore, the number of mixture components (clusters) in each class as well as the initial location of cluster centers should be determined by users. For cluster centers, using standard techniquessuch as k-means clustering (MacQueen, 1967; Everitt et al., 2001) or learning vector quantization(Kohonen, 1989) are recommended. However, they are also iterative algorithms and have no guarantee that the global solution can be obtained. Furthermore, there seems to be no systematic methodfor determining the number of clusters.
On the other hand, the proposed LFDA contains no tuning parameters (given that the affinity
matrix is determined by the local scaling method, see Appendix D.4) and the global solution can1039

S UGIYAMA

be obtained analytically. However, it still lacks a probabilistic interpretation, which remains open
currently.
4.3 Neighborhood Component Analysis
Goldberger et al. (2005) proposed a supervised dimensionality reduction method called neighborhood component analysis (NCA). The NCA transformation matrix T NCA is defined as follows.
!
n

T NCA ≡ argmax
T ∈Rd×r

where
pi, j (U) ≡







∑ ∑ pi, j (T T > ) ,(21)

i=1 j:y j =yi
exp −(xi − x j )>U(xi − x j )∑k6=i exp {−(xi − xk )>U(xi − xk )}if i 6= j,0

if i = j.

(22)

The above definition corresponds to maximizing the expected number of correctly classified samples
by a stochastic variant of nearest neighbor classifiers. Therefore, NCA seeks a transformation matrix
T such that the between-class separability is maximized.
Eqs. (21) and (22) imply that nearby data pairs in the same class are made close, which issimilar to the proposed LFDA. Indeed, the simulation results in Section 5.2 show that NCA tends
to preserve the multimodal structure of the data very well. However, a crucial weakness of NCA
is optimization: the optimization problem (21) is non-convex. Therefore, there is no guarantee thatthe globally optimal solution can be obtained. Goldberger et al. (2005) proposed using a gradientascent method for optimization:
T ← T + ε∇JNCA (T ),(23)
where ε (> 0) is the step size and the gradient ∇JNCA (T ) is given by()(
n∇JNCA (T ) = 2T ∑∑ pi, j (T T > )j:y j =yi

i=1

−

n

∑ pi, j (T T > )(xi − x j )(xi − x j )>j=1

∑ pi, j (T T )(xi − x j )(xi − x j )>j:y j =yi

>

!

)

.

The gradient ascent iteration (23) is computationally rather inefficient. Also, the choice of the step
size ε is troublesome. If the step size is small enough, the convergence to one of the local optima
is guaranteed but such a choice makes the convergence very slow; on the other hand, if the step
size is too large, gradient flows oscillate and proper convergence properties may not be guaranteed
anymore. Furthermore, the choice of the termination condition in the iterative algorithm is often
cumbersome in practice.
Because of the non-convexity of the optimization problem, the quality of the obtained solution
depends on the initialization of the matrix T . A useful heuristic to alleviate the local optimum
problem is to employ the FDA (or LFDA) result as an initial matrix for optimization (Goldberger
et al., 2005). In the experiments in Section 5, using the LFDA result as an initial matrix appears tobe better than the random initialization. However, the local optima problem still remains even withthe above heuristic.
1040

L OCAL F ISHER D ISCRIMINANT A NALYSIS

When a dimensionality reduction technique is applied to classification tasks, we often want toembed the data samples into spaces with several different dimensions—the best dimensionality islater chosen by, for example, cross-validation (Stone, 1974; Wahba, 1990). In such a scenario,
NCA requires to optimize the transformation matrix individually for each dimensionality r of theembedding space. On the other hand, LFDA needs to compute the transformation matrix only once
for the largest r; its sub-matrices become the optimal solutions for smaller dimensions. Therefore,
LFDA is computationally more efficient than NCA in this scenario.
A simple MATLAB implementation of NCA is available.4 We use this software in Section 5.
4.4 Maximally Collapsing Metric Learning
In order to overcome the computational problem of NCA, Globerson and Roweis (2006) proposedan alternative method called maximally collapsing metric learning (MCML).
Let p∗i, j be the ‘ideal’ value of pi, j (U) defined by Eq. (22):p∗i, j ∝where p∗i, j is normalized so that(

1 if yi = y j ,0 if yi 6= y j ,∑ p∗i, j = 1.
j6=ip∗i, j can be attained if all samples in the same class collapse into a single point while samples
in other classes are mapped to other locations. In reality, however, any U may not be able toattain pi, j (U) = p∗i, j exactly; instead the optimal approximation to p∗i, j under the Kullback-Leibler
divergence (Kullback and Leibler, 1951) is obtained. This is formally defined as
!
∗
n
p
i, j
U MCML ≡ argmin ∑ p∗i, j log
p
(U)
d×d
i,
j
U∈R
i, j=1subject to U ∈ PSD(r),(24)

where PSD(r) is the set of all positive semidefinite matrices of rank r (i.e., r eigenvalues are positive and others are zero). Once U MCML is obtained, the MCML transformation matrix T MCML is
computed by
T MCML = (φ1 |φ2 | · · · |φr ),(25)
where {φk }rk=1 are the eigenvectors associated with the positive eigenvalues η 1 ≥ η2 ≥ · · · ≥ ηr > 0of the following eigenvalue problem:U MCML φ = ηφ.
One of the motivations of MCML is to alleviate the difficulty of optimization in NCA. However,
MCML still has a weakness in optimization: the optimization problem (24) is convex only whenr = d, that is, the dimensionality is not reduced but only the distance metric of the original space ischanged. This means that if r < d (which is our primal focus in this paper), we may not be able to
4. Implementation available at ‘http://www.cs.berkeley.edu/˜fowlkes/software/nca/’.

1041

S UGIYAMA

obtain the globally optimal solution. Globerson and Roweis (2006) proposed the following heuristic
algorithm to approximate T MCML .
First, the optimization problem (24) with r = d is solved:
!
∗
n
p
i,
j
b MCML ≡ argmin ∑ p∗i, j log
Upi, j (U)U∈Rd×d i, j=1
subject to U ∈ PSD(d).

(26)

b MCML is not known
Although Eq. (26) is convex, an analytic form of the unique optimal solution U
yet. Globerson and Roweis (2006) proposed using the following alternate iterative procedure forb MCML .
obtaining U
U ← U − ε∇JMCML (U),

(27)

d

>
U ← ∑ max(0, ηbk )b
φk b
φk ,

(28)

k=1

bk and b
where ε (> 0) is the step size, η
φk are eigenvalues and eigenvectors of U, and the gradient∇JMCML (U) is given byn∇JMCML (U) = ∑ (p∗i, j − pi, j (U))(xi − x j )(xi − x j )> .
i, j=1b MCML is carried out and eigenvalues η
b1 ≥ η
b2 ≥ · · · ≥ ηbdThen the eigenvalue decomposition of Udband associated eigenvectors {φk }k=1 are obtained:b MCML b
bb
U
φ=η
φ.

φk }rk=1 , which yields
Finally, {φk }rk=1 in Eq. (25) are replaced by {bT MCML ≈ (bφ1 |bφ2 | · · · |b
φr ).

(29)

This approximation is shown to be practically useful (Globerson and Roweis, 2006), although there
seems to be no theoretical analysis for this approximation.
MCML may have an advantage over NCA in computation: there exists the analytic approximation (29) that can be computed efficiently using the solution of another convex optimization
problem (26). However, MCML still relies on the gradient-based alternate iterative algorithm (27)–(28) to solve the convex optimization problem (26), which is computationally very expensive since
the eigenvalue decomposition of a d-dimensional matrix should be carried out in each iteration (seeEq. 28). Furthermore, the difficulty of appropriately choosing the step size and the termination
condition in the iterative procedure still remains.
Since MCML requires all the samples in the same class to collapse into a single point, it is
not necessarily useful in dimensionality reduction of multimodal data samples. Furthermore, the
MCML results can be significantly influenced by outliers since the outliers are also required tocollapse into the same single point together with other samples. This phenomenon is illustrated in
Figure 3, where a single outlier significantly changes the MCML result.
Globerson and Roweis (2006) showed that the sufficient statistics of the MCML algorithm arepointwise scatter matrices (cf. Section 3.3). Since LFDA also has an interpretation in terms ofpointwise scatter matrices, there may be a link between LFDA and MCML and this needs to beinvestigated in the future work.
1042

L OCAL F ISHER D ISCRIMINANT A NALYSIS

10

10
LFDA
MCML

8
6

6

4

4

2

2

0

0

−2

−2

−4

−4

−6

−6

−8

−8

−10
−10

0

10

20

30

40

50

LFDA
MCML

8

−10
−10

(a) Toy data set 1

outlier↑

0

10

20

30

40

50

(b) Toy data set 1’

Figure 3: Toy examples of dimensionality reduction. The toy data set 1 is equivalent to the oneused in Figure 1(a). The data set 1’ includes a single outlier.

4.5 Remark on Rank Constraint
The optimization problem of MCML (see Eq. 24) is not generally convex since the rank constraint
is non-convex (Boyd and Vandenberghe, 2004). The non-convexity induced by the rank constraint
seems to be a universal problem in dimensionality reduction. NCA eliminates the rank constraint
by decomposing U into T T > (see Eqs. 21 and 22). However, even with this decomposition, theoptimization problem is still non-convex. On the other hand, FDA, LDI, and LFDA cast the optimization problem in the form of the Rayleigh quotient. This is computationally very advantageous
since it allows us to analytically determine the range of the embedding space. However, we cannot
determine the distance metric in the embedding space since the Rayleigh quotient is invariant underlinear transformations. For this reason, an additional criterion is needed to determine the distancemetric (see also Section 3.3).

5. Numerical Examples
In this section, we numerically evaluate the performance of LFDA and existing methods.
5.1 Exploratory Data Analysis
Here we use the Thyroid disease data set available from the UCI machine learning repository (Blake
and Merz, 1998) and illustrate how LFDA can be used for exploratory data analysis.
The original data consists of 5-dimensional input vector x of the following laboratory tests.
1. T3-resin uptake test.
2. Total Serum thyroxin as measured by the isotopic displacement method.
1043

S UGIYAMA

8Hyperthyroidism
Hypothyroidism

8

Hyperthyroidism
Hypothyroidism

6
6
4

4

2

2
0

−25

−20

−15
First Feature

−10

−5

20

Euthyroidism

30

0

25

3

4

5
First Feature

6

7

5First Feature

6

7

Euthyroidism

15

20
10

15
10

5
5
0

−25

−20

−15
First Feature

−10

−5

0

(a) FDA34(b) LFDA

Figure 4: Histograms of the first feature values obtained by FDA and LFDA for the Thyroid disease
data set. The top row corresponds to the sick patients and the bottom row corresponds tothe healthy patients.

3. Total Serum triiodothyronine as measured by radioimmuno assay.
4. Basal thyroid-stimulating hormone (TSH) as measured by radioimmuno assay.
5. Maximal absolute difference of TSH value after injection of 200 micro grams of thyrotropinreleasing hormone as compared to the basal value.
The task is to predict whether patients’ thyroids are euthyroidism, hypothyroidism, or hyperthyroidism (Coomans et al., 1983), that is, whether patients’ thyroids are normal, hypo-functioning,or hyper-functioning (Blake and Merz, 1998). The diagnosis (the class label) is based on a complete medical record, including anamnesis, scan etc. Here we merge the hypothyroidism class andthe hyperthyroidism class into a single class and create binary labeled data (whether thyroids arenormal or not). Our goal is to predict whether patients’ thyroids are normal, hypo-functioning, or
hyper-functioning from the binary labeled data samples.
Figure 4 depicts the histograms of the first feature values obtained by FDA and LFDA—the
top row corresponds to the sick patients and the bottom row corresponds to the healthy patients.
This shows that both FDA and LFDA separate the patients with normal thyroids from sick patients
reasonably well. In addition to between-class separability, LFDA clearly preserves the multimodal
structure among sick patients (i.e., hypo-functioning and hyper-functioning), which is lost by ordinary FDA. Another interesting finding from the figure is that the first feature values obtained by
LFDA has a strong negative correlation to the functioning level of thyroids—this could be used forpredicting the functioning level of thyroids.
1044

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Data Set
Letter recognition
Iris

d
16
4

‘◦’-and-‘•’ class‘A’ & ‘C’
‘Setosa’ & ‘Virginica’‘×’ class‘B’‘Versicolour’

Table 1: Two-class data sets used for visualization experiments (r = 2).
5.2 Data Visualization
Here we apply the proposed and existing dimensionality reduction methods to benchmark data sets
and investigate how they behave in data visualization tasks.
We use the Letter recognition data set and the Iris data set available from the UCI machinelearning repository (Blake and Merz, 1998). Table 1 describes the specifications of the data sets.
Each data set contains three types of samples specified by ‘◦’, ‘•’, and ‘×’. We merged ‘◦’ and‘•’ into a single class and created two-class problems. We test LFDA, FDA, LPP, LDI, NCA, and
MCML and evaluate the between-class separability (i.e., ‘◦’ and ‘•’ are well separated from ‘×’)and the within-class multimodality preservation capability (i.e., ‘◦’ and ‘•’ are well grouped). For
LPP and LFDA, we determined the affinity matrix by the local scaling method (see Appendix D.4).
For NCA, we used the LFDA result as an initial matrix since this initialization scheme appears towork better than the random initialization. FDA allows us to extract only one meaningful feature in
two-class classification problems (see Section 2.2), so we choose the second feature randomly here.
Figures 5 and 6 depict the samples embedded in the two-dimensional space found by each method.
The horizontal axis is the first feature found by each method, while the vertical axis is the second
feature.
First, we compare the embedding results of LFDA with those of FDA and LPP. For the Letter recognition data set (see the top row of Figure 5), LFDA nicely separates samples in different
classes from each other, and at the same time, it clearly preserves within-class multimodality. FDA
separates ‘◦’ and ‘•’ from ‘×’ well, but within-class multimodality is lost, that is, ‘◦’ and ‘•’ aremixed. LPP gives two separate clusters of samples, but samples in different classes are mixed in oneof the clusters. For the Iris data set (see the top row of Figure 6), LFDA simultaneously achieves
between-class separation and within-class multimodality preservation. On the other hand, FDA
tends to mix samples in different classes, which would be caused by within-class multimodality.
LPP also works well for this data set because three clusters are well separated from each other inthe original high-dimensional space. Overall, LFDA is found to be more appropriate for embedding labeled multimodal data samples than FDA and LPP, implying that our primal goal has beensuccessfully achieved.
Next, we compare the results of LFDA with those of LDI, NCA, and MCML. For the Letter
recognition data set (see Figure 5), LFDA, LDI, NCA, and MCML separate the samples in different
classes from each other very well. However, LDI and MCML collapse ‘◦’ and ‘•’ into a single
cluster, while LFDA and NCA preserve the multimodal structure clearly. The NCA result is almost
identical to the LFDA result (i.e., the initial value of the NCA iteration), but the result may vary if theinitial value for the gradient ascent algorithm is changed. For the Iris data set (see Figure 6), LFDA,
LDI, and NCA work excellently in both between-class separation and within-class multimodality
preservation. On the other hand, MCML mixes the samples in different classes. Overall, LDI works
fairly well, but the within-class multimodal structure is sometimes lost since LDI only partiallytakes within-class multimodality into account (see Section 4.1). NCA also works very well, which
1045

S UGIYAMA

LFDA

FDA

30

A
C
B

28

LPP

12
6

11

26

10

24

9

22

8

4
2
0

20

−2

7

18

−4

6

16

−6

5

14

−8

4

−10

3

−12

12
10

2
−10

0

10

20

−14

−10

−5

0

5

−40

−35

NCA

LDI
0.2

30

0.1

25

−30

−25

0

1

MCML
−0.4
−0.5
−0.6

0

20
−0.7

−0.1

15

−0.8

−0.2

−0.9

10

−0.3

−1
−0.6

−0.4

−0.2

0

0.2

−20

−10

0

10

20

−3

−2

−1

Figure 5: Visualization of the Letter recognition data set.
LFDA

LPP

FDA

55

3

−1.5

2
50
1

−2

0

45
−2.5

40

−1
−2

−3

−3
35

SetosaVirginica
Verisicolour
0

50

100

150

−4

−3.5

−5

200

−4

−3.5

−3

−2.5

−2

−1.5

−14

NCA

LDI

−12

−10

−8

−6

MCML

0.5

0.2

0

0.45

0

0.4
0.35

−5

−0.2

0.3
−0.4

0.25

−10

0.2

−0.6

0.15

−15

−0.8

0.1
−1

0.05

−20
0

−1.2

0

0.05

0.1

0.15

0.2

0.25

0.3

0

50

100

150

0

Figure 6: Visualization of the Iris data set.
1046

0.5

1

1.5

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Data name∗ banana

breast-cancer
diabetes
flare-solar
germanheart
image
ringnorm
splice∗ thyroidtitanic
twonorm∗ waveform∗ USPS-eo∗ USPS-sl

Input
dimensionality
2
9
8
9
20
13
18
20
60
5
3
20
21
256
256

# of
training
samples
400
200
468
666
700
170
1300
400
1000
140
150
400
400
1000
1000

# of
test
samples
4900
77
300
400
300
100
1010
7000
2175
75
2051
7000
4600
1000
1000

# of
realizations
100
100
100
100
100
100
20
100
20
100
100
100
100
20
20

Table 2: List of binary classification data sets. Data sets indicated by ‘∗’ contain intrinsic withinclass multimodal structures.

implies that the heuristic to use the LFDA result as an initial value is useful. However, NCA does
not provide significant performance improvement over LFDA in the above simulations. The MCML
results have similar tendencies to FDA.
Based on the above simulation results, we conclude that LFDA is a promising method in thevisualization of multimodal labeled data.
5.3 Classification
Here we apply the proposed and existing dimensionality reduction techniques to classification tasks,
and objectively evaluate the effectiveness of LFDA.
There are several measures for quantitatively evaluating separability of data samples in different
classes (e.g., Fukunaga, 1990; Globerson et al., 2005). Here we use a simple one: misclassification rate by a one-nearest-neighbor classifier. As explained in Section 3.3, the LFDA criterion
is invariant under linear transformations, while the misclassification rate by a one-nearest-neighborclassifier depends on the distance metric. This means that the following simulation results are highly
dependent on the normalization scheme (15).
We employ the IDA data sets,5 which are standard binary classification data sets originally usedin Rätsch et al. (2001). In addition, we use two binary classification data sets created from the USPS
handwritten digit data set. The first task (USPS-eo) is to separate even numbers from odd numbers
and the second task (USPS-sl) is to separate small numbers (‘0’ to ‘4’) from large numbers (‘5’ to‘9’). For training and testing, 100 samples are randomly chosen for each digit. Table 2 summarizes
5. Data sets available at http://ida.first.fraunhofer.de/projects/bench/benchmarks.htm.

1047

S UGIYAMA

Data set
LFDA
LDI
NCA
MCML
LPP
PCA
∗ banana◦ 13.7 ± 0.8 ◦ 13.6 ± 0.8 14.3 ± 2.0◦◦
39.4 ± 6.7 13.6 ± 0.8 13.6 ± 0.8
breast-cancer ◦ 34.7 ± 4.3 36.4 ± 4.9 34.9 ± 5.0 ◦ 34.0 ± 5.8 ◦ 33.5 ± 5.4 ◦ 34.5 ± 5.0◦ 31.2 ± 2.1 31.5 ± 2.5 ◦ 31.2 ± 3.0
diabetes
32.0 ± 2.5 ◦ 30.8 ± 1.9—◦ 39.2 ± 4.9 ◦ 39.1 ± 5.1
flare-solar ◦ 39.2 ± 5.0 ◦ 39.3 ± 4.8——◦◦
german29.9 ± 2.8 30.7 ± 2.4 29.8 ± 2.6 31.3 ± 2.4 30.7 ± 2.4 ◦ 30.2 ± 2.4◦ 21.9 ± 3.7 23.9 ± 3.1
heart
23.0 ± 4.3 23.3 ± 3.8 23.3 ± 3.8 24.3 ± 3.5◦ 3.2 ± 0.8◦ 3.0 ± 0.6◦ 3.4 ± 0.5
image—
4.7 ± 0.8
3.6 ± 0.7◦
ringnorm
21.1 ± 1.3 17.5 ± 1.0 21.8 ± 1.3 22.0 ± 1.2 20.6 ± 1.1 21.6 ± 1.4◦◦ 17.3 ± 0.9 23.2 ± 1.2
16.9 ± 0.9 17.9 ± 0.8
splice—22.6 ± 1.3∗ thyroid◦ 4.6 ± 2.6◦
4.5 ± 2.2 18.5 ± 3.8 ◦ 4.2 ± 2.9 ◦ 4.9 ± 2.6
8.0 ± 2.9
◦ 33.1 ± 11.9 ◦ 33.1 ± 11.9 ◦ 33.0 ± 11.9 ◦ 33.1 ± 11.9 ◦ 33.0 ± 11.9 ◦ 33.0 ± 12.0
titanic◦ 3.5 ± 0.4◦ 3.5 ± 0.4
twonorm
4.1 ± 0.6
3.7 ± 0.6
3.7 ± 0.7
3.6 ± 0.6
∗ waveform◦ 12.5 ± 1.0 20.7 ± 2.5 ◦ 12.6 ± 0.8 17.9 ± 1.5 ◦ 12.4 ± 1.0 12.7 ± 1.2∗ USPS1◦ 3.5 ± 0.7
9.0 ± 0.8
12.5 ± 0.9——
7.2 ± 1.0
∗ USPS2◦ 3.9 ± 0.8
12.9 ± 1.2 25.9 ± 1.7—
11.7 ± 1.3
7.5 ± 0.8
Computation
1.00
1.11
97.23
70.61
1.04
0.91
time (ratio)
Table 3: Means and standard deviations of the misclassification rate when the embedding dimensionality is chosen by cross validation. For each data set, the best method and comparable
ones based on the t-test at the significance level 5% are marked by ‘◦’. Data sets indicatedby ‘∗’ contain the intrinsic within-class multimodal structure.

the specifications of the data sets. The ringnorm, twonorm, and waveform data sets contain featureswith only noise. The thyroid, waveform, USPS-eo, and USPS-sl data sets contain intrinsic withinclass multimodal structures since they are converted from multi-class problems by merging some of
the classes. The banana data set is also multimodal.
We test LFDA, LDI, NCA, MCML, LPP, and principal component analysis (PCA). Note that
LPP and PCA are unsupervised dimensionality reduction methods, while others are supervised
methods. NCA is not tested for the diabetes, flare-solar, image, splice, USPS-eo, and USPS-sl
data sets and MCML is not tested for the flare-solar and USPS-eo data sets since the execution timeis too long.
Figure 7 depicts the mean misclassification rate by a one-nearest-neighbor classifier as functionsof the dimensionality r of the reduced space. The error bars are omitted for clear visibility. Instead,
we plotted the results of the following significance test: for each dimensionality r, the mean misclassification rate by the best method and comparable ones based on the t-test (Henkel, 1979) at thesignificance level 5% are marked by ‘◦’. The results show that LFDA works quite well, but overall
there is no single best method that consistently outperforms the others.
Table 3 describes the mean and standard deviation of the misclassification rate by each method
when the embedding dimensionality r is chosen by 5-fold cross validation (Stone, 1974; Wahba,
1990); for the USPS-eo and USPS-sl data sets, we used 20-fold cross validation since this wasmore accurate. For each data set, the best method and comparable ones based on the t-test at
the significance level 5% are indicated by ‘◦’. The table shows that overall LFDA has excellent
1048

L OCAL F ISHER D ISCRIMINANT A NALYSIS

breast−cancer

bananadiabetes

0.38

0.4
0.38

0.3

0.25
LFDA
LDI
NCA
MCML
LPP
PCA

0.2

0.15

Mean Misclassification Rate

Mean Misclassification Rate

Mean Misclassification Rate

0.37

0.35

0.36
0.35
0.34
0.33

1

2

2

3

0.34
0.33
0.32

4
5
6
Reduced Dimension r

7

8

9

1

0.27

0.4

0.26

0.393
0.392
0.391

Mean Misclassification Rate

Mean Misclassification Rate

0.394

0.38
0.36
0.34
0.32
0.3

0.39
4
5
6
Reduced Dimension r

3

7

8

4
5
6
Reduced Dimension r

7

8

heart

0.42
0.395

3

2

germanflare−solar

Mean Misclassification Rate

0.35

0.3
1

Reduced Dimension r

2

0.36

0.31
0.32

1

0.37

0.25
0.24
0.23
0.22
0.21

5

9

image

10
15
Reduced Dimension r

20

2

4

ringnorm

6
8
10
Reduced Dimension r

12splice

0.35
0.45

0.36

0.25
0.2
0.15
0.1

Mean Misclassification Rate

0.34
Mean Misclassification Rate

Mean Misclassification Rate

0.3

0.32
0.3
0.28
0.26
0.24
0.22
0.2

0.05
4

6

8
10
12
Reduced Dimension r

14

16

0.2

18

5

thyroid10
15
Reduced Dimension r

20

0

10

20
30
40
Reduced Dimension rtitanic

50

60

twonorm

0.345

0.2

0.065

0.16
0.14
0.12
0.1
0.08

Mean Misclassification Rate

0.18
Mean Misclassification Rate

Mean Misclassification Rate

0.3
0.25

0.18
2

0.34

0.335

0.06
0.33

0.04
1

2

3
Reduced Dimension r

4

0.06
0.055
0.05
0.045
0.04
0.035

5

1

2
Reduced Dimension r

waveform

35

USPS−eo

Mean Misclassification Rate

0.24
0.22
0.2
0.18
0.16

0.5

0.45

0.45

0.4
0.35
0.3
0.25
0.2
0.15

0.4

0.3
0.25
0.2
0.15

0.1

0.1

0.12

0.05

0.05

5

10
15
Reduced Dimension r

20

0

50

100
150
Reduced Dimension r

20

0.35

0.14

0

10
15
Reduced Dimension r

USPS−sl

0.5

Mean Misclassification Rate

0.26

Mean Misclassification Rate

0.4
0.35

200

250

0

50

100
150
Reduced Dimension r

200

250

Figure 7: Mean misclassification rates by a one-nearest-neighbor method as functions of the dimensionality of the embedding space. For each dimension, the best method and comparable
ones based on the t-test at the significance level 5% are marked by ‘◦’.
1049

S UGIYAMA

Data set∗ banana
breast-cancer
diabetes
flare-solar
germanheart
image
ringnorm
splice∗ thyroidtitanic
twonorm∗ waveform∗ USPS-eo∗ USPS-sl

LFDA

EUCLIDFDA◦ 13.7 ± 0.8◦ 13.6 ± 0.8◦ 13.6 ± 0.8

34.7 ± 4.3 ◦ 32.7 ± 4.8 ◦ 32.9 ± 4.5
32.0 ± 2.5 ◦ 30.1 ± 2.1 ◦ 30.6 ± 2.2◦ 39.2 ± 5.0 ◦ 39.2 ± 5.1 ◦ 39.0 ± 4.9◦ 29.9 ± 2.8 ◦ 29.5 ± 2.5
30.5 ± 2.8◦ 21.9 ± 3.7
23.2 ± 3.7 24.0 ± 3.7◦ 3.2 ± 0.8◦ 3.4 ± 0.5
6.5 ± 1.7
◦ 21.1 ± 1.3
35.0 ± 1.4 31.2 ± 1.6◦ 16.9 ± 0.9
28.9 ± 1.5 33.7 ± 1.5◦ 4.6 ± 2.6◦ 4.4 ± 2.2
5.3 ± 2.5
◦ 33.1 ± 11.9 ◦ 33.0 ± 11.9 ◦ 33.1 ± 11.9
◦ 3.5 ± 0.4
6.7 ± 0.7
5.0 ± 0.7
◦ 12.5 ± 1.0
15.8 ± 0.7 17.6 ± 1.4◦ 3.6 ± 0.7
9.0 ± 0.8
15.1 ± 2.4
12.9 ± 1.2 ◦ 3.8 ± 0.8 13.3 ± 2.9

Table 4: Means and standard deviations of the misclassification rate. The LFDA results are copiedfrom Table 3. ‘EUCLID’ denotes a naive one-nearest-neighbor classification without dimensionality reduction. ‘FDA’ denotes a naive one-nearest-neighbor classification after
the samples are projected onto a one-dimensional FDA subspace.

performance. LDI and MCML also work quite well, but they tend to perform rather poorly for themultimodal data sets specified by ‘∗’. NCA also works well, but it does not compare favorably with
LFDA. Note that NCA with random initialization was slightly worse; therefore our heuristic to use
the LFDA results for initialization would be reasonable. LPP and PCA perform well, despite thefact that they are unsupervised dimensionality reduction methods. In particular, PCA has excellent
performance for the USPS data sets since the projection onto the two-dimensional PCA subspace
already gives reasonably separate embedding (He and Niyogi, 2004).
The computation time of each method summed over 9 data sets for which NCA is tested is
described in the bottom of Table 3. For better comparison, we normalized the values by the computation time of LFDA. This shows that LFDA is much faster than NCA and MCML, 6 and iscomparable to LDI, LPP, and PCA.
The misclassification rates by a naive one-nearest-neighbor classification without dimensionality reduction (‘EUCLID’) are described in Table 4. The table shows that, on the whole, the performance of LFDA is comparable to EUCLID. This implies that the use of LFDA is advantageouswhen the dimensionality of the original data is very high since the computation time in the test phasecan be reduced. Table 4 also includes the misclassification rates by a naive one-nearest-neighborclassification after the samples are projected onto a one-dimensional FDA subspace, showing that
LFDA tends to outperform FDA.
6. In our implementation of MCML, we used a constant step size for the gradient descent. The computation time could
be improved if, for example, an Armijo like step size rule (Bertsekas, 1976) is employed.

1050

L OCAL F ISHER D ISCRIMINANT A NALYSIS

Based on the above simulation results, we conclude that the proposed LFDA is a promisingdimensionality reduction technique also in classification scenarios.

6. Conclusions
We discussed the problem of supervised dimensionality reduction. FDA (Fisher, 1936; Fukunaga,
1990; Duda et al., 2001) works well for this purpose, given that data samples in each class are unimodal Gaussian. However, samples in a class are often multimodal in practice, for example, when
multi-class classification problems are solved by a set of two-class ‘one-versus-rest’ problems. LPP(He and Niyogi, 2004) can work well in dimensionality reduction of multimodal data. However, it
is an unsupervised method and does not necessarily useful in supervised dimensionality reduction
scenarios. In this paper, we proposed a new method called LFDA, which effectively combines theideas of FDA and LPP. LFDA allows us to reduce the dimensionality of multimodal labeled dataappropriately by maximizing between-class separability and preserving the within-class local structure at the same time. The derivation of LFDA is based on a novel pairwise interpretation of FDA(see Section 3.1). The original FDA provides a meaningful result only when the dimensionality ofthe embedding space is smaller than the number of classes because of the rank deficiency of the
between-class scatter matrix. On the other hand, LFDA does not share this limitation and can be
employed for dimensionality reduction into any dimensional spaces (see Section 3.3). This is asignificant improvement over the original FDA.
As discussed in Section 3.3, the LFDA criterion is invariant under linear transformations. Thismeans that the range of the transformation matrix can be uniquely determined, but the distancemetric in the embedding space cannot be determined. In this paper, we determined the distancemetric in a heuristic manner. Although this normalization scheme is shown to be reasonable inexperiments, there is still room for further improvement. An important future direction is to develop
a computationally efficient method of determining the distance metric of the embedding space,for example, following the lines of Goldberger et al. (2005), Globerson and Roweis (2006), andWeinberger et al. (2006).
We showed in Section 3.4 that a non-linear variant of LFDA can be obtained by employing thekernel trick. FDA, LPP, and MCML can also be kernelized similarly (Baudat and Anouar, 2000;
Mika et al., 2003; Belkin and Niyogi, 2003; He and Niyogi, 2004; Globerson and Roweis, 2006).
As shown in these papers, the performance of the kernelized methods heavily depend on the choiceof the family and parameters of kernel functions. Therefore, how to optimally determine the kernel
function for supervised dimensionality reduction needs to be explored.
The performance of LFDA depends on the choice of the affinity matrix. In this paper, we simply
employed a standard definition as it is (see Appendix D.4). Although this standard choice appeared
to be reasonable in experiments, it is important to find the optimal way to define the affinity matrix
in the context of supervised dimensionality reduction.
MDA (Hastie and Tibshirani, 1996b) provides a solid probabilistic framework for supervised
dimensionality reduction with multimodality (see Section 4.2). On the other hand, LFDA still lacks
a probabilistic interpretation. An interesting future direction is to analyze the behavior of LFDA interms of density models.

1051

S UGIYAMA

Acknowledgments
The author would like to thank Klaus-Robert Müller, Hideki Asoh, Motoaki Kawanabe, and Stefan
Harmeling for fruitful discussions. His gratitude also goes to anonymous reviewers for their valuable comments; particularly, the weighting scheme of the eigenvectors in Eq. (15) was suggested by
one of the reviewers. He also acknowledges financial support from MEXT (Grant-in-Aid for Young
Scientists 17700142 and Grant-in-Aid for Scientific Research (B) 18300057). A part of this work
has been done when the author was staying at University of Edinburgh, which is supported by EU
Erasmus Mundus Scholarship. He would like to thank Sethu Vijayakumar for the warm hospitality.

Appendix A. Proof of Lemma 1
It follows from Eq. (1) that
S

(w)

c

`=1 i:yi =`

1
xi −
xj
n` j:y∑
=`
j

n

c

=∑ ∑

!

1
xi −
xj
n` j:y∑=`j!>

1xi x>j∑n`=1 ` i, j:yi =y j =`!

= ∑ xi x>
i −∑
i=1
n

=∑

i=1

n

n

>
xi x>
i − ∑ Wi, j xi x j∑ Wi, j

j=1

(w)

(w)

i, j=1

1 n
(w)
>
>
>
= ∑ Wi, j (xi x>
i + x j x j − xi x j − x j xi ),2 i, j=1
which yields Eq. (5). Let S(m) be the mixture scatter matrix (Fukunaga, 1990):
S(m) ≡ S(w) + S(b)n= ∑ (xi − µ)(xi − µ)> .
i=1

Then we have
1 n
∑ xi x>j − S(w)ni, j=1
i=1
!
n
n
n
1
1 >
=∑ ∑
xi x>
xi x j − S(w)
i − ∑
n
n
i=1 j=1i, j=1


1
1 n
(w)
>
>
>
= ∑
−Wi, j (xi x>
i + x j x j − xi x j − x j xi ),2 i, j=1 nnS(b) = ∑ xi x>i −which yields Eq. (6).

Appendix B. Interpretation of FDA
In Section 3.1, we claimed that FDA tries to keep data pairs in the same class ‘close’ and data pairs
in different classes ‘apart’. Here we show this claim more formally.
1052

L OCAL F ISHER D ISCRIMINANT A NALYSIS

For
vi, j ≡ T > (xi − x j ),let us investigate the change in the Fisher criterion (3) when v i, j yields αvi, j with α > 0. Note that
there does not generally exist a transformation T 0 that keeps all vi, j and only changes a particularpair. For this reason, the following analysis may be regarded as comparing the values of the Fisher
criterion for two different data sets. This analysis will give an insight into what kind of transformation matrices the Fisher criterion favors.
Let
W ≡ T > S(w) T ,B ≡ T > S(b) T ,(w)

W α ≡ W − βWi, j vi, j v>
i, j ,(b)Bα ≡ B − βWi, j vi, j v>
i, j ,β≡1 − α2.
2Note that W α and Bα correspond to the within-class and between-class scatter matrices for αv i, j ,respectively. We assume that W and W α are positive definite and B and Bα are positive semi
definite. Thenthe values of the Fisher criterion (3) for vi, j and αvi, j are expressed as tr W −1 B and
tr W −1α Bα , respectively.
The standard matrix inversion lemma (e.g., Albert, 1972) yields(w)−1
W −1 = (W α + βWi, j vi, j v>
i, j )= W −1
α −
(w)

−1
>
W −1
α vi, j (W α vi, j )(w)(βWi, j )−1 + hW −1
α vi, j , vi, j i.

(b)

If yi = y j , we have Wi, j > 0 and Wi, j < 0. Then we have



(w)
(b)
−1
>
tr W −1 B = tr (W α + βWi, j vi, j v>
)
(B
+
βW
v
v
)
α
i, j
i, j i, j i, j(b)−1= tr W −1
α Bα + βWi, j hW α vi, j , vi, j i
(b)

−

= tr

−1
−1
2
hBαW −1
α vi, j ,W α vi, j i + βWi, j hW α vi, j , vi, j i(w)(βWi, j )−1 + hW −1
α vi, j , vi, j i
W −1
α Bα



(b)

−

(w) −1−1hBαW −1
α vi, j ,W α vi, j i −Wi, j Wi, jhW −1
α vi, j , vi, j i(w)(βWi, j )−1 + hW −1
α vi, j , vi, j i

If α < 1, we have β > 0 since α > 0 by definition. Therefore, Eq. (30) yields

−1
tr W −1B ,α Bα > tr W

.

(30)

where we used the facts that W α is positive definite and Bα is positive semi-definite. This implies
that the value of the Fisher criterion increases if a data pair in the same class is made close.
1053

S UGIYAMA

(w)

(b)

Similarly, if yi 6= y j , we have Wi, j = 0 and Wi, j > 0. Then we have



(w)
(b)
> −1
tr W −1
B
=
tr
(W
−
βW
v
v
)
(B
−
βW
v
v
)
α
α
i, j i, j i, j
i, j i, j i, j(b)= tr W −1 B − βWi, j hW −1 vi, j , vi, j i.

(31)

If α > 1, we have β < 0 and hence Eq. (31) yields


−1
B .
tr W −1α Bα > tr W

This implies that the value of the Fisher criterion increases if a pair of samples in different classes
are separated from each other.

Appendix C. Efficient Computation of T LFDA
As shown in Eq. (15), the LFDA transformation matrix T LFDA can be computed analytically using
the generalized eigenvectors and generalized eigenvalues of the following generalized eigenvalue
problem.
(b)(w)
ee = eλe
e.
S ϕ
S ϕ
(b)

(w)

, the computational complexity of calculating T LFDA is O (rd 2 ). Here, we provide(b)(w)
an efficient computing method of e
S and eS .

Given e
S

and e
S

(m)

Let e
S

be the local mixture scatter matrix defined by

(m)
(b)
(w)
e
S ≡e
S +e
S .

(m)
From Eqs. (9)–(11), we can immediately show that e
S is expressed in the following pairwise form.
newhere W

Since

(m)

(m)
1
e
e (m) (xi − x j )(xi − x j )> ,S = ∑W
2 i, j=1 i, j

is the n-dimensional matrix with (i, j)-th element being(
Ai, j /n if yi = y j ,e (m) ≡
Wi, j
1/n if yi 6= y j .
n

(m)
1
>
>
>
e
e (m) (xi x>
S = ∑Wi + x j x j − xi x j − x j xi )2 i, j=1 i, j
!
n
n
n
(m)
e (m) xi x>= ∑ ∑W−∑ We xi x>j ,
i
i, j

i=1

i, j

j=1

i, j=1(m)e
S can be expressed in a matrix form as(m)ee(m) X > ,S = XL1054(32)

L OCAL F ISHER D ISCRIMINANT A NALYSIS

where
e
and D

(m)

e(m) ≡ De (m) − We (m) ,
L(33)

is the n-dimensional diagonal matrix with i-th diagonal element beingne (m) .
e (m) ≡ ∑ W
D
i, j
i,i
j=1

(w)

Similarly, eS

can be expressed in a matrix form as(w)ee(w) X > ,S = XL

where
e
and D

(w)

(34)

e(w) ≡ De (w) − We (w) ,
L(35)

is the n-dimensional diagonal matrix with i-th diagonal element beingne (w) .
e (w) ≡ ∑ W
D
i, j
i,i
j=1

e(m) and Le(w) are n-dimensional matrices and could be very high dimensional. However, Le(w)L
can be made block-diagonal if the samples {xi }ni=1 are sorted according to the labels {yi }ni=1 . Fure(w) can be sparse if the affinity matrix A is sparsely defined
thermore, diagonal sub-matrices of L(w)(see Appendix D for detail). Therefore, directly calculating e
S by Eq. (34) may be already computationally efficient.
(m)e (m) is dense.
On the other hand, computing e
S directly by Eq. (32) is not so efficient since WeThis problem can be alleviated as follows. W(m)can be decomposed ase (m) = 1 1n 1>
e (m1 ) + W
e (m2 ) ,
W
n +W
n
(m1 )

e
e
where 1n is the n-dimensional vector with all ones and W
and W
matrices with (i, j)-th element being(
Ai, j /n if yi = y j ,
(m
)
1
e
W
≡
i, j
0
if yi 6= y j ,(−1/n if yi = y j ,e (m2 ) ≡Wi, j0if yi 6= y j .
(m)Then eS(m2 )are the n-dimensionalcan be expressed as(m)ee (m1 ) X > − X We (m2 ) X > ,e (m) X > − 1 X1n (X1n )> − X W
S = XDnewhere the diagonal matrix D(m)e
is expressed in terms of W
n

(m1 )

as

e (m1 ) .
e (m) = 1 − nyi + ∑ W
D
i, j
i,i
n
j=1
1055

(36)

S UGIYAMA

Note that nyi in the above equation is the number of samples in the class which the sample x i
e (m2 ) is a constant block-diagonal matrix if the samples {x i }n are sorted according
belongs to. W
i=1

e (m2 ) X > in the right-hand side of Eq. (36) can be computed
Therefore, X W
e (m1 ) can also be made block-diagonal, so X W
e (m1 ) X > in the right-hand side
efficiently. Similarly, W
of Eq. (36) may also be computed efficiently; if the affinity matrix A is sparse, the computational
efficiency can be further improved. The first two terms in the right-hand side of Eq. (36) can also(m)
be computed efficiently. Therefore, computing e
S by Eq. (36) may be more efficient than directly
to the labels {yi }ni=1 .

(b)

by Eq. (32). Finally, we can compute eS(m)efficiently by using e
S

as

(b)
(m)
(w)
e
S =e
S −e
S .

To further improve computational efficiency, the affinity matrix A may be computed in a classwise manner since we do not need the affinity values for sample pairs in different classes. This
speeds up the nearest neighbor search which is often carried out when defining A (see Appendix D).
The nearest neighbor search itself could also be a bottleneck, but this may be eased by incorporating the prior knowledge of the data structure or by approximation (see Saul and Roweis, 2003, andreferences therein).
The above efficient implementation of LFDA is summarized as a pseudo code in Figure 2.

Appendix D. Definitions of Affinity Matrix
Here, we briefly review typical choices of the affinity matrix A.
D.1 Heat Kernel
A standard choice of the affinity matrix A iskxi − x j k2
Ai, j = exp −
σ2



,

(37)

where σ (> 0) is a tuning parameter which controls the ‘decay’ of the affinity (e.g., Belkin andNiyogi, 2003).
D.2 Euclidean Neighbor
The heat kernel gives a non-sparse affinity matrix. It would be computationally advantageous if theaffinity matrix is made sparse. A sparse affinity matrix can be obtained by assigning positive affinity
values only to neighboring samples. More specifically, x i and x j are said to be neighbors if
kxi − x j k ≤ ε,where ε (> 0) is a tuning parameter. Then Ai, j is defined by Eq. (37) for two neighboring samples
and Ai, j = 0 for non-neighbors (Tenenbaum et al., 2000).
This definition includes two tuning parameters (ε and σ), which are rather troublesome to determine in practice. To ease the problem, we may simply let A i, j = 1 if xi and x j are neighbors andAi, j = 0 otherwise. This corresponds to setting σ = ∞.
1056

L OCAL F ISHER D ISCRIMINANT A NALYSIS

D.3 Nearest Neighbor
Tuning the distance threshold ε is practically rather cumbersome since the relation between thenumber of neighbors and the value of ε is not intuitively clear. Another option to determine theneighbors is to directly specify the number of neighbors (Roweis and Saul, 2000; Tenenbaum et al.,(K)2000). Let NNi be the set of K nearest neighbor samples of xi under the Euclidean distance, where(K)(K)
K is a tuning parameter. If x j ∈ NNi or xi ∈ NN j , xi and x j are regarded as neighbors; otherwise
they are regarded as non-neighbors. Then the affinity matrix is defined by the heat kernel or in thesimple zero-one manner.

D.4 Local Scaling
A drawback of the above definitions could be that the affinity is computed globally in the same way.
The density of data samples may be different depending on regions. Therefore, it would be more
appropriate to take the local scaling of the data into account. Following this idea, Zelnik-Manor andPerona (2005) proposed defining the affinity matrix askxi − x j k2
Ai, j = exp −σi σ j.

σi represents the local scaling of the data samples around x i , which is determined by(K)

σi = kxi − xi k,(K)where xi is the K-th nearest neighbor of xi . The parameter K is a tuning parameter, but ZelnikManor and Perona (2005) demonstrated that K = 7 works well on the whole. This would be aconvenient heuristic for those who do not have any subjective/prior preferences. We employed the
local scaling method with this heuristic all through the paper.
For computational efficiency, we may further sparsify the above affinity matrix based on, forexample, the nearest neighbor idea, although this is not tested in this paper.

Appendix E. Pairwise Expression of S
A pairwise expression of S

(b)

(b)

can be derived as1057

S UGIYAMA

S

(b)

n

1 c [k] [k][k]n (µ` − µ[k] )(µ` − µ[k] )>[k] ∑ `
n
k=1
`=1

=∑
n

1
= ∑ [k]k=1 nn1= ∑ [k]
k=1 n
=

1 n 1
∑ n[k]
2 k=1

c

∑

`=1

>
[k] [k] [k] >n` µ` µ` − n[k] µ[k] µ[k]!

nn[k] [k] [k] >[k] [k] [k] >n` µ` µ` − Ai,k xi x>+
Ai,k xi x>ii −n µ µ
i=1
i=1
`=1
c

∑

∑

c

−∑

1

∑

∑

Ai,k A j,k (xi − x j )(xi − x j )[k]`=1 n` i, j:yi =y j =`>1 n
+ [k] ∑ Ai,k A j,k (xi − x j )(xi − x j )>n i, j=1
=

!

!

1 n
(b)
W i, j (xi − x j )(xi − x j )> ,
2 i,∑
j=1

which yields Eqs. (19) and (20).

References
A. Albert. Regression and the Moore-Penrose Pseudoinverse. Academic Press, New York andLondon, 1972.
N. Aronszajn. Theory of reproducing kernels. Transactions of the American Mathematical Society,68:337–404, 1950.
G. Baudat and F. Anouar. Generalized discriminant analysis using a kernel approach. NeuralComputation, 12(10):2385–2404, 2000.
M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation.
Neural Computation, 15(6):1373–1396, 2003.
D. P. Bertsekas. On the Goldstein-Levitin-Polyak gradient projection method. IEEE Transactions
on Automatic Control, 21(2):174–184, 1976.
C. L. Blake and C. J. Merz. UCI repository of machine learning databases, 1998.
http://www.ics.uci.edu/˜mlearn/MLRepository.html.

URL

S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge,
2004.
F. R. K. Chung. Spectral Graph Theory. American Mathematical Society, Providence, R.I., 1997.
D. Coomans, M. Broeckaert, M. Jonckheer, and D. L. Massart. Comparison of multivariate discriminant techniques for clinical data—Application to the thyroid functional state. Methods ofInformation in Medicine, 22:93–101, 1983.
1058

L OCAL F ISHER D ISCRIMINANT A NALYSIS

A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the
EM algorithm. Journal of the Royal Statistical Society, series B, 39(1):1–38, 1977.
R. O. Duda, P. E. Hart, and D. G. Stor. Pattern Classification. Wiley, New York, 2001.
N. Duffy and M. Collins. Convolution kernels for natural language. In T. G. Dietterich, S. Becker,
and Z. Ghahramani, editors, Advances in Neural Information Processing Systems 14, Cambridge,
MA, 2002. MIT Press.
B. S. Everitt, S. Landau, and M. Leese. Cluster Analysis. Arnold, London, 2001.
R. A. Fisher. The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7(2):179–188, 1936.
J. H. Friedman. Regularized discriminant analysis. Journal of the American Statistical Association,
84:165–175, 1989.
K. Fukunaga. Introduction to Statistical Pattern Recognition. Academic Press, Inc., Boston, second
edition, 1990.
T. Gärtner. A survey of kernels for structured data. SIGKDD Explorations, 5(1):S268–S275, 2003.
T. Gärtner, P. Flach, and S. Wrobel. On graph kernels: Hardness results and efficient alternatives.
In Proceedings of the Sixteenth Annual Conference on Computational Learning Theory, 2003.
A. Globerson, G. Chechik, F. Pereira, and N. Tishby. Euclidean embedding of co-occurrence data.
In L. K. Saul, Y. Weiss, and L. Bottou, editors, Advances in Neural Information Processing
Systems 17, pages 497–504. MIT Press, Cambridge, MA, 2005.
A. Globerson and S. Roweis. Metric learning by collapsing classes. In Y. Weiss, B. Sch ölkopf,
and J. Platt, editors, Advances in Neural Information Processing Systems 18, pages 451–458,
Cambridge, MA, 2006. MIT Press.
J. Goldberger, S. Roweis, G. Hinton, and R. Salakhutdinov. Neighbourhood components analysis.
In L. K. Saul, Y. Weiss, and L. Bottou, editors, Advances in Neural Information Processing
Systems 17, pages 513–520. MIT Press, Cambridge, MA, 2005.
J. Ham, D. D. Lee, S. Mika, and B. Schölkopf. A kernel view of the dimensionality reduction of
manifolds. In Proceedings of the Twenty-First International Conference on Machine Learning,
New York, NY, 2004. ACM Press.
T. Hastie and R. Tibshirani. Discriminant adaptive nearest neighbor classification. IEEE Transactions on Pattern Analysis and Machine Intelligence, 18(6):607–615, 1996a.
T. Hastie and R. Tibshirani. Discriminant analysis by Gaussian mixtures. Journal of the Royal
Statistical Society, Series B, 58(1):155–176, 1996b.
X. He and P. Niyogi. Locality preserving projections. In S. Thrun, L. Saul, and B. Sch ölkopf,
editors, Advances in Neural Information Processing Systems 16. MIT Press, Cambridge, MA,
2004.
1059

S UGIYAMA

R. E. Henkel. Tests of Significance. SAGE Publication, Beverly Hills, 1979.
G. E. Hinton and R. R. Salakhutdinov. Reducing the dimensionality of data with neural networks.
Science, 313(5786):504–507, 2006.
H. Kashima and T. Koyanagi. Kernels for semi-structured date. In Proceedings of the Nineteenth International Conference on Machine Learning, pages 291–298, San Francisco, CA, 2002. Morgan
Kaufmann.
H. Kashima, K. Tsuda, and A. Inokuchi. Marginalized kernels between labeled graphs. In Proceedings of the Twentieth International Conference on Machine Learning, San Francisco, CA, 2003.
Morgan Kaufmann.
T. Kohonen. Self-Organization and Associative Memory. Springer-Verlag, Berlin, 1989.
R. I. Kondor and J. Lafferty. Diffusion kernels on graphs and other discrete input spaces. In
Proceedings of the Nineteenth International Conference on Machine Learning, pages 315–322,
2002.
S. Kullback and R. A. Leibler. On information and sufficiency. Annals of Mathematical Statistics,
22:79–86, 1951.
H. Lodhi, C. Saunders, J. Shawe-Taylor, N. Cristianini, and C. Watkins. Text classification using
string kernels. Journal of Machine Learning Research, 2:419–444, 2002.
J. B. MacQueen. Some methods for classification and analysis of multivariate observations. In Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability, volume 1,
pages 281–297, Berkeley, CA., USA, 1967. University of California Press.
S. Mika, G. Rätsch, J. Weston, B. Schölkopf, A. Smola, and K.-R. Müller. Constructing descriptive and discriminative nonlinear features: Rayleigh coefficients in kernel feature spaces. IEEE
Transactions on Pattern Analysis and Machine Intelligence, 25(5):623–628, 2003.
G. Rätsch, T. Onoda, and K.-R. Müller. Soft margins for adaboost. Machine Learning, 42(3):287–320, 2001.
S. Roweis and L. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science,
290(5500):2323–2326, 2000.
L. K. Saul and S. T. Roweis. Think globally, fit locally: Unsupervised learning of low dimensional
manifolds. Journal of Machine Learning Research, 4(Jun):119–155, 2003.
B. Schölkopf, A. Smola, and K.-R. Müller. Nonlinear component analysis as a kernel eigenvalue
problem. Neural Computation, 10(5):1299–1319, 1998.
B. Schölkopf and A. J. Smola. Learning with Kernels. MIT Press, Cambridge, MA, 2002.
M. Stone. Cross-validatory choice and assessment of statistical predictions. Journal of the Royal
Statistical Society, Series B, 36:111–147, 1974.
1060

L OCAL F ISHER D ISCRIMINANT A NALYSIS

J. B. Tenenbaum, V. de Silva, and J. C. Langford. A global geometric framework for nonlinear
dimensionality reduction. Science, 290(5500):2319–2323, 2000.
V. N. Vapnik. Statistical Learning Theory. Wiley, New York, 1998.
G. Wahba. Spline Model for Observational Data. Society for Industrial and Applied Mathematics,
Philadelphia and Pennsylvania, 1990.
K. Weinberger, J. Blitzer, and L. Saul. Distance metric learning for large margin nearest neighborclassification. In Y. Weiss, B. Schölkopf, and J. Platt, editors, Advances in Neural Information
Processing Systems 18, pages 1473–1480, Cambridge, MA, 2006. MIT Press.
L. Zelnik-Manor and P. Perona. Self-tuning spectral clustering. In L. K. Saul, Y. Weiss, and
L. Bottou, editors, Advances in Neural Information Processing Systems 17, pages 1601–1608.
MIT Press, Cambridge, MA, 2005.

1061