# Trace Ratio Criterion for Feature Extraction in Classification

Authors: Guoqi Li; Changyun Wen; Wei Wei; Yi Xu; Jie Ding; Guangshe Zhao; Luping Shi
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A generalized linear discriminant analysis based on trace ratio criterion algorithm (GLDA-TRA) is derived to extract features for classification. With the proposed GLDA-TRA, a set of orthogonal features can be extracted in succession. Each newly extracted feature is the optimal feature that maximizes the trace ratio criterion function in the subspace orthogonal to the space spanned by the previous extracted features.


---

## Body

## 1. Introduction

Linear discriminant analysis (LDA) [1–3] has been proposed as a class separatory measure, which has been intensively used to reduce dimensionality of a classification problem as well as improve the generalization capability of a pattern classifier. Generally speaking, LDA method is to optimize the ratio criterion of the between-class distance and within-class distance constructed based on the available learning data. Such optimization can be realized by solving a generalized eigenvalue problem of the between-class and within-class scatter matrices [4].In our opinion, there are three main problems for LDA methods. The first problem is its difficulty of dealing with high-dimensional data, where the number of observed samples is much lower than the samples’ feature dimension [5]. Many methods have been studied and proposed to address this problem; see, for example, the linear discriminant feature selection (LDFS) [6], Sparse Discriminant Analysis [5, 7], and Sparse Tensor Discriminant Analysis [8].The second problem is the well-known undersampled problem [9] in LDA method, in which scatter matrices may become singular due to insufficient samples. The solutions to this problem have also been well investigated in the following works such as the Regularized LDA [10, 11] using regularization techniques [12, 13] and the Penalized LDA [14], the Pseudo Fisher Linear Discriminant [15], the Generalized Singular Value Decomposition [16], and the Uncorrelated LDA [17] and the Orthogonal LDA [17].Basically the above two problems are quite similar and they can be unified as the same problem, which has also been extensively investigated in the above schemes. However, the third problem due to the LDA method can only extract quite limited features for classification problems [4]. For example, in two-class classification, one can only find one nonzero eigenvalue (extracted feature), as the between-class scatter matrix is a rank-one matrix. To the best of our knowledge, currently there is no good way to deal with this problem yet.We focus on the third problem in this paper. A generalized LDA based on trace ratio criterion [18–23] is proposed to overcome such a problem and an algorithm called generalized linear discriminant analysis based on trace ratio criterion algorithm (GLDA-TRA) is derived to extract features from the input feature space. The algorithm first extracts a feature which maximizes the trace ratio criterion by solving a generalized eigenvalue problem. It is shown that such a generalized eigenvalue problem is the same as the generalized eigenvalue problem of LDA. Then, the learning data are projected to a subspace orthogonal to the space spanned by the extracted features. In that orthogonal subspace, the algorithm continues to extract a feature which maximizes the proposed trace ratio criterion. This process continues and, in this way, a set of orthogonal features is obtained iteratively. It is proven that each newly extracted feature is the optimal feature that maximizes the trace ratio criterion in the subspace orthogonal to the space spanned by the previous extracted features. Finally the extracted features are shown to give a sequence of trace ratios with magnitudes monotonically decreasing.


## 2. Problem Formulation

Let(x,y)∈Rd×Y be a sample, where Rd denotes a d-dimensional feature space and Y={1,2,…,C} is a label set. Let xij denote the ith sample in the jth class. The within-class scatter matrix SW and the between-class scatter matrix SB [24] are, respectively, defined as
(1)SW=1n∑j=1C∑i=1ni(xij-mj)(xij-mj)′,SB=1n∑j=1C∑i=1ni(mj-m)(mj-m)′,
where ni is the number of samples in the ith class and n is the total number of samples and mj is the sample mean of the jth class and m is the sample mean of all classes and the notation ′ means matrix transpose. Without loss of generality, we assume that n≫d≫C. Then SW can be constructed as a full rank matrix of rank d while SB is at most with rank C-1. So, in this paper, SW is considered as a symmetrical and positive definite matrix and SB is a symmetrical and nonnegative definite matrix. The LDA method extracts features, that is, the column vectors in matrix Wopt, in such a way that the ratio of the between-class scatter and the within-class scatter is maximized [25]. Consider
(2)Wopt=arg⁡max⁡W⁡|W′SBW||W′SWW|=[w1w2⋯wm],
where |·| is the determinant of a matrix and {wi∣i=1,…,m} is the set of generalized eigenvectors of SB and SW corresponding to the m largest generalized eigenvalues {λi∣i=1,…,m} such that
(3)SBwi=λiSWwi,i=1,…,m.
Unfortunately, there are at most C-1 nonzero generalized eigenvalues as the rank of SB is at most C-1. As such, for a two-class classification problem, LDA can only extract one feature.To overcome this problem, we hope that one can continue to extractw2 after w1 is extracted. Generally speaking, we hope to extract wi+1 after the features w1,…,wi are extracted with i starting from 0. We now use the trace ratio criterion function in [24, 26] to formulate the above-mentioned feature extraction problem. The optimization model we proposed is shown as follows:
(4)wi+1=arg⁡max⁡⁡wi+1⁡F(W)=12tr⁡(W′SBW)tr⁡(W′SWW)subjecttoW=Wi+1=[Wiwi+1]hhhhhhhhhwi+1isorthonormaltospan⁡{Wi},
where Wi=[w1w2⋯wi] is a matrix denoting all the features extracted with W0 being empty, wi+1 is the feature to be determined in the space orthogonal to span⁡{Wi}, and
(5)F(W)=12w1SBw1+⋯+wiSBw1+wi+1SBwi+1w1SWw1+⋯+w1SWw1+wi+1SWwi+1,
with tr⁡(·) denoting the trace of a matrix.Remark 1.
Later it can be shown that the first extracted vectorw1 which maximizes F(W) in the space Rd is the same w1 found by LDA in (2). Here it is also necessary to point out that the orthogonal constraint condition in W is needed in our formulated problem (4). As in (2), the numerator is the determinant of matrix W′SBW. Implicitly there is a constraint that w1,…,wi cannot be the same. Because when w1=⋯=wi, the numerator |W′SBW| would be zero. But in (4), the numerator is the trace of W′SBW, which does not exclude the possibility that w1,…,wi are the same. With the constraint in (4), such a possibility can be avoided.


