Top-k Supervise Feature Selection via ADMM for Integer Programming

1 Introduction

In modern machine learning and pattern recognition applications, data are commonly represented by high dimensional feature vectors, such as image classification [Chatfield et al., 2011; Chang et al., 2016] and video recognition [Lan et al., 2015; Chang et al., 2017]. High dimensional data are not suitable for directly learning because the time cost and storage requirement will be very high. Furthermore, they usually contain many noise and redundant features [Peng et al., 2005] that could degrade the generalization capability of the learning algorithms. Feature Selection (FS) [Guyon and Elisseeff, 2003] is designed to identify the most relevant and important features that can give a compact and accurate data representation for learning. Compared with other feature analysis methods, such as feature extraction [Belhumeur et al., 1997], FS has better interpretability because it keeps the semantic meaning of the features. Also, the cost of feature collection can be reduced because that one only needs to collect the selected features in FS rather than all the features as in feature extraction does. As a result, FS has become a hot topic in machine learning and pattern recognition [Chang et al., 2014].

∗Corresponding author.

Based on the rule of the learning algorithm, there are roughly three types of FS methods in literature: the filtertype, the wrapper-type, and the embedded-type methods. The filter-type method evaluates all data features based on certain criteria, where no learning algorithm is involved. Representative filter-type feature selection methods include the reliefF [Kira and Rendell, 1992], mRMR [Peng et al., 2005], Fisher score [Duda et al., 2000], and Laplace score [He et al., 2006] methods. The wrapper-type method applies a classifier as a black box to score the features. The widely used wrappertype methods include the Support Vector Machine Recursive Feature Elimination (SVM-RFE) [Guyon et al., 2002] and the Correlation-based Feature Selection (CFS) [Hall and Smith, 1999]. The embedded-type method [Wang et al., 2007; Argyriou et al., 2007] embeds the feature selection procedure in a classifier algorithm and only a single optimization problem is involved. Because the embedded and wrapper methods interact with a learning algorithm, they tend to achieve better classification results than the filter-type method when some specific learner is involved for evaluation.

Recently, empirical studies of sparse representation and compressed sensing [Elad, 2010] indicate that sparsity is one of the basic and intrinsic properties of real world data. FS, which finds a sparse attributes to represent the input data, can be regarded as a natural application of the sparse representation theory. A large number of FS methods resort to the sparsity-inducing regularization terms/constraints, such as the ℓ0, ℓ1-norm, ℓ0,2-norm and ℓ1,2-norm based penalty terms/constraints, to achieve feature evaluation and selection (The appearance of the notations may differ with those used in previous papers, but they are essentially equivalent). From the sparsity perspective, ℓ0-norm and ℓ0,2-norm are more desirable to select the features because that they can induce the sparsest solution, i.e., each feature should be associated with either the zero score or a large score. However, ℓ0-norm and ℓ0,2-norm regularized/constrained optimization problems have been proved to be NP-hard [Amaldi and Kann, 1998] and are very difficult to solve. Fortunately, theoretical results show that, under mild conditions, ℓ1 and ℓ0 are essentially equivalent [Donoho, 2004], i.e., ℓ1-norm and ℓ1,2-norm can be regarded as efficient approximations to ℓ0norm and ℓ0,2-norm respectively. The ℓ1-norm based problem, also known as Lasso [Tibshirani, 1994], is commonly used in FS on binary class data sets. Destrero et al. utilize the Lagrangian form of Lasso for feature selection in face recognition [Destrero et al., 2009]. In [Zou and Hastie, 2005], the elastic net regularization is proposed to handle features with strong correlations. To remove the feature redundancy, group Lasso is introduced to integrate the feature structure and then evaluate the importance of features, where the structures include the disjoint groups [Zhang et al., 2012], the overlapping groups [Jenatton et al., 2011], and so on. To address multi-class problems, Nie et al. [Nie et al., 2010] propose to apply ℓ1,2-norm instead of ℓ1-norm as the penalty and have shown promising results. Many recent FS methods are proposed in the form of ℓ1,2-norm regularized/constrainted optimization problems [Xiang et al., 2012; Du and Shen, 2015; Han et al., 2015]. The matrix norm has been extended to ℓp,2, (p ∈ (0, 1]) [Wang et al., 2014; Tao et al., 2016] and ℓp,r, (1 < r) norms for robust FS. Because p = 0 is more desirable than any p > 0, an exact top-k FS via a optimization problem with the ℓ0,2-norm constraint is proposed in [Cai et al., 2013].

