Joint Feature Selection and Subspace Learning

1 Introduction

High-dimensional data in the input space is usually not good for classification due to the curse of dimensionality. A common way to resolve this problem is dimensionality reduction, which has attracted much attention in machine learning community in the past decades. Generally speaking, dimensionality reduction techniques can be classified into two categories: (1) feature selection [Guyon and Elisseeff, 2003]: to select a subset of most representative or discriminative features from the input feature set, and (2) subspace learning [Belhumeur et al., 1997][He and Niyogi, 2003][He et al., 2005][Yan et al., 2007] (a.k.a feature transformation): to transform the original input features to a lower dimensional subspace.

The most popular subspace learning methods include Principal Component Analysis (PCA) [Belhumeur et al., 1997], Linear Discriminant Analysis (LDA) [Belhumeur et al., 1997], Locality Preserving Projection (LPP) [He and Niyogi, 2003] and Neighborhood Preserving Embedding (NPE) [He et al., 2005]. Despite different motivations of these methods, they can all be interpreted in a unified Graph Embedding framework [Yan et al., 2007].

One major disadvantage of the above methods is that the learned projection is a linear combination of all the original features, thus it is often difficult to interpret the results. Sparse subspace learning methods attempted to solve this problem. For example, [Zou et al., 2004] proposed a sparse PCA algorithm based on ℓ2-norm and ℓ1-norm regularization. [Moghaddam et al., 2006] proposed both exact and greedy algorithms for binary class sparse LDA as well as its spectral bound. [Cai et al., 2007] proposed a unified sparse subspace learning (SSL) framework based on ℓ1-norm regularized Spectral Regression.

However, the selected features by sparse subspace methods are independent and generally different for each dimension of the subspace. See Figure 1 (a) for an illustrative toy example of the projection matrix learned by SSL. Each row of the projection matrix corresponds to a feature, while each column corresponds to a dimension of the subspace. We can see that for the first dimension of the subspace, the 3rd and 6th features are not selected, while for the second dimension of the subspace, the selected features are all except the 1st and 4th one. Hence it is still unclear which features are really useful as a whole. Our goal is to learn a projection matrix like Figure 1 (b), which has row-sparsity (elements in a row are all zero). Hence it is able to discard the irrelevant features (e.g., the 1st, 5th and 7th features) and transform the relevant ones simultaneously. One intuitive way is performing feature selection [Guyon and Elisseeff, 2003] before subspace learning. However, since these two sub-processes are conducted individually, the whole process is likely suboptimal.

Based on the above motivation, in this paper, we aim to jointly perform feature selection and subspace learning. To achieve this goal, we reformulate subspace learning as solving a linear system equation, during which we use L2,1norm on the projection matrix, encouraging row-sparsity. It is worth noting that L2,1-norm has already been successfully applied in Group Lasso [Yuan et al., 2006], multi-task feature learning [Argyriou et al., 2008], joint covariate selection and joint subspace selection [Obozinski et al., 2010]. The resulted optimization problem includes two situations, for each of which we present a very simple algorithm, that is theoretically guaranteed to converge. Experiments on benchmark face recognition data sets demonstrate the effectiveness of the proposed framework.

The remainder of this paper is organized as follows. In Section 2, we briefly introduce the graph embedding view of subspace learning. In Section 3, we present a framework for joint feature selection and subspace learning. In Section 4, we review some related works. Experiments on benchmark face recognition data sets are demonstrated in Section 5. Finally, we draw a conclusion and point out future work in Section 6.

1.1 Notations

Given a data matrix X = [x1, . . . , xn] ∈ Rd×n, we aim to learn a projection matrix A ∈ Rd×m, projecting the input data into an m-dimensional subspace. For a matrix A ∈ Rd×m, we denote the ith row of A by ai , and the jth column of A by aj . The Frobenius norm of A is defined as ||A||F = √∑di ||ai ||22 , and the L2,1-norm of A is defined as ||A||2,1 = ∑di ||ai ||2.

2 Graph Embedding View of Subspace Learning

Many dimensionality reduction methods have been proposed to find low-dimensional representation of xi. Despite different motivations of these methods, they can be nicely interpreted in a general graph embedding framework [Yan et al., 2007]. In graph embedding, we construct a data graph G whose vertices correspond to {x1, . . . , xn}. Let W ∈ Rn×n be a symmetric adjacency matrix with Wij characterizes the favorite relationship among the data. The purpose of graph embedding is to find the optimal low-dimensional vector representation for the vertices of graph G that best preserves the relationship between the data points. In this paper, we focus on linear dimensionality reduction. That is, XT A. The optimal A is given by the following optimization problem,