## 3. Proposed Algorithm and Analysis

In this section, we present and analyze the proposed feature extraction algorithm. Our idea is summarized as follows. We first extract a featurew1 by maximizing the trace ratio criterion function involving SW and SB in (4). When the current extracted features become Wi=[w1,…,wi], let span⁡{Wi} denote a space spanned by the linear combination of all the columns of Wi and let span⁡{Wi}⊥ denote the space orthogonal to span⁡{Wi}. Then, SW and SB are projected onto the subspace span⁡{Wi}⊥ by using projection operators (I-WiWi+) and WiWi+, respectively, where Wi+=(Wi′Wi)-1Wi′ is the generalized matrix inverse of a column full rank matrix Wi. We continue the process to find wi+1 by optimizing (4) until all m(m≤d) features are extracted.We first present the algorithm in Section3.1. In Section 3.2., we will show how this algorithm is derived and then analyze its properties.


### 3.1. Proposed Algorithm

For convenience, we present the definition of a generalized eigenvalue as follows.Definition 2.
A numberλ is called a generalized eigenvalue of matrix B with respect to A if λ satisfies that Bx=λAx for a nonzero vector x, where A is a positive definite symmetrical matrix. When A=I, λ is a normal eigenvalue of B.Now the algorithm is given as follows.Initialization Step: construct symmetrical matricesSW1=SW∈Rd×d and SB1=SB∈Rd×d as shown in (1) based on the available learning data with W0 as an empty matrix.Stepi  (i=1,…,m) is as follows.(a)
Calculation stage: find a unit vector in the direction of a generalized eigenvector which corresponds to the maximum eigenvalue of the generalized eigenvalue problem ofSBi with respect to SWi. This can be achieved by the following process:do the Cholesky decompositionSWi=GiGi′. Let Si=Gi-1Bi-1(Gi-1)′ and obtain its maximum eigenvalue λi together with its corresponding eigenvector xi. Choose wi=(Gi-1)′xi/∥(Gi-1)′xi∥2 and Wi=[Wi-1wi].(b)
Update stage is as follows:(6)SWi+1=(I-WiWi+)SW(I-WiWi+)+μWiWi+SWWiWi+,SBi+1=(I-WiWi+)SB(I-WiWi+),where u is a sufficiently small positive number.Replacei by i+1 until m features are extracted in W.Remark 3.
(1) AsW is orthonormal, (6) can be rewritten as the following recursive from:
(7)SWi+1=(I-wiwi+)SWi(I-wiwi+)+μ∑j=1i‍wjwj+SWwjwj+,SBi+1=(I-wiwi+)SBi(I-wiwi+).
Later in Lemma 5, it will be shown that SWi for i=1,…,m can still be positive definite for a sufficiently small positive u. At each step i, SWi is positive definite and the generalized eigenvalue of SBi with respect to SWi exists.
(2) The reason why we need to find the maximum eigenvalue ofSi in the Calculation stage is shown in Lemma 8.
(3) Suppose thati features w1,…,wi,(i<d) have been extracted. Then Theorem 12 shows that wi+1 is found to ensure the trace ratio criterion function in (4) to attain its maximum value.
(4) Whenm=1, the LDA and GLDA-TRA extract the same feature.
(5) GLDA-TRA extractsm≤d features one by one. When i=m=d, SBi+1 becomes a zero matrix, and the algorithm will not extract any more features.


### 3.2. Derivation and Analysis of GLDA-TRA