In this paper, we propose a novel efficient and robust supervised FS method, which has the following properties.

1. The proposed method is formulated as a problem with two group of variables and an ℓ0,2(or ℓ0)-norm constraint. The sub-problem with respect to one group of the variables in the model is a 0-1 integer programming and the sub-problem of the other group of the variables admits a closed-form solution.
2. The 0-1 integer programming of our method is firstly transformed into an equivalent optimization problem with two additional continuous constraints. The equivalent problem is shown able to be efficiently solved by the Alternating Direction Method of Multipliers (ADMM). The Matlab code is published online 1 .
3. We provide an efficient algorithm to solve ℓ0,2(or ℓ0)- constrained supervised FS method. The proposed method guarantees to select the top k features instead of k top features under the proposed criterion. Experiments show that the proposed method has superior performance than the compared state-of-the-art supervised FS methods.

The rest of this paper is structured as follows: in Section 2, we review some related sparsity-inducing FS methods. The proposed supervised FS method is described in Section 3. Experimental comparisons with state-of-the-art supervised FS methods on benchmark data sets are presented in Section 4. Finally, the conclusion is drawn in Section 5.

1 https://github.com/cxj273/IJCAI2017_1274

2 Sparsity-Inducing Supervised Feature Selection Background

2.1 Notations and Definitions

The ℓp-norm of a vector v ∈ RD is defined as ‖v‖p = (∑Di=1 |vi|p)1p , where vi denotes the i-th entry in v. The ℓ0-norm of a vector v is defined as ‖v‖0 = ∑Di=1 |vi|0 , i.e., the counts of nonzero entries in v. diag(v) (v ∈ RD ) is a diagonal matrix whose diagonal elements are the entries of vector v and diag(Θ) (Θ ∈ RD×D ) is a D-dimensional vector consists of the diagonal elements of the matrix Θ. The ℓp,r-norm of a matrix A ∈ RC×D is defined as

where Aj (j = 1, · · · ,D) denotes the j-th column of A. Consequently, the ‖A‖0,2 is naturally defined as ‖A‖0,2 = ∑Dj=1 ‖Ai‖02 , which counts the number of nonzero columns in A.

Let X = [x1, · · · , xN] ∈ RD×N be the input data matrix, where D is the input dimensionality, N is the number of data points, the i-th column, xi, denotes a data vector. For presentation clarity and simplicity, we assume that the constant value 1 has been added to the bottom of each data vector as an additional dimension and thus the bias term can be omitted throughout this paper.

2.2 Sparsity-Inducing Supervised Feature Selection

Given a binary class data set, many sparsity-inducing FS methods can be interpreted as the approximation or relaxed version of the following problem

where A ∈ RD , y ∈ RN and its i-th entry yi ∈ {0, 1} is the class label of xi. The ℓ0-norm constraint means only k entries are non-zeros and thus only k features of X are used. The purpose of FS is to find out which k features are effective. We can also write the problem (1) in the equivalent Lagrangian form:

where nonnegative γ is a given tradeoff parameter. However, the problems (1) and (2) have been proven to be NP-hard problems and are computationally infeasible. Many FS methods [Destrero et al., 2009; Cai et al., 2010] replace ℓ0-norm with its convex surrogate ℓ1-norm and have shown promising results. Generally, the values of entries in A indicate the importance of data features. The features with high values in A are then selected.

For multi-class data set, the ideal top-k FS [Cai et al., 2013] is to optimize

or its equivalent form

where A ∈ RC×D be the projection matrix, Y ∈ RC×N be the label matrix with the i-th column yi ∈ RC . If xi is in the k-th class, the k-th entry of yi is 1 and the rest entries are 0s. However, ℓ0,2-norm regularized/constrained optimization problem is still NP hard and is computationally infeasible. To address this problem, matrix ℓp,r-norm has been proposed for inducing column sparsity to achieve feature selection, where p ∈ (0, 1] and r > 1. The objective function of ℓ1,2 based Robust Feature Selection (RFS) [Nie et al., 2010] is proposed as