where Dii = ∑j Wij is a diagonal matrix, and L = DW is called graph Laplacian [Chung, 1997], I is the identity matrix with proper size. The above problem can be solved by generalized eigen-decomposition XLXT A = ΛXDXT A, where Λ is a diagonal matrix whose diagonal elements are eigenvalues.

With different choices of W, the linear graph embedding framework leads to many popular linear dimensionality reduction methods, e.g. PCA [Belhumeur et al., 1997], LDA [Belhumeur et al., 1997], LPP [He and Niyogi, 2003] and NPE [He et al., 2005]. We briefly give two examples below.

LDA: Suppose we have c classes and the k-th class have nk samples, n1 + · · · + nc = n. Define

LPP [He and Niyogi, 2003]: Define

where Nk(xi) denotes the set of k nearest neighbors of xi, d(xi, xj) measures the similarity between xi and xj , which can be chosen as Gaussian kernel e‖xi-xj ‖2 2σ2 or cosine distance xTi xj ‖xi‖‖xj‖ . For more examples and other extensions of graph embedding, please refer to [Yan et al., 2007].

3 Joint Feature Selection and Subspace Learning

Since each row of the projection matrix corresponds to a feature in the original space, in order to do feature selection, it is desirable to have some rows of the projection matrix be all zeros. This motivates us to use L2,1-norm on the projection matrix, which leads to row-sparsity of the projection matrix. As a result, based on Eq. (1), we formulate joint feature selection and subspace learning as follows,

where μ is a regularization parameter. Although the objective function is convex, the constraint is not. Hence it is difficult to optimize. In the following, we will reformulate the problem to make it easy to be solved.

Theorem 3.1. Let Y ∈ Rn×m be a matrix of which each column is an eigenvector of eigen-problem Wy = λDy. If there exists a matrix A ∈ Rd×m such that XT A = Y, then each column of A is an eigenvector of eigen-problem XWXT a = λXDXT a with the same eigenvalue λ.

Theorem 3.1 shows that instead of solving the eigenproblem XWXT A = ΛXDXT A, A can be obtained by the following two steps:

1. Solve the eigen-problemWY = ΛDY to get Y;
2. Find Awhich satisfies XT A=Y.

Note that only the second step involves A. XT A = Y is a linear system problem, which may behave in any one of three possible ways: (1) The system has infinitely many solutions; (2) The system has a single unique solution; and (3) The system has no solution. In the rest of this section, we will discuss (1) in one situation and (2) (3) in another situation.

3.1 Situation 1

When the linear system has infinitely many solutions, we formulate the proposed method as

The optimization problem in Eq. (5) is similar to the problem appeared in Multiple Measurement Vector model in signal processing [Sun et al., 2009]. We derive a very simple algorithm in the sequel.

The lagrangian function of the problem in Eq. (5) is

Taking the derivative of L(A) with respect to A1 , and setting the derivative to zero, we get

where G is a diagonal matrix with the i-th diagonal element equal to

Left multiplying the two sides of Eq. (7) by XT G-1 , and using the constraint XT A = Y, we have

Note that when gii = 0, G-1 cannot be computed. This is handled by adding a smoothing term ϵ to gii. It can be justified that the algorithm converges to the global solution if decreasing value of ϵ is used.

Substituting Eq. (9) into Eq. (7), we get

In summary, we present the algorithm for optimizing Eq.(5) in Algorithm 1.

Algorithm 1 Joint Feature Selection and Subspace Learning (Situation 1)

Initialize:G0 = I, t = 0; Compute Y based on WY = ΛDY;

repeat

Compute At+1 = G-1

t X(XT G-1

Compute Gt+1 based on At+1;

t = t + 1;

until convergence

The convergence of this algorithm was proved in [Nie et al., 2010]. In addition, there is no additional parameter introduced besides the parameters that are needed to construct the graph as in traditional linear graph embedding. This is a very appealing property.

1 For ||A||2,1, we compute its sub-gradient because it is not smooth.

t X)-1Y;

3.2 Situation 2

When the linear system has a single unique solution, then it is also the solution of Eq. (5). However, it does not have row sparsity. When the linear system has no solution, then the formulation in Eq. (5) will have no solution either. In both of these cases, we turn to solve the constrained problem as follows,

Or equivalently the regularized problem,