Lemma 4.
For an arbitrary full column rank matrixWi∈Rd×i, WiWi+=Wi(Wi′Wi)-1Wi′ and I-WiWi+ are projection operators which project a vector onto span⁡{Wi} and span⁡{Wi}⊥, respectively.Proof.
Suppose thatv∈Rm×1 belongs to span⁡{Wi}; there exists a vector x∈Ri×1 such that v=Wix. It can be obtained that WiWi+v = Wi(Wi′Wi)-1Wi′Wix = Wix=v and (I-WiWi+)v = (I-Wi(Wi′Wi)-1Wi′)v = (I-Wi(Wi′Wi)-1Wi′)Wix=0. This lemma holds.Lemma 5.
S
W
i
+
1 in (6) for i=1,…,m are positive definite for a sufficiently small positive μ.Proof.
Note thatxSW1x>0 for any nonzero x∈Rm×1 as SW1=SW. Let x1=(I-WiWi+)x and x2=WiWi+x. Then x=(I-WiWi+)x + WiWi+x=x1+x2. As x≠0, x1 and x2 cannot be zero at the same time. Thus
(8)x′SWi+1x=x′(I-WiWi+)SW(I-WiWi+)x+μx′(WiWi+)SW(WiWi+)x=((I-WiWi+)x)′SW(I-WiWi+)x+μ(WiWi+x)′SWWiWi+x=x1′SWx1+μx2′SWx2>0,
for a sufficiently small positive number μ.Lemma 6.
For matricesW∈Rd×m, X∈Rd×d, define g(W)=tr⁡(W′XW). Then dg/dW=(X′+X)W.Proof.
LetW=[w1⋯wi⋯wm] and X=[x1⋯xd], where wi∈Rd×1 for i=1,…,m, and xi∈Rd×1 for i=1,…,d. We have
(9)g(W)=tr⁡(W′XW)=tr⁡(∑i=1dw1′xiwi1+⋯+∑i=1dwk′xiwik+⋯+∑i=1dwm′xiwi1)=∑i=1d∑j=1dwj1xjiwi1+⋯+∑i=1d∑j=1dwjkxjiwik+∑i=1d∑j=1dwjmxjiwim.
Then, for k0=1,…,m, j0=1,…,d, we get
(10)dg(W)dwj0k0=∑i=1dxj0iwik0+∑j=1dwjk0xjj0=∑i=1dxj0iwik0+∑i=1dxij0wik0.
So (10) gives dg/dW=(X′+X)W.Lemma 7.
Letw∈Rd×1 and the trace ratio function f(w)=(1/2)(tr⁡(w′SWw)/tr⁡(w′SBw)). The gradient ∇f(w)=df(w)/dw = (wSW·tr⁡(w′SBw)-wSB·tr⁡(w′SWw))/(tr⁡(w′SBw))2.Proof.
Letg1(w)=tr⁡(w′SWw) and g2(w)=tr⁡(w′SBw). From Lemma 6, ∇f(w)=df(w)/dw = ((dg1(w)/dw)g2(w)-(dg2(w)/dw)g1(w))/(g2(w))2 = (wSW·tr⁡(w′SBw)-wSB·tr⁡(w′SWw))/(tr⁡(w′SBw))2.DefineR(w)=tr⁡(w′SBiw)/tr⁡(w′SWiw), and we have the following results mentioned in the Calculation stage of step i.Lemma 8.
Assume that the maximum eigenvalue ofSi is λ* with an eigenvector x*. Then the unit vector w*=(Gi-1)′x*/∥(Gi-1)′x*∥ is an extracted feature ensuring that R(w) attains its maximum value which is equal to λ*.Proof.
Note thatf(w)=(1/2)(1/R(w))=(1/2)(tr⁡(w′SWiw)/tr⁡(w′SBiw)). Thus maximizing R(w) is equivalent to minimizing f(w). From Lemma 7, we have
(11)∇f(w(k))=df(w)dw|w=w(k)=wSWi·tr⁡(w′SBiw)-wSB·tr⁡(w′SWiw)(tr⁡(w′SBiw))2|w=w(k).
Then one can minimize f(w) by using an iterative method. Let w(k+1)-w(k)=-η∇f(w(k)), where k denotes the kth iteration and η is a small positive constant. Then f(w(k+1)) = f(w(k))-η(∇f(w(k)))′(∇f(w(k))) + o(∥w(k+1)-w(k)∥22)≤f(w(k)). Thus {f(w(k))} is a nonincreasing positive sequence and its limit, denoted as f*(w), exists. We have lim⁡k→∞f(w(k+1))=f(w(k))=f*(w). This implies that lim⁡k→∞η(∇f(wk))′(∇f(w(k)))=0. As η>0, then {w(k)} converges to an accumulation point w¯ that satisfies ∇f(w)|w=w¯=0. This gives
(12)df(w)dw|w=w¯=w¯SWi·tr⁡(w¯′SBiw¯)-w¯SBi·tr⁡(w¯′SWiw¯)[tr⁡(w¯′SBiw¯)]2=0.
Let λ=tr⁡(w¯′SBiw¯)/tr⁡(w¯′SWiw¯). From (12), we have SBiw¯=λSWiw¯. That is, λ is the generalized eigenvalue obtained from
(13)SBiw=λSWiw.
Note that each eigenvector w in (13) is an accumulation point of f(w), since it satisfies (12). Now we hope to convert the generalized eigenvalue problem of (13) to a normal eigenvalue problem. This is achieved by Cholesky decomposition of SWi, which is given as SWi=GiGi′, where Gi is a full column rank lower triangular matrix. By doing this, (13) becomes
(14)SBiw=λGiGi′w.
Defining x=Gi′w and substituting x into (14), it can be obtained that Six=λx, where Si=Gi-1SBi(Gi-1)′ is a symmetrical positive definite matrix. Let the maximum eigenvalue and its corresponding eigenvector of Si be λ* and x*. Finally the optimal w*=(Gi-1)′x*/∥(Gi-1)′x*∥ is a unit vector and λ*=max⁡⁡(R(w))=tr⁡(w*′SBw*)/tr⁡(w*′SWw*).Theorem 9.
Suppose thatwi is a feature extracted at step i; the trace ratio (1/2)(tr⁡(wi′SBi+1wi)/tr⁡(wi′SWi+1wi))=0 at step i+1.Proof.
Note thattr⁡(wi′SBi+1wi)=wi(I-WiWi+)SB(I-WiWi+)wi. We have tr⁡(wi′SWi+1wi)>0 since SWi+1 is positive definite by Lemma 5. Also, (I-WiWi+)wi=0 as (I-WiWi+)wi is a projection operator which projects a vector onto span⁡{Wi}⊥ based on Lemma 4. Thus, (1/2)(tr⁡(wi′SBi+1wi)/tr⁡(wi′SWi+1wi))=0.Remark 10.
As seen in Theorem9, if the feature wi has been extracted at step i, then wi is supposed to make the trace ratio (1/2)(tr⁡(wSBiw)/tr⁡(wSWiw)) attain its maximum in this step. While at step i+1, after SWi and SBi are updated to SWi+1 and SBi+1, respectively, we have (1/2)(tr⁡(wSBi+1w)/tr⁡(w′SWi+1w))=0 at w=wi. This means that the algorithm needs to find wi+1 which can maximize (1/2)(tr⁡(wSBi+1w)/tr⁡(w′SWi+1w)), and, obviously, wi+1 must be different from wi.Theorem 11.
Sequence{F(wi)}i=1m produced by GLDA-TRA is a decreasing sequence.Proof.
Assume thatVi is a space which the ith feature wi belongs to. That is, wi is the feature extracted from Vi such that F(wi) attains its maximum. Now we consider wi∈Vi at step i and wi+1∈Vi+1 at step i+1 with the respective corresponding maximum eigenvalues λi and λi+1. As SWi+1 is positive definite; from Lemma 4 it can be obtained that
(15)∥SWi+1w∥2=∥((I-WiWi+)SW(I-WiWi+)+μWiWi+SWWiWi+)w∥2=∥(I-WiWi+)SW(I-WiWi+)w∥2+μ∥WiWi+SWWiWi+w∥2⟶∥(I-WiWi+)SW(I-WiWi+)w∥2,
when μ→0. If w∈span⁡{Wi}, (I-WiWi+)w=0. Then, ∥SWi+1w∥2 = ∥(I-WiWi+)SW(I-WiWi+)w∥2=0. If w∈span⁡{Wi}⊥, (I-WiWi+)w=w and then ∥SWi+1w∥2>0. Also, as λi and λi+1 correspond to the maximum eigenvalues at steps i and i+1, it is obvious that
(16)∥λiSWiwi∥2≥μ∥wi∥2>0,∥λi+1SWi+1wi+1∥2≥μ∥wi+1∥2>0.
Thus, it can be concluded that wi∈span⁡{Wi} while wi+1∈Vi+1⊆span⁡{Wi}⊥. Similarly, wi∈Vi⊆span⁡{Wi-1}⊥. As Wi=[w1,…,wi], span⁡{Wi} increases as i increases. Thus, Vi decreases as i increases. In addition, we have V1=Rd=V2⊕span⁡{w1} = V3⊕span⁡{w1}⊕span⁡{w2} = ⋯ = Vm⊕span⁡{w1}⊕⋯⊕span⁡{wm-1}. Thus, F(w1)≥F(w2)≥⋯≥F(wm).Theorem 12.
The featurewi+1 extracted by GLDA-TRA is the optimal feature maximizing F(Wi+1) in (4) when Wi=[w1⋯wi] is extracted.Proof.
Note thatwi+1∈span⁡{Wi}⊥ as wi+1 is orthonormal to span⁡{Wi}. Consider an arbitrary w~i+1≠wi+1∈span⁡{Wi}⊥. Construct Wi+1=[Wiwi+1] and W~i+1=[Wiw~i+1]. Note that Wi is an orthonormal matrix. Then wi0′wj0=0 as long as i0≠j0 and we have
(17)F(Wi+1)-F(W~i+1)=wi+1SBwi+1+⋯+wi+1SBwi+1+wi+1SBwi+1wi+1SWwi+1+⋯+wi+1SWwi+1+wi+1SWwi+1-wi+1SBwi+1+⋯+wi+1SBwi+1+w~i+1SBw~i+1wi+1SWwi+1+⋯+wi+1SWwi+1+w~i+1SWw~i+1=(wi+1SBwi+1w~i+1SWw~i+1-w~i+1SBw~i+1wi+1SWwi+1)×((w1SBw1)2+⋯+(wiSBw1)2+(wi+1SBwi+1)(w~i+1SWw~i+1)((w1SBw1)2+⋯+(wiSBw1)2)-1.
Denote that SB=(I-WiWi+)SB+WiWi+SB and SW = (I-WiWi+)SW + WiWi+SW. Then, WiWi+wi+1=0 and (I-WiWi+)wi+1=wi+1. Thus, we have
(18)wi+1SBwi+1=wi+1((I-WiWi+)SB+WiWi+SB)wi+1=wi+1((I-WiWi+)SB(I-WiWi+))wi+1=wi+1SBi+1wi+1,
where SBi+1=((I-WiWi+)SB(I-WiWi+)). Similarly, it can be obtained that
(19)wi+1SWwi+1=wi+1SWi+1wi+1,w~i+1SBw~i+1=w~i+1SBi+1w~i+1,w~i+1SWw~i+1=w~i+1SWi+1w~i+1,
where SWi+1=((I-WiWi+)SW(I-WiWi+)). From Lemma 8, it is noted that tr⁡(wSBi+1w)/tr⁡(wSWi+1w) attains its maximum if and only if w is parallel to the eigenvector corresponding to the maximum generalized eigenvalue of SBi+1 with respect to SWi+1. Since wi+1 is parallel to the generalized eigenvector, as long as wi+1≠w~i+1, we have
(20)tr⁡(wi+1SBi+1wi+1)tr⁡(wi+1SWi+1wi+1)>tr⁡(w~i+1SBi+1w~i+1)tr⁡(w~i+1SWi+1w~i+1).
That is to say, wi+1SBwi+1w~i+1SWw~i+1>w~i+1SBw~i+1wi+1SWwi+1. So we have F(Wi+1)>F(W~i+1). Thus, this theorem holds.Remark 13.
Our main theoretical conclusions of this paper are shown in Theorems9–12. The implication of Theorem 9 is summarized in Remark 10. In Theorem 11, the decreasing of the trace ratio sequence {F(wi)} means that the separability of the data set on the each newly extracted features decreases compared with the previous extracted features. This is the main result and contribution of this paper. Note that the trace ratio cost function is to evaluate whether a feature is good or not. By employing our proposed method, one can obtain m(m≤d and m can be greater than C-1) new extracted orthogonal features in which the wi is better than or at least equal to wi+1. Theorem 12 gives us more concrete conclusions: if we already have i features denoted as w1,…,wi, the i+ith feature wi+1 extracted by GLDA-TRA is the optimal feature in the space orthogonal to the space spanned by w1,…,wi. We present and prove Lemmas 4–8 as they are needed in proving Theorems 9–12. The theorems will be illustrated and verified in Section 4.


## 3.1. Proposed Algorithm

For convenience, we present the definition of a generalized eigenvalue as follows.Definition 2.
A numberλ is called a generalized eigenvalue of matrix B with respect to A if λ satisfies that Bx=λAx for a nonzero vector x, where A is a positive definite symmetrical matrix. When A=I, λ is a normal eigenvalue of B.Now the algorithm is given as follows.Initialization Step: construct symmetrical matricesSW1=SW∈Rd×d and SB1=SB∈Rd×d as shown in (1) based on the available learning data with W0 as an empty matrix.Stepi  (i=1,…,m) is as follows.(a)
Calculation stage: find a unit vector in the direction of a generalized eigenvector which corresponds to the maximum eigenvalue of the generalized eigenvalue problem ofSBi with respect to SWi. This can be achieved by the following process:do the Cholesky decompositionSWi=GiGi′. Let Si=Gi-1Bi-1(Gi-1)′ and obtain its maximum eigenvalue λi together with its corresponding eigenvector xi. Choose wi=(Gi-1)′xi/∥(Gi-1)′xi∥2 and Wi=[Wi-1wi].(b)
Update stage is as follows:(6)SWi+1=(I-WiWi+)SW(I-WiWi+)+μWiWi+SWWiWi+,SBi+1=(I-WiWi+)SB(I-WiWi+),where u is a sufficiently small positive number.Replacei by i+1 until m features are extracted in W.Remark 3.
(1) AsW is orthonormal, (6) can be rewritten as the following recursive from:
(7)SWi+1=(I-wiwi+)SWi(I-wiwi+)+μ∑j=1i‍wjwj+SWwjwj+,SBi+1=(I-wiwi+)SBi(I-wiwi+).
Later in Lemma 5, it will be shown that SWi for i=1,…,m can still be positive definite for a sufficiently small positive u. At each step i, SWi is positive definite and the generalized eigenvalue of SBi with respect to SWi exists.
(2) The reason why we need to find the maximum eigenvalue ofSi in the Calculation stage is shown in Lemma 8.
(3) Suppose thati features w1,…,wi,(i<d) have been extracted. Then Theorem 12 shows that wi+1 is found to ensure the trace ratio criterion function in (4) to attain its maximum value.
(4) Whenm=1, the LDA and GLDA-TRA extract the same feature.
(5) GLDA-TRA extractsm≤d features one by one. When i=m=d, SBi+1 becomes a zero matrix, and the algorithm will not extract any more features.


## 3.2. Derivation and Analysis of GLDA-TRA

Lemma 4.
For an arbitrary full column rank matrixWi∈Rd×i, WiWi+=Wi(Wi′Wi)-1Wi′ and I-WiWi+ are projection operators which project a vector onto span⁡{Wi} and span⁡{Wi}⊥, respectively.Proof.
Suppose thatv∈Rm×1 belongs to span⁡{Wi}; there exists a vector x∈Ri×1 such that v=Wix. It can be obtained that WiWi+v = Wi(Wi′Wi)-1Wi′Wix = Wix=v and (I-WiWi+)v = (I-Wi(Wi′Wi)-1Wi′)v = (I-Wi(Wi′Wi)-1Wi′)Wix=0. This lemma holds.Lemma 5.
S
W
i
+
1 in (6) for i=1,…,m are positive definite for a sufficiently small positive μ.Proof.
Note thatxSW1x>0 for any nonzero x∈Rm×1 as SW1=SW. Let x1=(I-WiWi+)x and x2=WiWi+x. Then x=(I-WiWi+)x + WiWi+x=x1+x2. As x≠0, x1 and x2 cannot be zero at the same time. Thus
(8)x′SWi+1x=x′(I-WiWi+)SW(I-WiWi+)x+μx′(WiWi+)SW(WiWi+)x=((I-WiWi+)x)′SW(I-WiWi+)x+μ(WiWi+x)′SWWiWi+x=x1′SWx1+μx2′SWx2>0,
for a sufficiently small positive number μ.Lemma 6.
For matricesW∈Rd×m, X∈Rd×d, define g(W)=tr⁡(W′XW). Then dg/dW=(X′+X)W.Proof.
LetW=[w1⋯wi⋯wm] and X=[x1⋯xd], where wi∈Rd×1 for i=1,…,m, and xi∈Rd×1 for i=1,…,d. We have
(9)g(W)=tr⁡(W′XW)=tr⁡(∑i=1dw1′xiwi1+⋯+∑i=1dwk′xiwik+⋯+∑i=1dwm′xiwi1)=∑i=1d∑j=1dwj1xjiwi1+⋯+∑i=1d∑j=1dwjkxjiwik+∑i=1d∑j=1dwjmxjiwim.
Then, for k0=1,…,m, j0=1,…,d, we get
(10)dg(W)dwj0k0=∑i=1dxj0iwik0+∑j=1dwjk0xjj0=∑i=1dxj0iwik0+∑i=1dxij0wik0.
So (10) gives dg/dW=(X′+X)W.Lemma 7.
Letw∈Rd×1 and the trace ratio function f(w)=(1/2)(tr⁡(w′SWw)/tr⁡(w′SBw)). The gradient ∇f(w)=df(w)/dw = (wSW·tr⁡(w′SBw)-wSB·tr⁡(w′SWw))/(tr⁡(w′SBw))2.Proof.
Letg1(w)=tr⁡(w′SWw) and g2(w)=tr⁡(w′SBw). From Lemma 6, ∇f(w)=df(w)/dw = ((dg1(w)/dw)g2(w)-(dg2(w)/dw)g1(w))/(g2(w))2 = (wSW·tr⁡(w′SBw)-wSB·tr⁡(w′SWw))/(tr⁡(w′SBw))2.DefineR(w)=tr⁡(w′SBiw)/tr⁡(w′SWiw), and we have the following results mentioned in the Calculation stage of step i.Lemma 8.
Assume that the maximum eigenvalue ofSi is λ* with an eigenvector x*. Then the unit vector w*=(Gi-1)′x*/∥(Gi-1)′x*∥ is an extracted feature ensuring that R(w) attains its maximum value which is equal to λ*.Proof.
Note thatf(w)=(1/2)(1/R(w))=(1/2)(tr⁡(w′SWiw)/tr⁡(w′SBiw)). Thus maximizing R(w) is equivalent to minimizing f(w). From Lemma 7, we have
(11)∇f(w(k))=df(w)dw|w=w(k)=wSWi·tr⁡(w′SBiw)-wSB·tr⁡(w′SWiw)(tr⁡(w′SBiw))2|w=w(k).
Then one can minimize f(w) by using an iterative method. Let w(k+1)-w(k)=-η∇f(w(k)), where k denotes the kth iteration and η is a small positive constant. Then f(w(k+1)) = f(w(k))-η(∇f(w(k)))′(∇f(w(k))) + o(∥w(k+1)-w(k)∥22)≤f(w(k)). Thus {f(w(k))} is a nonincreasing positive sequence and its limit, denoted as f*(w), exists. We have lim⁡k→∞f(w(k+1))=f(w(k))=f*(w). This implies that lim⁡k→∞η(∇f(wk))′(∇f(w(k)))=0. As η>0, then {w(k)} converges to an accumulation point w¯ that satisfies ∇f(w)|w=w¯=0. This gives
(12)df(w)dw|w=w¯=w¯SWi·tr⁡(w¯′SBiw¯)-w¯SBi·tr⁡(w¯′SWiw¯)[tr⁡(w¯′SBiw¯)]2=0.
Let λ=tr⁡(w¯′SBiw¯)/tr⁡(w¯′SWiw¯). From (12), we have SBiw¯=λSWiw¯. That is, λ is the generalized eigenvalue obtained from
(13)SBiw=λSWiw.
Note that each eigenvector w in (13) is an accumulation point of f(w), since it satisfies (12). Now we hope to convert the generalized eigenvalue problem of (13) to a normal eigenvalue problem. This is achieved by Cholesky decomposition of SWi, which is given as SWi=GiGi′, where Gi is a full column rank lower triangular matrix. By doing this, (13) becomes
(14)SBiw=λGiGi′w.
Defining x=Gi′w and substituting x into (14), it can be obtained that Six=λx, where Si=Gi-1SBi(Gi-1)′ is a symmetrical positive definite matrix. Let the maximum eigenvalue and its corresponding eigenvector of Si be λ* and x*. Finally the optimal w*=(Gi-1)′x*/∥(Gi-1)′x*∥ is a unit vector and λ*=max⁡⁡(R(w))=tr⁡(w*′SBw*)/tr⁡(w*′SWw*).Theorem 9.
Suppose thatwi is a feature extracted at step i; the trace ratio (1/2)(tr⁡(wi′SBi+1wi)/tr⁡(wi′SWi+1wi))=0 at step i+1.Proof.
Note thattr⁡(wi′SBi+1wi)=wi(I-WiWi+)SB(I-WiWi+)wi. We have tr⁡(wi′SWi+1wi)>0 since SWi+1 is positive definite by Lemma 5. Also, (I-WiWi+)wi=0 as (I-WiWi+)wi is a projection operator which projects a vector onto span⁡{Wi}⊥ based on Lemma 4. Thus, (1/2)(tr⁡(wi′SBi+1wi)/tr⁡(wi′SWi+1wi))=0.Remark 10.
As seen in Theorem9, if the feature wi has been extracted at step i, then wi is supposed to make the trace ratio (1/2)(tr⁡(wSBiw)/tr⁡(wSWiw)) attain its maximum in this step. While at step i+1, after SWi and SBi are updated to SWi+1 and SBi+1, respectively, we have (1/2)(tr⁡(wSBi+1w)/tr⁡(w′SWi+1w))=0 at w=wi. This means that the algorithm needs to find wi+1 which can maximize (1/2)(tr⁡(wSBi+1w)/tr⁡(w′SWi+1w)), and, obviously, wi+1 must be different from wi.Theorem 11.
Sequence{F(wi)}i=1m produced by GLDA-TRA is a decreasing sequence.Proof.
Assume thatVi is a space which the ith feature wi belongs to. That is, wi is the feature extracted from Vi such that F(wi) attains its maximum. Now we consider wi∈Vi at step i and wi+1∈Vi+1 at step i+1 with the respective corresponding maximum eigenvalues λi and λi+1. As SWi+1 is positive definite; from Lemma 4 it can be obtained that
(15)∥SWi+1w∥2=∥((I-WiWi+)SW(I-WiWi+)+μWiWi+SWWiWi+)w∥2=∥(I-WiWi+)SW(I-WiWi+)w∥2+μ∥WiWi+SWWiWi+w∥2⟶∥(I-WiWi+)SW(I-WiWi+)w∥2,
when μ→0. If w∈span⁡{Wi}, (I-WiWi+)w=0. Then, ∥SWi+1w∥2 = ∥(I-WiWi+)SW(I-WiWi+)w∥2=0. If w∈span⁡{Wi}⊥, (I-WiWi+)w=w and then ∥SWi+1w∥2>0. Also, as λi and λi+1 correspond to the maximum eigenvalues at steps i and i+1, it is obvious that
(16)∥λiSWiwi∥2≥μ∥wi∥2>0,∥λi+1SWi+1wi+1∥2≥μ∥wi+1∥2>0.
Thus, it can be concluded that wi∈span⁡{Wi} while wi+1∈Vi+1⊆span⁡{Wi}⊥. Similarly, wi∈Vi⊆span⁡{Wi-1}⊥. As Wi=[w1,…,wi], span⁡{Wi} increases as i increases. Thus, Vi decreases as i increases. In addition, we have V1=Rd=V2⊕span⁡{w1} = V3⊕span⁡{w1}⊕span⁡{w2} = ⋯ = Vm⊕span⁡{w1}⊕⋯⊕span⁡{wm-1}. Thus, F(w1)≥F(w2)≥⋯≥F(wm).Theorem 12.
The featurewi+1 extracted by GLDA-TRA is the optimal feature maximizing F(Wi+1) in (4) when Wi=[w1⋯wi] is extracted.Proof.
Note thatwi+1∈span⁡{Wi}⊥ as wi+1 is orthonormal to span⁡{Wi}. Consider an arbitrary w~i+1≠wi+1∈span⁡{Wi}⊥. Construct Wi+1=[Wiwi+1] and W~i+1=[Wiw~i+1]. Note that Wi is an orthonormal matrix. Then wi0′wj0=0 as long as i0≠j0 and we have
(17)F(Wi+1)-F(W~i+1)=wi+1SBwi+1+⋯+wi+1SBwi+1+wi+1SBwi+1wi+1SWwi+1+⋯+wi+1SWwi+1+wi+1SWwi+1-wi+1SBwi+1+⋯+wi+1SBwi+1+w~i+1SBw~i+1wi+1SWwi+1+⋯+wi+1SWwi+1+w~i+1SWw~i+1=(wi+1SBwi+1w~i+1SWw~i+1-w~i+1SBw~i+1wi+1SWwi+1)×((w1SBw1)2+⋯+(wiSBw1)2+(wi+1SBwi+1)(w~i+1SWw~i+1)((w1SBw1)2+⋯+(wiSBw1)2)-1.
Denote that SB=(I-WiWi+)SB+WiWi+SB and SW = (I-WiWi+)SW + WiWi+SW. Then, WiWi+wi+1=0 and (I-WiWi+)wi+1=wi+1. Thus, we have
(18)wi+1SBwi+1=wi+1((I-WiWi+)SB+WiWi+SB)wi+1=wi+1((I-WiWi+)SB(I-WiWi+))wi+1=wi+1SBi+1wi+1,
where SBi+1=((I-WiWi+)SB(I-WiWi+)). Similarly, it can be obtained that
(19)wi+1SWwi+1=wi+1SWi+1wi+1,w~i+1SBw~i+1=w~i+1SBi+1w~i+1,w~i+1SWw~i+1=w~i+1SWi+1w~i+1,
where SWi+1=((I-WiWi+)SW(I-WiWi+)). From Lemma 8, it is noted that tr⁡(wSBi+1w)/tr⁡(wSWi+1w) attains its maximum if and only if w is parallel to the eigenvector corresponding to the maximum generalized eigenvalue of SBi+1 with respect to SWi+1. Since wi+1 is parallel to the generalized eigenvector, as long as wi+1≠w~i+1, we have
(20)tr⁡(wi+1SBi+1wi+1)tr⁡(wi+1SWi+1wi+1)>tr⁡(w~i+1SBi+1w~i+1)tr⁡(w~i+1SWi+1w~i+1).
That is to say, wi+1SBwi+1w~i+1SWw~i+1>w~i+1SBw~i+1wi+1SWwi+1. So we have F(Wi+1)>F(W~i+1). Thus, this theorem holds.Remark 13.
Our main theoretical conclusions of this paper are shown in Theorems9–12. The implication of Theorem 9 is summarized in Remark 10. In Theorem 11, the decreasing of the trace ratio sequence {F(wi)} means that the separability of the data set on the each newly extracted features decreases compared with the previous extracted features. This is the main result and contribution of this paper. Note that the trace ratio cost function is to evaluate whether a feature is good or not. By employing our proposed method, one can obtain m(m≤d and m can be greater than C-1) new extracted orthogonal features in which the wi is better than or at least equal to wi+1. Theorem 12 gives us more concrete conclusions: if we already have i features denoted as w1,…,wi, the i+ith feature wi+1 extracted by GLDA-TRA is the optimal feature in the space orthogonal to the space spanned by w1,…,wi. We present and prove Lemmas 4–8 as they are needed in proving Theorems 9–12. The theorems will be illustrated and verified in Section 4.


## 4. Simulation and Discussion

To simply illustrate our idea, two experiments are done. The first one is on Iris data set and the second one is on an artificial data set.Example 14.
Iris data set is a standard data set to verify the performance of classification algorithms. There are150 data points belonging to three classes (C=3): Setosa, Versicolor, and Virginica, respectively. Each class has 50 samples with four features (d=4): Sepal length, Sepal width, Petal length, and Petal width.
One can constructSW and SB based on the Iris data set. It can be checked that rank⁡(SW)=4 and rank⁡(SB)=C-1=2, so LDA can only extract at most two features, while by employing GLDA-TRA to the Iris data set, we can extract m(m≤4) features one by one, which are denoted as W=[w1w2⋯wm]. When m=4, the obtained W=[w1w2w3w4] is given as
(21)W=[-0.08280.17170.6691-0.7183-0.4612-0.00120.62220.63250.4446-0.84340.30040.02700.76340.50900.27350.2886].
Obviously,W′W=I. Figure 1 shows the distribution of the three classes of the data points after projecting them onto each extracted feature w1,w2,w3, and w4. As seen in Figure 1, the separability of the data decreases in the directions of w1,w2,w3, and w4. One can also notice that the sequence {F(wi)}i=14 is strictly decreasing as shown in Figure 2.
When using the extracted features to do classification via support vector machine(SVM), both LDA and GLDA-TRA can obtain good results in this example as the data points are well separated even if they are projected to one-dimensional space. When we extract one feature, the accuracy is 98% (147/150) for both LDA and GLDA-TRA, since they extract the same feature when m=1  (see Remarks 1 and 3). When we extract m features with our methods for m=2,3, and 4, the accuracies are, respectively, 98% (147/150) when m=2, 98.6667% (148/150) when m=3, and 98% (147/150) when m=4. More details are shown in Table 1 and even when the data points are well separable, it is still observed that GLDA-TRA can be better than LDA in this case.Table 1
Comparison of LDA and GLDA-TRA in the Iris data set.
Name
Number of extracted features
Accuracy rate
Name
Accuracy rate
LDA
1
98%
GLDA-TRA
98%
LDA
2
98%
GLDA-TRA
98%
LDA
3
Not available
GLDA-TRA
98.6667%
LDA
4
Not available
GLDA-TRA
98%Illustration of projecting data points onto the extracted feature in order.
(a)
(b)
(c)
(d)Figure 2
The trace ration for each extracted feature in order.Example 15.
This radar data was a 34-dimensional data set(x∈R34, y∈Y={1,2}, C=2, d=34) collected by a system in Goose Bay, Labrador (http://archive.ics.uci.edu/ml/datasets/Ionosphere). This system consists of a phased array of 16 high-frequency antennas with a total transmitted power on the order of 6.4 kilowatts. The targets were free electrons in the ionosphere. “Good” radar returns are those showing evidence of some type of structure in the ionosphere. “Bad” returns are those that do not; their signals pass through the ionosphere.
The comparisons of LDA and GLDA-TRA in this example are shown in Figure3. It is seen that LDA and GLDA-TRA are the same when m=1. But LDA cannot continue to extract more features while the proposed GLDA-TRA still extracts m≤d features one by one. When m>12, the classification accuracy attains its maximum, which means that the dimension of the radar data set can be reduced to as low as 12 dimensions without losing any classification information.
In Theorem11, it is shown that the trace ratio sequence {F(wi)} on each extracted feature in order is decreasing; this can also be verified in Figure 4, where both trace ratio and its logarithm are plotted. And from Theorem 12, we know that each newly extracted feature is the optimal feature that maximizes the trace ratio function in the subspace orthogonal to the space spanned by the previous extracted features, which is exactly what we claimed in the abstract.Figure 3
Classification accuracy rate for LDA and GLDA-TRA in Example15.The trace ratio (a) and logarithm of the trace ratio (b) on each extracted feature in order.
(a)
(b)Example 16.
In the previous examples, both LDA and GLDA-TRA perform quite well even when only one feature is extracted. In this example, we will show that the data points are inseparable if they are projected to one-dimensional space. The classification problem is given by(22)(x,y)∈R7×Y={∥x∥2=|1+vi|ify=1∥x∥2=|4+vi|ify=2,
where y∈Y={1,2} and vi is a variable following normal distribution N(0,1). It can be known that most data points with label 1 locate around the surface of a sphere with radius y=1 while data points with label y=2 mostly locate around the surface of a sphere with radius 2.
To do the experiment, two hundred data points that are equally distributed in above two classes are generated. It is observed that LDA does not perform well if it is used to extract features in this problem. This is because LDA can only extract one feature and the data points are inseparable or well separable in an arbitrary one-dimensional feature space. By using GLDA-TRA, when we extract two features (m=2), we obtain a projection matrix W∈R34×2. By employing x~=WTx, one can then visualize the data points by projecting the date onto the two extracted features x~. Figure 5 shows the classification using SVM on this two-dimensional extracted feature space. With GLDA-TRA, the data points can be better separated if we extract more (m>1) features. The comparisons of LDA and GLDA-TRA regarding the classification accuracy rate are shown in Table 2 and Figure 6.Table 2
Comparison of LDA and GLDA-TRA in Example16.
Name
Number of extracted features
Accuracy rate
Name
Accuracy rate
LDA
1
69%
GLDA-TRA
69%
LDA
2
Not available
GLDA-TRA
77.5%
LDA
3
Not available
GLDA-TRA
88%
LDA
4
Not available
GLDA-TRA
92.5%
LDA
5
Not available
GLDA-TRA
96%
LDA
6
Not available
GLDA-TRA
97%Figure 5
The data points projected onto the first two extracted features in Example16.Figure 6
Classification accuracy rate for LDA and GLDA-TRA in Example16.


## 5. Conclusion

In this paper, a generalized linear discriminant analysis based on trace ratio criterion (GLDA-TRA) algorithm has been proposed. This is to overcome the problem that linear discriminant analysis (LDA) can only extract limited features in classification. It is shown that, in GLDA-TRA, a set of orthogonal features can be extracted one by one. Each newly extracted feature is the optimal feature that maximizes the trace ratio criterion function in the subspace orthogonal to the space spanned by the previous extracted features. Finally the extracted features are such that the trace ratio sequence of these features is decreasing in order. Experimental results also show the effectiveness of our proposed algorithm.



---

*Source: 725204-2014-05-29.xml*