To enforce further sparsity, [Wang et al., 2014] proposes the ℓp,2-norm based FS method.

where p ∈ (0, 1] is a given parameter.

3 Top-k Supervised Feature Selection

3.1 Formulation

Both ℓp,2-norm and ℓp-norm are surrogates of their original sparsity-inducing constraints, ℓ0,2-norm and ℓ0-norm respectively. Therefore, it is more desirable to optimize the original sparsest problem instead of its relaxed formulation. In this paper, we propose a novel method that directly solves the sparsest feature selection model. The proposed method does not score every features, but selects the optimal feature set with k features under the proposed criteria. For convenience, the classical multi-class least squares regression model is utilized to learn the linear projection matrix A as

As can be seen, the proposed model is a classical regularized least squares regression with an additional ℓ0,2-norm based constraint. The ℓ0,2-norm constraint requires k and only k columns of the matrix A are not zero vectors. Therefore, the k corresponding features of X are active and the rest features hibernate in the regression (6). Instead of directly solving the ℓ0,2 constrained problem, we first transform it to an equivalent ℓ0-norm constrained problem as follows:

Once let ˆ A = Adiag(v) = [v1A1, · · · , vDAD], we can see that the FS problems (6) and (7) are essentially equivalent. With respect to the variable A, the problem (7) has the closedform least squares solution. However, the problem with respect to the variable v is a 0-1 integer programming and is generally very difficult to solve. Here we utilize the ℓ2-box method to efficiently address this problem [Wu and Ghanem,

2016]. The binary constraint {0, 1}D can be replaced with an equivalent set of continuous constraints, i.e., the intersection of a box and a shifted ℓ2-sphere. The result is presented in Proposition 3.1.

Proposition 3.1 [Wu and Ghanem, 2016] Let 1D ∈ RD be the vector whose entries are all 1s, we have

An illustrative example of the equivalence between the binary constraint and the continuous constraints in 2-D space is given in Figure 3.1. Based on Proposition 3.1, we can obtain the following problem which is equivalent to (7).

where the two sets Sb = {v : v ∈ [0, 1]D } and Sp = {v : ‖v - 1D 2 ‖22 = D4 }. In (8), the two continuous constraints in Proposition 3.1 are separated by two additional variables v1 and v2. The problem (8) can now be efficiently optimized by the ADMM method.

3.2 Optimization

ADMM method has been widely used in convex optimization, and there is also growing interests and applications on the advantage of ADMM in non-convex optimization. In this section, we study an efficient solution for our FS method (8) using the ADMM method, which solves the following subproblems iteratively:

1. Fix the other variables, optimize A through solving a classic least squares regression problem.
2. Given A, compute v through solving an unconstraint quadratic optimization problem.

1. Project v on to Sb and Sp to obtain v1 and v2 respectively.
2. Update the Lagrange multipliers.

Using a parameter ρ > 0, the augmented Lagrangian function of (8) is obtained as

where y1 ∈ RD , y2 ∈ RD and y3 ∈ R are Lagrange multipliers for the three equality constraints. The ADMM approach then iteratively optimize the variables individually. Denote by (A(t) , v(t) , v(t) 1 , v(t) 2 ) the optimization variables at iteration t, and by (y(t) 1 , y(t) 2 , y(t) 3 ) the Lagrange multipliers at iteration t.

Step 1: Solving the linear projection matrix A when other variables are fixed. The (9) with respect to A is a classical least squares regression problem and the solution can be directly provided as (at the t + 1-th iteration)

Step 2: Optimize v when other variables are fixed. After some mathematical deductions, we can obtain the following unconstraint quadratic optimization problem with respect to the variable v as

where Φ = XXT , Ψ(t+1) = (A(t+1) )T A(t+1) , Θ(t+1) = XY T A(t+1) . Imposing the derivative of the objective function with respect to v to zero, we obtain the closed-form solution as

Step 3: Update variables v1 and v2 through projections onto Sb and Sp. The variables are updated as follows

