Direct Sparsity Optimization Based Feature Selection for Multi-Class Classification

1 Introduction

Feature selection has been an important component in machine learning models (Guyon & Elisseeff, 2003). Feature selection approaches in general can be divided into three groups: filter methods (Kira & Rendell, 1992; Lewis, 1992; Peng et al., 2005), wrapper methods (Guyon et al., 2002), and embedded methods (Cawley et al., 2006; Wang et al., 2008; Xiang et al., 2012). The filter methods use proxy measures that are independent on the machine learning models to rank features according to their relevance to the learning problem. The wrapper methods search subsets of features to optimize a given learning model's performance, and typically have higher computational cost than the filter methods. The embedded methods integrate the feature selection task into the model learning, and are able to achieve good performance with a moderate computational cost. Particularly, sparse linear model based feature selection methods have achieved promising performance. The sparse linear model based methods typically adopt ℓ1-norm regularization, and many variants have been proposed with different sparsity regularization terms, such as Lasso (Tibshirani, 1996) and sparse support vector machine (SVM) (Bradley & Mangasarian, 1998; Mangasarian, 2006).

In multi-task learning, various ℓ2,1-norm (Liu et al., 2009; Nie et al., 2010; Obozinski et al., 2006) and ℓ∞,1-norm (Liu et al., 2009) based regularization models have been investigated for selecting features with joint sparsity across different tasks. Moreover, group Lasso based methods (Kong & Ding, 2013; Kong et al., 2014) have also been proposed in recent years. In fact, the multi-task feature selection algorithms share similarities with group lasso based methods. Since non-convex ℓ𝑝-norm or ℓ𝑟,𝑝-norm (0 < 𝑝 < 1) based regularization models can yield sparser solutions than ℓ1 -norm or ℓ𝑟,1 -norm based models, they have gained increasing attention in recent studies (Chartrand & Staneva, 2008; Liu et al., 2007; Zhang et al., 2014).

The sparse linear model based feature selection algorithms typically take a trade-off between a data-fitting loss function term and a sparsity term, therefore there inevitably exists residual in the loss function. However, little is known about such a residual's impact on the feature selection. Most of these algorithms are developed for classification problems. However, their data-fitting loss functions adopted are typically based on the regression of class labels, rather than surrogates of 0-1 loss that are more appropriate for classification problems (Bartlett et al., 2004). To overcome the limitations of existing sparse model based feature selection methods, we present a novel feature selection method via directly optimizing a linear model's sparsity with ℓ2,𝑝 -norm (0 < 𝑝 ≤ 1) , subject to data-fitting inequality constraints, instead of adopting the sparsity as a regularization term. Our constrained optimization formulation circumvents the difficulty of regularization parameter setting, and the inequality constraint based data-fitting loss function enables large between-class margins, similar to SVMs, but capable of handling multi-class problems. We propose an efficient algorithm to solve the optimization problem associated with the direct sparsity optimization, which is non-convex and non-smooth when 0 < 𝑝 < 1, by transforming it to a Frobenius-norm induced problem at each iteration step, which has been proved to converge to Karush-Kuhn-Tucker (KKT) points.

The proposed algorithm has been evaluated based on 6 publicly available datasets, and extensive comparison experiments have demonstrated that our algorithm could achieve feature selection performance competitive to state-of-the-art algorithms, including ReliefF (Kira & Rendell, 1992), Min-Redundancy Max-Relevance (mRMR) (Peng, et al., 2005), ℓ1 -norm based ℓ1 -SVM (Bradley & Mangasarian, 1998; Mangasarian, 2006), and ℓ2,1 -norm based Robust Feature Selection (RFS) (Nie, et al., 2010; Xiang, et al., 2012).

2 Direct Sparsity Optimization Based Feature Selection (DSO-FS)