This is a Group Lasso problem [Yuan et al., 2006]. When μ = ∞, the optimization problem in Eq. (12) degenerates to that in Eq. (5). Note that it is difficult to give an analytical relationship between δ and μ. However, such a relationship is not crucial for our problem. The objective function in Eq. (12) is a non-smooth but convex function. In the following, we will derive an algorithm which is similar to Algorithm 1 for solving Eq. (12).

Taking the derivative of Eq. (12) with respect to A, and setting the derivative to zero, we get

which leads to

where G is defined in Eq. (8). Since G also depends on A, the above closed-form expression of the optimal A is fundamentally a fixed-point iteration.

According to the Woodbury matrix identity [Golub and Loan, 1996]

we can further get

It is worth noting that if μ = ∞, Eq. (16) reduces to Eq. (10). This shows the relation between Eq. (5) and Eq. (12) again.

In summary, we present the algorithm for optimizing Eq.(12) in Algorithm 2. The convergence of this algorithm can be proved similarly to [Nie et al., 2010].

4 RelatedWork

In this section, we discuss some approaches which are closely related to our method.

In order to pursue interpretability, [Cai et al., 2007] proposed Sparse Subspace Learning (SSL), which is based on ℓ1-norm regularization on each column of the projection, i.e., a,

Algorithm 2 Joint Feature Selection and Subspace Learning (Situation 2)

Initialize:G0 = I, t = 0 and μ; Compute Y based on WY = ΛDY; repeat Compute At+1 = G-1 t X(XT G-1 t X+ 1 2μ I)-1Y; Compute Gt+1 based on At+1; t = t + 1; until convergence

where y is the eigenvector of Wy = λDy. Due to the nature of the ℓ1 penalty, some entries in a will be shrunk to exact zero if μ is large enough, which results in a sparse projection. However, SSL does not lead to feature selection, because each column of the projection matrix is optimized one by one, and their sparsity patterns are independent.

Most recently, [Masaeli et al., 2010] proposed Linear Discriminant Feature Selection, which modifies LDA to admit feature selection as follows,

where Sb is the between-class scatter matrix, and Sw is the with-in class scatter matrix [Belhumeur et al., 1997], ∑di=1 ||ai ||∞ is the ℓ1/ℓ∞ norm of A. Note that ℓ1/ℓ∞ has the similar nature of L2,1-norm, which leads to row-sparsity. The optimization problem is convex and solved by quasiNewton method [Boyd and Vandenberghe, 2004]. The main disadvantage of this method is that the evaluation of the gradient of tr((AT SwA)-1 (AT SbA)) is computationally very expensive. Hence it is limited to small-scale data.

5 Experiments

In this section, we evaluate two instances of our framework, FSSL with LDA-type adjacency matrix defined in Eq. (2), referred to as FSSL(LDA), FSSL with LPP-type adjacency matrix defined in Eq. (3), referred to as FSSL(LPP), and compare them with the state of the art subspace learning methods, e.g. PCA, LDA, LPP [He and Niyogi, 2003], and sparse subspace learning approaches, e.g. SSL with LDA-type adjacency matrix, denoted by SSL(LDA) [Cai et al., 2007], SSL with LPP-type adjacency matrix, denoted by SSL(LPP) [Cai et al., 2007]. We also compare it with a feature selection method, e.g., Fisher score (FS), and feature selection followed with subspace learning (FS+SL). In detail, we use Fisher score to do feature selection and LDA (or LPP) to do subspace learning, which are referred to as FS+SL(LDA) and FS+SL(LPP) respectively. We use 1-Nearest Neighbor classifier as baseline. All the experiments were performed in Matlab on a Intel Core2 Duo 2.8GHz Windows 7 machine.

5.1 Data Sets

We use two standard face recognition databases which are used in [Cai et al., 2007].

Extended Yale-B database contains 16128 face images of 38 human subjects under 9 pose and 64 illumination conditions. In our experiment, we choose the frontal pose and use all the images under different illumination, thus we get 2414 image in total. All the face images are manually aligned and cropped. They are resized to 32 × 32 pixels, with 256 gray levels per pixel. Thus each face image is represented as a 1024-dimensional vector.

CMU PIE face database [Sim et al., 2003] contains 68 individuals with 41368 face images as a whole. The face images were captured by 13 synchronized cameras and 21 flashes, under varying pose, illumination and expression. In our experiment, one near frontal poses (C27) are selected under different illuminations, lighting and expressions which leaves us 49 near frontal face images for each individual.

5.2 Parameter Settings