For any x, the projection on a box PSb is a element-wise function, which is given by

For any vector x ∈ RD , one first compute the two candidates of the projection on Sp as

Then, the projection of x on PSp can be computed as

Step 4: Update variables y1, y2, y3 and ρ. Having variables (A(t+1) , v(t+1) , v(t+1) 1 , v(t+1) 2 ) fixed, the Lagrange multipliers and step size ρ are updated as follows

where µ > 1 is a given parameter.

These steps are repeated until convergence is achieved or the number of iterations exceeds a maximum iteration number. Convergence is achieved when we have ‖v(t+1) - v(t+1) 1 ‖∞ ≤ ε, ‖v(t+1) - v(t+1) 2 ‖∞ ≤ ε, and |1DT v(t+1) - k| ≤ ε. The updates for the details of ADMM implementation are summarized in Algorithm 1.

Algorithm 1 ADMM for solving problem (7)

Input: Data matrix X, label matrix Y , γ;

A is initialized as the identity matrix I, v = 1D, v1 = v2 = 0D, ρ = 1, and µ = 1.05

Output: Projection matrix A and vector v

• 1: while not converged do
• 2: Update A(t+1) as in (10);
• 3: Update v(t+1) as in (11);
• 4: Update v(t+1) 1 and v(t+1) 2 through projections onto Sb and Sp as in (12);
• 5: Update y(t+1) 1 , y(t+1) 2 , y(t+1) 3 and ρ as Eq. (13).
• 6: If not converged, set t ← t + 1.
• 7: end while

3.3 Computational Complexity Analysis

To optimize the objective function of the proposed method, the step 1 and step 2 are time consuming operations. At step 1, the algorithm needs to inverse an D × D matrix and the time complexity is O(D3 ). At step 2, an matrix inversion is computed on a D × D matrix whose time complexity is O(D3 ). Assuming there are T iterations before the algorithm stops, the total cost of the proposed method is O(T(D3 )). As can be seen, the computation complexity is comparably higher than most of the sparsity-inducing supervised FS methods. We plan to reduce the computation complexity of the proposed method in the future work.

4 Experiments

In this section, the proposed method is compared with stateof-the-art supervised feature selection methods on benchmark image data sets. The experiments include the supervised classification by the Nearest Neighbor classifier (1-NN) and the Support Vector Machine (SVM) under various experimental settings. The numerical convergence analysis of the proposed method is also included.

4.1 Datasets Description

Three real world data sets are used in our experiments. The important statistics of these data sets are briefly summarized as bellow:

• The Coil-20 data set 2 contains 1440 image samples from 20 classes and each image is transformed into a 1024dimensional data point. There are 72 samples in each class.
• The MNIST handwritten digital image data set 3 has 6996 data points of digits '0' - '9'. Each sample is a 784 dimensional feature vector.
• There are 2114 frontal-face images of 38 individuals in the Yale-B face image data set 4 . Each image is stacked to a 1024-dimensional data vector.

2 http://www.cs.columbia.edu/CAVE/software/softlib/coil-

20.php

3 http://www.escience.cn/people/fpnie/

4 http://www.cad.zju.edu.cn/home/dengcai/Data/FaceData.html

4.2 Experimental Setup

The following state-of-the-art supervised feature selection algorithms are compared in this paper. The Fisher score method [Duda et al., 2000] evaluates each features independently by computing the score using the Fisher criterion. The Spectral method [Zhao and Liu, 2007] for supervised feature selection. The ReliefF method [Kira and Rendell, 1992] for multi-class supervised feature selection. Robust Feature Selection (RFS) [Nie et al., 2010] selects features by solving an ℓ1,2-norm regularized regression problem. The ℓ0,2-FS method [Cai et al., 2013] which exactly selects the top k features in the supervised scenarios. The Discriminative Least Squares Regression (DLSR) [Xiang et al., 2012] which takes the ℓ1,2-norm regularized least squares regression formulation. Also, the original data with all features for classification is compared as the baseline.