Throughout this paper, matrices are written in bold uppercase, vectors are written in bold lowercase, and all the scalars are denoted by normal letters. 𝑰 denotes an identity matrix and 𝟏 denotes a vector or matrix with all the elements equal to 1. Given a matrix 𝑨 ∈ ℝ𝑚×𝑛 , the ℓ𝑟,𝑝-norm(𝑟 > 0, 𝑝 > 0) 1 of 𝑨 is defined as:

where ‖𝒂𝑖‖𝑟 denotes ℓ𝑟-norm of the 𝑖-th row vector of 𝑨.

Given 𝑚 training samples {𝒙𝑖 , 𝑦𝑖 }𝑖=1 𝑚 where 𝒙𝑖 ∈ ℝ𝑛 is a data point and 𝑦𝑖 is its associated class label in 𝑐(𝑐 ≥ 2) classes, the multiclass classification problem can be modeled as a linear learning problem. For simplicity, the bias of the standard linear regression is absorbed into 𝑾 as an additional dimension with all elements equal to 1.

where 𝑿 = [𝒙1 ; 𝒙2 ;…; 𝒙𝑚; 1] , 𝑾 ∈ ℝ(𝑛+1)×𝑐 is the weight matrix to be learned, and 𝒀 = [𝒇1 ; 𝒇2 ;… , 𝒇𝑖 ;…; 𝒇𝑚] ∈ ℝ𝑚×𝑐 is a class label matrix with labels rearranged using a one-versus-rest model, i.e., 𝒇𝑖 = [-1, …,1, … , -1] ∈ ℝc (the 𝑗-th element is 1 and others

For selecting a subset of features, 𝑾 should have sparse columns and share joint sparsity along its row direction since each row of 𝑾 corresponds to the same feature. Therefore, we model the feature selection problem as

where ‖𝐖‖2,0 is the number of non-zero rows in 𝑾, of which not all the elements are zero.

It is NP-hard to solve the optimization problem of Eqn. (3). Therefore ℓ2,𝑝-norm (0 < 𝑝 ≤ 1) can be adopted instead, resulting in a relaxed sparsity optimization problem:

1 If 0 < 𝑟 < 1 or 0 < 𝑝 < 1, ℓ𝑟,𝑝-norm does not satisfy triangle inequality. However, this does not affect the proposed algorithm.

where 0 < 𝑝 ≤ 1.

Many studies (Liu, et al., 2009; Nie, et al., 2010; Obozinski, et al., 2006) assumed that the ℓ2,1-norm based problems is equivalent to ℓ2,0-norm based problems under certain conditions (Candes & Tao, 2005). In practice, ℓ2,𝑝-norm (0 < 𝑝 < 1) can lead to sparser solutions in most cases although it is non-convex (Chartrand & Staneva, 2008; Liu, et al., 2007; Zhang, et al., 2014).

It is desired that the classification model's margins between classes are as large as possible for obtaining improved generalization performance. Accordingly, the equality constraint in Eqn. (4) is relaxed to be inequality constraints, i.e., 𝑚𝑖𝑛

where 𝟏 ∈ ℝm×c , ⊙ is a Hadamard product operator for element-wise multiplication between two matrices of the same dimensions, and ≽ denotes that the elements of the left matrix are greater than or equal to their corresponding ones of the right matrix.

The optimization problem of Eqn. (5) can be reformulated by introducing a slack variable 𝑬 ∈ ℝ𝑚×𝑐 whose elements have the same positive or negative sign as the corresponding elements of 𝒀, as

where 𝑬 ∈ ℝ𝑚×𝑐 .

To solve the optimization problem of Eqn. (6), we first solve the linear equation 𝑿𝑾= 𝒀+𝑬 to obtain the solution space of 𝑾, and then directly search the solution space to find a solution with the minimum of ‖𝑾‖2,𝑝. Note that if 𝑿𝑾= 𝒀+𝑬 is inconsistent, especially when the number of data samples 𝑚 is greater than the number of features 𝑛, least-square solution space of the equation can be used as a substitute. Actually, we only need to solve 𝑿𝑾= 𝑿𝑿+(𝒀+ 𝑬), where 𝑿+ is pseudo-inverse. This equation is compatible when 𝑿 is row full rank, since 𝑿𝑿+ = 𝑰 on this occasion. Gaussian Elimination is a simple and efficient way to obtain the solution space of 𝑾. Without loss the generality, we assume that the first 𝑚0 column vectors of 𝑿 are linearly independent, i.e.

where the rank of 𝑿 is 𝑚0 , 𝑿1 ∈ ℝ𝑚×𝑚0 , 𝑿2 ∈ ℝ𝑚×𝑛0 , 𝑰 ∈ ℝ𝑚0×𝑚0 , 𝑴 ∈ ℝ𝑚0×𝑛0 , 𝑵 ∈ ℝ𝑚0×𝑐 , and 𝑫 ∈ ℝ𝑚0×𝑚0 is product matrix of a series of elemental matrices, 𝑳 = 𝑫𝑿𝑿+ ∈ ℝ𝑚0×𝑚0 , and 𝑛0 = 𝑛 + 1-𝑚0. Thus, the solution space of 𝑾 is

where 𝑷 ∈ ℝ𝑛×𝑛0 , 𝑸 ∈ ℝ𝑛×𝑐 , 𝑭 ∈ ℝ𝑛×𝑐 , 𝑰 ∈ ℝ𝑛0×𝑛0 , and 𝑼 ∈ ℝ𝑛0×𝑐 is an arbitrary matrix.

Finally, the optimization problem of Eqn. (5) can be reformulated as

where 0 < 𝑝 ≤ 1.

3 An Iterative Algorithm for DSO-FS

We propose an iterative algorithm to solve the optimization problem of Eqn. (9) due to that no analytical solution is available. At each iteration step, we alternately optimize variables 𝑼 and 𝑬. An objective function with ℓ2,p -norm (0 < 𝑝 ≤ 1) is non-smooth and non-convex when 0 < 𝑝 < 1. To efficiently solve this problem, ℓ2,p-norm is reformulated by Frobenius-norm (ℱ-norm) that is smooth and convex, as

where 𝜮 is a diagonal matrix with the 𝑖-th diagonal element 𝜮𝑖𝑖 = 1/‖𝒂𝑖‖2 1-𝑝 2 / , and ‖𝒂i‖2 is defined in Eqn.(1) 2 . Note that 𝒁∗ = argmin 𝒁‖𝑨(𝒁)‖2,𝑝 𝑝 = argmin 𝒁‖𝑨(𝒁)‖2,𝑝, where 𝑨(𝒁) is a function of 𝒁. The similar strategies has also been adopted in FOcal Underdetermined System Solver (FOCUSS) (Cotter et al., 2005; Gorodnitsky & Rao, 1997), and various optimization techniques for ℓ2,1-norm regularization based feature selection algorithms (Hou et al., 2011; Krishnapuram et al., 2005; Merchante et al., 2012; Nie, et al., 2010; Xiang, et al., 2012; Yi Yang et al., 2011) can be considered as its special case when 𝑝 = 1. Our solution is summarized in Algorithm 1.

Algorithm 1. Feature Selection via Direct Sparsity Optimization (DSO-FS)

Input: data points {𝒙𝑖 }𝑖=1 𝑚 (𝒙𝑖 ∈ ℝ𝑛 ) and their corresponding label {𝑦𝑖 }𝑖=1 𝑚 ; norm power 𝑝; number of features 𝑑 to be selected.

Construct 𝑿 and 𝒀 according to Eqn. (2), and 𝑳, 𝑴, 𝑵, 𝐏 and 𝑸 according to Eqn. (7) and Eqn. (8).

Set 𝑘 = 0 and initialize 𝑬0 and 𝜮0 ∈ ℝ𝑛×𝑛 with a zero matrix and an identity matrix respectively.

2 𝜮 = 𝚷-𝟏 , where 𝚷 is a diagonal matrix with the 𝑖-th diagonal element 𝚷𝑖𝑖 = ‖𝒂𝒊‖2 1-𝑝 2 / . If any row vector 𝒂𝑖 = 0, we define 𝜮 = 𝚷+.

Update 𝜮𝑘 , where the 𝑖-th diagonal element is 𝟏 ‖𝒘𝑖𝑘‖21-𝑝 2 /

until convergence

Sort all features according to ‖𝒘𝑖 ‖2 and select the largest 𝑑 features.

Note that at each iteration step, for solving 𝑼𝑘+1 = argmin 𝑼‖𝜮𝑘 (𝑷𝑼+𝑮)‖𝐹2 , an analytical solution is available, i.e.,

where 𝑺 = (𝜮𝑘 )2 .

According to Eqn. (8), Eqn. (11) can be reformulated as

where 𝑺1 ∈ ℝ𝑚0×𝑚0 is a diagonal matrix, of which the diagonal elements are the first m0 elements of 𝑺 , 𝑺𝟐 ∈ ℝ𝑛0×𝑛0 is a diagonal matrix, of which the diagonal elements are the last 𝑛0 elements of 𝑺, and 𝑲 = 𝑵+𝑳𝑬𝑘 (∈ ℝ𝑚0×𝑐 ).

If 𝑚0 < 𝑛0 , the optimization problem's computational cost can be further reduced by a simple matrix operation

where 𝑻 = 𝑺2-𝟏𝑴𝑡 𝑺1 and 𝑰 ∈ ℝ𝑚0×𝑚0 .

For feature selection problems with the number of samples less than the number of features, i.e., 𝑚0 ≤ 𝑚 ≪ 𝑛0 , updating Eqn. (13) has reduced computational cost than Eqn. (11). It is worth noting that it is not necessary to directly calculate 𝑪 = (𝑴𝑻+ 𝑰)-𝟏 𝑲 , which is computationally more expensive than solving the following linear equation

Although no analytical solution is available for 𝑬𝑘+1 = argmin 𝑬‖𝜦𝒌 (𝑳𝑬+𝑯)‖𝐹2 (s. t. 𝒀⊙𝑬 ≽ 𝟎), it is a smooth and convex optimization problem. Therefore, it can be efficiently solved by existing tools, such as CVX (CVX Research, 2011).

4 Convergence proof

Algorithm 1 makes ‖𝑾‖2,𝑝 to monotonically decrease at every iteration step and finally converges.

Lemma 1. Given any two vectors 𝒂 and 𝒃, we have

where 0 < 𝜃 < 1 and the equality holds if and only if 𝑎 = 𝑏. Proof. Since ln(𝑥 2) is concave, we have

where 0 < 𝜃 < 1. The equality holds if and only if 𝑥1 = 𝑥2, indicating that (1 -𝜃)𝑥12 +𝜃𝑥22 ≥ 𝑥12-2𝜃 𝑥22𝜃 . (17)

Then, we have

where the equality holds if and only if 𝒂 = 𝒃. □

Lemma 2. Given an optimization problem:

where 𝜱(𝒁) is a function of 𝒁, 𝓕 is the feasible region, and 𝑺 is a diagonal matrix whose 𝑖 -th diagonal element is 1/‖𝜱(𝒁0)𝑖‖2 1-𝑝 2 / (𝒁0 could be any object in 𝓕, 𝜱(𝒁0)𝑖 is the 𝑖-th row vector of 𝜱(𝒁0) and 0 < 𝑝 ≤ 2 ), we have

where 𝒁∗ is the optimal solution of Eqn. (19) and the equality holds if and only if 𝜱(𝒁∗) = 𝜱(𝒁0)

Proof. Since 𝒁∗ is the optimal solution, we have

Then

where 𝜱(𝒁𝟎)𝑖 and 𝜱(𝒁∗ )𝑖 are the 𝑖-th row vector of 𝜱(𝒁0) and 𝜱(𝒁∗ ), respectively.

According to Lemma 1, we have

Then dividing the both sides by ‖𝜱(𝒁0)𝑖‖22-𝑝 , we have

It indicates that

Combining Eqns. (21) and (24), we obtain

Therefore,

where the equality holds if and only if 𝜱(𝒁∗ ) = 𝜱(𝒁𝟎).□

Theorem 1. The sequence {𝑾𝑘 } produced via Algorithm 1 has the following properties: ‖𝑾𝑘 ‖2,𝑝 is non-increasing at successive iteration steps and {‖𝑾𝑘 ‖2,𝑝} converges to a limited value.

Proof. Supposing we have obtained the solution 𝑼𝑘 , 𝑬𝑘 , and the objective function 𝑾𝑘 at the (𝑘 + 1)-th iteration step, we solve the optimization problem min 𝑼‖𝜮𝒌(𝑷𝑼+𝑮)‖𝑭2 to obtain 𝑼𝑘+1 by fixing 𝑬𝑘 . According to Lemma 2, we have

Then we fix 𝑼𝑘+1 ， and solve the optimization problem min 𝑬‖𝜦𝒌(𝑳𝑬+𝑯)‖𝑭2 to obtain 𝑬𝑘+1 . According to Lemma

Combining (27) and (28), we obtain

where the equality holds if and only if 𝐖k+1 = 𝐕𝐤 =𝐖k . Since the lower bound of ‖𝑾‖2,p is limited, {‖𝑾𝑘 ‖2,𝑝} will converge. □

Theorem 2. If sequences {𝑾𝑘 } and {𝑬𝑘 } produced in Algorithm 1 have limit points, the limit points satisfy the Karush-Kuhn-Tucker (KKT) conditions of Eqn. (6). When 𝑝 ≥ 1, the limited points are globally optimal.

Due to space limitations, the proof is provided in supplementary document (Peng & Fan, 2016)..

It is worth noting that Lemma 2 works for sparsity regularization based feature selection algorithms, i.e.,

where λ is the regularization coefficient.

Eqn. (30) can be reformulated as

The Robust Feature Selection (RFS) (Nie, et al., 2010; Xiang, et al., 2012) is a special case of Eqn. (30) when 𝑖 = 1 and 𝑝 = 1.

5 Experiments

5.1 Experimental datasets and settings

The proposed algorithm has been evaluated based on 6 publicly available datasets. In particular, 2 datasets were obtained from UCI, including ISOLET and SEMEION. ISOLET is a speech recognition data set with 7797 samples from 26 classes, and each sample has 617 features; SEMEION contains 1593 handwritten images from ~80 people, stretched in a rectangular box of 16× 16. For these 2 datasets, the number of features is less than the number of their samples. Another 2 datasets were microarray data, including LUNG and CLL-SUB-111. In particular, LUNG consists of 203 samples, each having 3312 genes with standard deviations larger than 50 expression units (Bhattacharjee et al., 2001); CLL-SUB-111 consists of 111 samples with 11340 features (featureselection.asu.edu). Our algorithm has also been validated based on 2 image datasets, including UMIST and AR. In particular, UMIST includes face images with a resolution of 56× 46 from 20 different people, and AR has 130 samples with 2400 features.

We compared our method with sparsity regularization based feature selection methods, including ℓ1-SVM(Bradley & Mangasarian, 1998; Mangasarian, 2006) and ℓ2,1-norm based Robust Feature Selection (RFS) (Nie, et al., 2010; Xiang, et al., 2012), with an objective function 𝒥(𝑾) = ‖𝑿𝑾-𝒀‖2,1 + 𝜆‖𝑾‖2,1. We also compared our algorithm with well-known filter methods, including ReliefF (Kira & Rendell, 1992) and mRMR (Peng, et al., 2005).

The feature selection methods were evaluated based on their classification accuracy. Particularly, linear SVM (Chang & Lin, 2011) was adopted to build classifiers based on the selected features. The parameter 𝐶 of linear SVM classifiers were tuned using a cross-validation strategy by searching a candidate set [10 -4 , 10 -3 , 10 -2 , 10 -1 , 1, 10 1 , 10 2 ]. The regularization parameter of ℓ1 -SVM and RFS were tuned using the same cross-validation strategy by searching a candidate set [10 -3 , 10 -2 , 10 -1 , 1, 10 1 , 10 2 , 10 3 ].

In our experiments, we first normalized all the data to have

0 mean and unit standard deviation for each feature. 10 trials were performed on each dataset. In each trial, the samples of each dataset were randomly spitted into training and testing subsets with a ratio of 6:4. For tuning parameters, a 3-fold was used for datasets with less than 200 training samples, and an 8-fold cross-validation was used for other datasets.

5.2 Effect of parameter 𝒑

To investigate how the classification performance is affected by the parameter 𝑝 of DSO-FS, we performed experiments based on datasets ISOLET, LUNG, and UMIST. We obtained solutions of 𝑾 with different settings of 𝑝 ∈ [0.1, 0.3, 0.5, 0.7, 0.9, 1] using Algorithm 1 based on the above 3 datasets, and then selected top ranked features according to ℓ2 -norm ‖𝒘𝑖‖2 to build classifiers. Figure 1 shows the classification accuracy as a function of the number of the selected features and the parameter 𝑝. The results shown in Figure 1 indicated that 𝑝 played an important role in the classification.

5.3 Comparisons with state-of-the-art methods

Since the classification performance of linear SVM classifiers built on the features selected by our method is hinged on the parameter 𝑝 , we used a cross-validation strategy to select an optimal value from [0.1, 0.3, 0.5, 0.7, 0.9, 1] for our method. Figure 2 shows the average classification performance of classifiers built on the features selected by different methods in 10 trials. In particular, the average classification accuracy is shown as a function of the number of features used in the classification model. Compared with other methods, the proposed method achieved higher classification accuracy on most datasets, indicating that our method had overall better performance than other algorithms.

Table 1 summarizes mean and standard deviation of the classification rates in 10 trails for classifiers built on the top 100 features selected by the algorithms under comparison. These results demonstrated that our algorithm had the overall best classification accuracy on all the 6 datasets.

6 Conclusions and Discussions

In this paper, a novel feature selection algorithm via direct sparsity optimization was proposed. Our method directly optimizes a large margin linear classification model's sparsity with ℓ2,𝑝-norm (0 < 𝑝 ≤ 1) subject to a data-fitting constraint, different from the sparse regularization based algorithms that typically take a trade-off between a data-fitting loss function term and a sparsity term. Since little is known about such a trade-off's impact on the feature selection, empirical parameter tuning is typically adopted for choosing the trade-off parameter. Our feature selection method adopts data-fitting inequality constraints to obtain large between-class margins, which could improve the generalization ability of the selected features, rather than data-fitting loss functions built upon the regression of class labels (Bartlett, et al., 2004). We also proposed an efficient algorithm to solve the non-convex

( 0 < 𝑝 < 1 ) and non-smooth optimization problem associated with the feature selection problem. Extensive experiments based on 6 datasets have demonstrated that the proposed method could achieve better performance than state-of-the-art feature selection algorithms. Furthermore, our algorithm can be easily extended for solving other sparsity regularization algorithms. In particular, our algorithm could be used to solve ℓ0 and ℓ2,0 based optimization problems subject to linear constraints by setting 𝑝 close to 0. Furthermore, similar to SVM, we could also add soft-margin to our algorithm for better performance in the future work.

Acknowledgments

This work was supported in part by National Key Basic Research and Development Program of China (No. 2015CB856404), NSFC (No. 81271514, 61473296), and NIH grants MH070365 and AG014971.