For both data sets, p = 10, 20, 30 images were randomly selected as training samples for each person, and the rest images were used for testing. The training set was used to learn a subspace, and the recognition was performed in the subspace by 1-Nearest Neighbor classifier. Since the training set was randomly chosen, we repeated each experiment 20 times and calculated the average recognition accuracy. In general, the recognition rate varies with the dimensionality of the subspace. The best average performance obtained as well as the corresponding dimensionality is reported in Table 1 and Table 2. We also report the computational time (in second) of subspace learning 2 .

For LDA, as in [Belhumeur et al., 1997], we first use PCA to reduce the dimensionality to n - c and then perform LDA to reduce the dimensionality to c - 1. This is also known as Fisher Face [Belhumeur et al., 1997]. For LPP, we use the cosine distance to compute the similarity between xi and xj. For FS+LDA and FS+LPP, we first use Fisher Score to select 50% features and then perform LDA (or LPP) to reduce the dimensionality. For SSL(LDA), SSL(LPP), we tune μ by searching the grid {10, 20, . . . , 100} according to [Cai et al., 2007]. For FSSL(LDA) and FSSL(LPP), when p = 10, the linear system has infinite solution, we run Algorithm 1. And when p = 20, 30, we run Algorithm 2, where we simply set μ = 0.1. The smoothing term ϵ is set to 0.01.

5.3 Results

The experimental results are shown in Table 1 and Table 2. We can observe that (1) SSL is better than the corresponding linear subspace learning method (e.g., SSL(LDA) > LDA), which implies sparse subspace learning is able to improve the classification performance of subspace learning method; (2) FSSL outperforms the corresponding SSL method overwhelmingly (e.g., FSSL(LDA) > SSL(LDA)), which indicates feature selection can improve the corresponding linear subspace learning method greatly; (3) Feature selection before subspace learning (FS+SL) sometimes achieves better results than subspace learning and even better than SSL. This implies the potential performance gain of doing feature selection during subspace learning. However, at more cases, FS+SL(LDA) and FS+SL(LPP) are worse than LDA and LPP. This is because feature selection and subspace learning are conducted individually in FS+SL. So the selected features are not necessary helpful for subspace learning. In contrast, FSSL perform feature selection and subspace learning simultaneously, the selected features are generally beneficial for subspace learning. Hence the proposed framework improves subspace learning consistently; and (4) The computational time of FSSL is much less than that of SSL.

2 It does not include the testing time of 1-Nearest Neighbor classifier.

5.4 Study on the Dimensionality of the Subspace

In this subsection, we get a closer look at the recognition accuracy with respect to the dimensionality of the learned subspace. Figure 2 shows the performance of all the subspace learning methods on the two databases with 20 training samples respectively, where the horizontal axis represents the dimensionality of the subspace, and the vertical axis denotes the average recognition accuracy of 20 independent runs.

It is shown that the improvement of our methods over sparse subspace learning methods and subspace space learning methods are consistent over a wide range of the dimensionality of the subspace, which strengthens the superiority of our approaches. Similar phenomenon can be observed when we use 10 and 30 training samples. For the space limit, we do not show them.

5.5 Projection Matrices & Selected Features

To get a better understanding of our approach, we plot the projection matrices of our method and related methods in Figure 3 (a),(b) and (c). Clearly, the projection matrix of LDA is not sparse, while each column of the projection matrix of SSL(LDA) is sparse. However, the sparsity patterns of each column are not coherent. In contrast, each row of the projection matrix of FSSL(LDA) tends to be zero or nonzero simultaneously, which benefits from the nature of L2,1-norm, and leads to feature selection.

From Figure 3 (d), we can see that the selected features (pixels) are asymmetric. In other word, if one pixel is selected, its axis symmetrical one will not be selected. This is because the face image is roughly axis symmetric, so one in a pair of axis symmetrical pixels is redundant given the other one is selected. Moreover, the selected pixels are mostly around the eyebrow, the corner of eyes, nose and cheek, which are discriminative for distinguishing face images of different people. This is consistent with our real-life expe-

��

rience.

6 Conclusion and Future Work

In this paper, we propose to do feature selection and subspace learning simultaneously in a joint framework, which is based on using L2,1-norm on the projection matrix, that achieves the goal of feature selection. Experiments on benchmark face recognition data sets illustrate the efficacy of the proposed framework. In our future work, we will study joint feature selection and nonlinear subspace learning in kernel space.

Acknowledgments

The work was supported in part by NSF IIS-09-05215, U.S. Air Force Office of Scientific Research MURI award FA9550-08-1-0265, and the U.S. Army Research Laboratory under Cooperative Agreement Number W911NF-09-2-0053 (NS-CTA).