The methods, Fisher, ReliefF and ℓ0,2-FS, are parameter free. The Spectral method requires the neighborhood size as a key parameter, which is tuned in the range {4, 6, 8, 10}. The regularization parameter λA for the RFS and DLSR methods is searched in the range {0.001, 0.05, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1}. The best results are reported with these parameters. To make our results reproducible, the regularization parameter γ = 0.2 is used for our method throughout the experiments.

4.3 Convergence Analysis

To solve the proposed formulation, we develop an iterative update algorithm. It is important to provide the experimental study of the convergence of the proposed method. The convergence rates on the Coil-20 data set are shown in Fig. 2. As can be seen from Fig. 2(b), most of the entries of vector v are either 0 or 1. Only a few entries lie in the range [0, 1]. This means the proposed method is effective in solving 0-1 integer programming. Fig. 2(a) presents that if a entry of v is 0, the corresponding column of matrix A would be a zero vector. We can see from Fig. 2(c) that the proposed method converges within 50 iterations, demonstrating that the proposed optimization algorithm is effective. There is a low-lying pit between the 3rd and the 10th iterations. This is because that at the early stage of the algorithm, the impact of the continuous constraints in our method is trivial, i.e., few of the entries of v are 0 or 1 or even in the range of [0, 1]. The problem is essentially equivalent to a classic least squares at the early steps. After some iterations, the two continuous constraints start exert significant influence on the objective function and then the objective value can reasonably rise.

4.4 Performance Evaluation

Given a data set, we randomly select p percents from each class to formulate the training data Xtrain and the remaining data are used as the test data. The supervised feature selection methods are performed on Xtrain to rank the features. The classifiers, both 1-NN and SVM, are trained on Xtrain represented by the selected features and then tested on the test data. For each setting, the experiments are repeated 10 times and both the mean of results and the deviation variances are reported.

Fig. 4 shows the plots of classification performance versus the number of selected features on the data sets, Coil-20, MNIST and Yale-B. The percentage of labeled training data is p = 30. As can be seen, our proposed method consistently outperforms or show comparable performance with the compared FS methods. The accuracy curves of the proposed method converge very fast, with typically around 80 features.

The filter-type FS methods (Fisher, Spectral and ReliefF) evaluate features individually and do not consider the relevance among them. Therefore, the accuracies of these methods are generally lower than those obtained by group sparsityinducing based FS methods (RFS, ℓ0,2-FS, DLSR and our method). Our proposed method can efficiently remove noise and redundant features and get comparably better performance with fewer features. By comparing with the ℓ0,2 FS method, which also selects the exact number of the top features, it is observed that our method is superior to ℓ0,2 FS on the data sets.

Due to the limited pages, the results on additional data sets and the results when the percentage of training data p = 50 are not shown here. The readers are encouraged to try other data sets with the provided Matlab code.

4.5 Parameter Sensitivity

The proposed method requires one parameter γ to be set in advance. In this subsection, we discuss the sensitivity of the proposed method over this parameter. The parameter γ is searched in the range of [0, 3]. 50 percents of data in each class are used as the training data (p=50) and the top 200 features are utilized. The results on Coil-20 and YaleB data sets are presented in Fig. 3. As can be seen, both large and small γ degrades the performance of the proposed method. On the other hand, any γ ∈ [0.2, 0.5] seems be able to give reasonable results. We set γ = 0.2 and do not tune the parameter in the experiments for different data sets under different experimental settings.

5 Conclusion

In this paper, we propose a novel supervised FS method to identify the best k features rather than the k top features. The proposed method is a classic least squares under the ℓ0 or ℓ0,2 constraint, where the non-smoothed constraint means the number of selected features. We transformed the original integer programming problem into an optimization problem with two continuous constraints. And then the problem is efficiently solved by ADMM method. It is shown that the proposed method outperforms the state-of-the-art supervised FS method on benchmark image data sets.

Acknowledgments

This work is supported by NSFC (Grants nos. 61472285, 61473212, 61503263, 61502289, 61511130084), Zhejiang Provincial Natural Science Foundation (Grants nos. LY15F030011, LY17F030004 and LR17F030001), Project of science and technology plans of Zhejiang Province (Grants no. 2015C31168). Project of science and technology plans of Wenzhou City (Grants Nos. G20150017, G20160002).