A Semisupervised Feature Selection with Support Vector Machine

Abstract

Feature selection has proved to be a beneficial tool in learning problems with the main advantages of interpretation and generalization. Most existing feature selection methods do not achieve optimal classification performance, since they neglect the correlations among highly correlated features which all contribute to classification. In this paper, a novel semisupervised feature selection algorithm based on support vector machine (SVM) is proposed, termed SENFS. In order to solve SENFS, an efficient algorithm based on the alternating direction method of multipliers is then developed. One advantage of SENFS is that it encourages highly correlated features to be selected or removed together. Experimental results demonstrate the effectiveness of our feature selection method on simulation data and benchmark data sets.

1. Introduction

Feature selection, with the purpose of selecting relevant feature subsets among thousands of potentially irrelevant and redundant features, is a challenging topic of pattern recognition research that has attracted much attention over the last few years. A good feature selection method has several advantages for a learning algorithm such as reducing computational cost, increasing its classification accuracy, and improving result comprehensibility [1].

Considering the usage of the class label information, feature selection methods can be classified into supervised methods, unsupervised methods, and semisupervised methods. Supervised feature selection methods usually use only information from labeled data to find the relevant feature subsets [2–4]. However, in many real, world applications, the labeled data are very expensive or difficult to obtain, which brings difficulty to create a large training data set. This situation arises naturally in practice, where large amount of data can be collected automatically and cheaply, when manual labeling of samples remains difficult, expensive and time consuming. Unsupervised feature selection methods could be an alternative in this case through exploiting the information conveyed by the large amount of unlabeled data [5, 6]. However, as these unsupervised algorithms ignore label information, important hints from labeled data are left out and this will generally downgrade the performance of unsupervised feature selection algorithms. The combination of both supervised methods and unsupervised methods is semisupervised approaches [7–10] which exploit the information of both labeled and unlabeled data. A good survey about semisupervised feature selection approaches can be found in [9].

The performances of the most existing semisupervised feature selection methods are insufficient when there are several highly correlated features, which are all relevant to classification and the way they interact can help with the interpretability of the objective problem [11]. Given these premises, this paper provides two main contributions as follows.(i)We present a novel semisupervised feature selection scheme based on support vector machine (SVM) and the elastic net penalty proposed by Zou and Hastie [12] combining 𝑙₁ and 𝑙₂ regularizations, termed SENFS.(ii)In order to solve SENFS with the nondifferentiability of both the loss function and the 𝑙₁-norm regularization term, an efficient algorithm based on the alternating direction method of multipliers (ADMM) [13] is developed.

Compared with other semisupervised feature selection algorithms, SENFS provides the following benefits.(i)It permits highly correlated features to be selected or removed together.(ii)It performs automatic feature selection as part of the training process and can achieve better classification performance using the selected features.

The effectiveness of SENFS is validated on simulated data and six benchmark semisupervised data sets. Our main finding is that SENFS can identify the features that are relevant to classification using the data set that consists of only a few labeled samples and many unlabeled samples.

This paper is organized as follows. Section 2 briefly introduces the methodology. In Section 3, we derive an iterative algorithm that yields the entire solution path based on ADMM to solve this proposed method. In Section 4, we evaluate the performance of this proposed method on both simulated and real-world data, followed by a summary in Section 5.

2. Methodology

Assume that all samples sampled from the same population generated by target concept consist of 𝑝 features. Given a set of samples 𝐗=(𝑥₁,…,𝑥_𝑛)^𝑇, in which 𝑛 is the number of samples, the 𝑖th sample or input vector 𝑥_𝑖 of original feature 𝐷 with 𝑝 features is denoted by 𝑥_𝑖={𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑝}. The set 𝐗 can be divided into two parts: labeled set 𝐗_𝑙=(𝑥₁,…,𝑥_𝑙)^𝑇 for which labels 𝐲_𝑙=(𝑦₁,…,𝑦_𝑙)^𝑇 are provided with 𝑦_𝑖∈{−1,1}^𝑙_𝑖=1 for binary problem and unlabeled set 𝐗_𝑢=(𝑥_𝑙+1,…,𝑥_𝑙+𝑢)^𝑇 whose labels are not given, where 𝑙 and 𝑢 are the number of labeled and unlabeled samples, respectively, and 𝑛=𝑙+𝑢. Then, the generic goal of semisupervised feature selection is to find a feature subset 𝐷^′ with 𝑑(𝑑<𝑝) features which contains the most informative features using both data information of 𝐗_𝑙 and 𝐗_𝑢. In other words, the samples (𝑥^′₁,…,𝑥^′_𝑛)^𝑇 represented in the 𝑑-dimensional space can well preserve the information of the samples 𝐗=(𝑥₁,…,𝑥_𝑛)^𝑇 represented in the original 𝑝-dimensional space.

We begin our discussion with the binary supervised feature selection based on the elastic net penalty. Wang et al. [11] proposed a supervised feature selection method based on SVM with the elastic net penalty term named doubly regularized support vector machine for binary classification problems, which solves the optimization of the following generic objective function over both the hyper plane parameters (𝛽,𝛽₀):min_𝛽,𝛽0𝐶^𝑛_𝑖=1𝑉𝑦_𝑖𝑔𝑥_𝑖+𝜆₁𝛽₁+𝜆₂𝛽²₂2, (1) where the decision function is defined as 𝑔(𝑥)=𝑥𝛽+𝛽₀ and both 𝜆₁ and 𝜆₂ are tuning parameters, and 𝐶 is the regularization parameter. 𝑉 is the margin loss function; for example, hinge loss 𝑉(𝑧)=𝐻₁(𝑧)=max(1−𝑧,0). The role of the 𝑙₁-norm penalty is to allow selection, and the role of the 𝑙₂-norm penalty is to help groups of highly correlated features get selected or removed together which is denoted by the grouping effect [11].

As for semisupervised feature selection, considering 𝐗_𝑙 and 𝐗_𝑢, inspired by the semisupervised learning algorithm TSVM [14], we apply the elastic net penalty for semisupervised feature selection (SENFS), which solves the following optimization task over both the hyper plane parameters (𝛽,𝛽₀) and the unlabeled vector 𝐲_𝑢=(𝑦_𝑙+1,…,𝑦_𝑛)^𝑇: min_{𝛽,𝛽0,𝐲𝑢∈{}−1,+1}𝐶^𝑙_𝑖=1𝑉𝑦_𝑖𝑔𝑥_𝑖+𝐶^∗𝑛_𝑖=𝑙+1𝑈𝑦_𝑖𝑔𝑥_𝑖+𝜆₁𝛽₁+𝜆₂𝛽²₂21s.t.𝑢^𝑛_𝑖=𝑙+1max0,𝑦_𝑖=𝑐. (2)

The constraint in (2) is called the balancing constraint and is necessary to avoid the trivial solutions where all unlabeled samples are assigned to the same class. This constraint enforces a manually chosen constant 𝑐 and an approximation of this constraint writes (1/𝑛)∑^𝑛_𝑖=1𝑔(𝑥_𝑖)=2𝑐−1 [15] with 𝑐. So the constraint can be rewritten as (1/𝑢)∑^𝑛_𝑖=𝑙+1𝑔(𝑥_𝑖)=(1/𝑙)∑^𝑙_𝑖=1𝑦_𝑖. 𝑉 and 𝑈 employ the same loss; for example, hinge loss 𝑉(𝑧)=𝐻₁(𝑧)=max(1−𝑧,0).

Obviously, the difficulty of the above optimization task consists in finding the optimal assignment for the unlabeled vector 𝑦 and the hyper plane parameters (𝛽,𝛽₀), which is a mixed-integer programming problem [16]. As described in [15], for a fixed (𝛽,𝛽₀), argmin_𝑦𝑉(𝑦𝑔(𝑥))=sign(𝑔(𝑥)). So the problem of (2) can be seen equivalently as min_𝛽,𝛽0𝐶^𝑙_𝑖=1𝑉𝑦_𝑖𝑔𝑥_𝑖+𝐶^∗𝑛_𝑖=𝑙+1𝑈𝑔𝑥_𝑖+𝜆₁𝛽₁+𝜆₂𝛽²₂21s.t.𝑢^𝑛_𝑖=𝑙+1𝑔𝑥_𝑖=1𝑙^𝑙_𝑖=1𝑦_𝑖. (3)

On the other hand, one effective approximation of the loss function 𝑈(|𝑧|) was a clipped variant [17] which can be expressed as 𝑈(|𝑧|)=𝑅_𝑠(𝑧)+𝑅_𝑠(−𝑧)−(1−𝑠), (4) where 𝑅_𝑠(𝑧) is the Ramp loss defined as 𝑅_𝑠(𝑧)=𝐻₁(𝑧)−𝐻_𝑠(𝑧) with 𝐻_𝑠(𝑧)=max(𝑠−𝑧,0),0≤𝑠<1. In our experiments, the typical value of 𝑠 is 0.3. The main reason to use the clipped symmetric hinge loss is the gain of sparsity in the number of support vectors yielded by the optimizer [17]. Then we can get min_𝛽,𝛽0𝐶^∗𝑛_𝑖=𝑙+1𝑈𝑔𝑥_𝑖=min_𝛽,𝛽0𝐶^∗𝑛_𝑖=𝑙+1𝑅_𝑠(𝑧)+𝑅_𝑠(−𝑧)−(1−𝑠)=min_𝛽,𝛽0𝐶^∗𝑛_𝑖=𝑙+1𝐻₁𝑔𝑥_𝑖+𝐻₁−𝑔𝑥_𝑖−𝐶^∗𝑛_𝑖=𝑙+1𝐻_𝑠𝑔𝑥_𝑖+𝐻_𝑠−𝑔𝑥_𝑖. (5)

From (5), we can know that solving the optimization problem (3) with the clipped symmetric hinge loss is equivalent to solving a classical SVM with the unlabeled samples counted twice with 𝑦_𝑖=1 when 𝑙+1≤𝑖≤𝑙+𝑢 and 𝑦_𝑖=−1 when 𝑙+𝑢+1≤𝑖≤𝑙+2𝑢, which are artificial labels. Therefore problem (3) can be rewritten as min_𝛽,𝛽0𝐶^𝑙_𝑖=1𝑉𝑦_𝑖𝑔𝑥_𝑖+𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑔𝑥_𝑖−𝐻_𝑠𝑦_𝑖𝑔𝑥_𝑖+𝜆₁𝛽₁+𝜆₂𝛽²₂21s.t.𝑢^𝑛_𝑖=𝑙+1𝑔𝑥_𝑖=1𝑙^𝑙_𝑖=1𝑦_𝑖. (6)

As we can seen from (6), when 𝐶≠0,𝐶^∗=0, SENFS evolves into a supervised feature selection algorithm, and when 𝐶=0,𝐶^∗≠0, it becomes an unsupervised model.

In the following, we will illustrate how SENFS has the grouping effect for correlated features. The following theorem describes this point.

Theorem 1. Denote the solution to (6) by (𝛽,𝛽₀) and 𝛽₀, the input 𝑗th feature by 𝐹_𝑗, and the input 𝑘th feature by 𝐹_𝑘. Then for any pair (𝑗,𝑘), one can have 𝛽_𝑗−𝛽_𝑘≤2𝑀𝜆₂^𝑙_𝑖=1𝑥_𝑖𝑗−𝑥_𝑖𝑘+2𝑀^′𝜆₂^𝑙+2𝑢_𝑖=𝑙+1𝑥_𝑖𝑗−𝑥_𝑖𝑘, (7) where 𝑀 and 𝑀^′are positive finite constants. Furthermore, if the input features 𝐹_𝑗 and 𝐹_𝑘 are centered and normalized, then 𝛽_𝑗−𝛽_𝑘≤2𝜆₂2𝑙1−𝜌_𝑙+2𝑢1−𝜌_𝑢, (8) where 𝜌_𝑙 and 𝜌_𝑢 are the sample correlations between 𝐹_𝑗 and 𝐹_𝑘,𝜌_𝑙=∑^𝑙_𝑖=1(𝑥_𝑖𝑗−𝑥_.𝑗)(𝑥_𝑖𝑘−𝑥_.𝑘)/∑^𝑙_𝑖=1(𝑥_𝑖𝑗−𝑥_.𝑗)²∑^𝑙_𝑖=1(𝑥_𝑖𝑘−𝑥_.𝑘)², 𝑥_.𝑗=(1/𝑙)∑^𝑙_𝑖=1𝑥_𝑖𝑗, 𝑥_.𝑘=(1/𝑙)∑^𝑙_𝑖=1𝑥_𝑖𝑘, and 𝜌_𝑘=∑^𝑙+2𝑢_𝑖=𝑙+1(𝑥_𝑖𝑗−𝑥^′_.𝑗)(𝑥_𝑖𝑘−𝑥^′_.𝑘)/∑^𝑙+2𝑢_𝑖=𝑙+1(𝑥_𝑖𝑗−𝑥^′_.𝑗)²∑^𝑙+2𝑢_𝑖=𝑙+1(𝑥_𝑖𝑘−𝑥^′_.𝑘)², 𝑥^′_.𝑗=(1/2𝑢)∑^𝑙+2𝑢_𝑖=𝑙+1𝑥_𝑖𝑗,𝑥^′_.𝑘=(1/2𝑢)∑^𝑙+2𝑢_𝑖=𝑙+1𝑥_𝑖𝑘.

The term |𝛽_𝑗−𝛽_𝑘| in (7) describes the difference between the coefficient paths of 𝐹_𝑗 and 𝐹_𝑙. If both features are highly correlated, that is, 𝜌=1, Theorem 1 says that difference between the coefficient paths of them is almost 0, in which case both features will be selected or removed together. The upper bound in (7) or (8) provides a quantitative description for the grouping effect of SENFS.

3. Algorithm for SENFS

The alternating direction method of multipliers (ADMM) developed in the 1970s and is well suited to distributed convex optimization and in particular to large-scale problems arising in statistics, machine learning, and related areas. The method is closely related to many other algorithms, such as the method of multipliers [18], Douglas-Rachford splitting [19], Bregman iterative algorithms [20] for 𝑙₁ problems, and others.

In this section, we first propose an efficient algorithm to solve SENFS based on ADMM by introducing auxiliary variables and reformulating the original problem. Then prove its convergence property and get the adjustment principle for penalty parameters. Finally describe the stopping criterion and computational cost.

3.1. Deriving ADMM for SENFS

It is hard to solve the model (6) directly due to the nondifferentiability of three loss functions and a 𝑙₁-norm term. In order to derive an ADMM algorithm, we introduce some auxiliary variables to handle these nondifferentiable terms.

Let 𝑋_𝑙={(𝑥_𝑖𝑗)^𝑙,𝑝_{𝑖=1,𝑗=1}} denote labeled data, let 𝑋_𝑢={(𝑥_𝑖𝑗)^{𝑙+2𝑢,𝑝}_{𝑖=𝑙+1,𝑗=1}} denote unlabeled data, and let 𝑌_𝑙,𝑌_𝑢 be diagonal matrixes with their diagonal elements to be the vector 𝐲_𝑙=(𝑦₁,…,𝑦_𝑙)^𝑇 and 𝐲_𝑢=(𝑦_𝑙+1,…,𝑦_𝑙+2𝑢)^𝑇, respectively. The constrained problem in (6) can be reformulated into an equivalent form min_𝛽,𝛽0𝐶^𝑙_𝑖=1ℎ_𝑖++𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝑎_𝑖+−𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝑏_𝑖++𝜆₁‖𝐭‖₁+𝜆₂2𝛽²₂s.t.𝐡=𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙𝐚=𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢𝐛=𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢1𝐭=𝛽𝑢^𝑛_𝑖=∑_𝑙+1𝑥_𝑖𝛽+𝛽₀=𝑟, (9) where 𝐡=(ℎ₁,…,ℎ_𝑙)^𝑇, 𝐚=(𝑎₁,…,𝑎_2𝑢)^𝑇, and 𝐛=(𝑏₁,…,𝑏_2𝑢)^𝑇, 𝟏_𝑙 is an 𝑙-column vector of 1s and 𝟏_𝑢 is an 2𝑢-column vector of 1s, and 𝑟=(1/𝑙)∑^𝑙_𝑖=1𝑦_𝑖. The Lagrangian function of (9) is 𝑙𝛾=𝐶^𝑙_𝑖=1ℎ_𝑖++𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝑎_𝑖+−𝑏_𝑖++𝜆₁‖𝐭‖₁+𝜆₂2𝛽²₂+𝐮^𝑇_ℎ𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙−𝐡+𝐮^𝑇_𝑎𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢−𝐚+𝐮^𝑇_𝑏𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢−𝐛+𝐯^𝑇1𝛽−𝐭+𝑞𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽+𝛽₀. (10)

In problem (10), 𝛾={𝛽,𝛽₀,𝐡,𝐚,𝐛,𝐭,𝐮_ℎ,𝐮_𝑎,𝐮_𝑏,𝐯,𝑞}, and 𝐮_ℎ, 𝐮_𝑎, and 𝐮_𝑏 are dual variables corresponding to the constraints 𝐡=𝟏_𝑙−𝑌_𝑙(𝑋_𝑙𝛽+𝛽₀𝟏_𝑙), 𝐚=𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢), and 𝐛=𝑠𝟏_𝑢−𝑌(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢), respectively, 𝐯 is corresponding to the constraint 𝐭=𝛽, and 𝑞 is a scalar corresponding to the balancing constrain. As in the method of multipliers, we form the augment Lagrangian 𝐿𝛾𝛾+𝜇=𝑙¹2𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙²₂+𝜇₂2𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢²₂+𝜇₃2𝑠𝟏_𝑢−𝑌𝑋_𝑢𝛽+𝛽₀𝟏_𝑢²₂+𝜇₄2𝛽−𝐭²₂+𝜇₅21𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽+𝛽₀²₂, (11) where 𝜇₁,𝜇₂,𝜇₃,𝜇₄>0 are parameters. Problem (11) is the form of ADMM, which consists of the following iterations: 𝛽^𝑘+1,𝛽₀^𝑘+1,𝐡^𝑘+1,𝐚^𝑘+1,𝐛^𝑘+1,𝐭^𝑘+1=argmin_{𝛽,𝛽0,𝐡,𝐚,𝐛,𝐭}𝐿𝛽,𝛽₀,𝐡,𝐚,𝐛,𝐭,𝐮^𝑘_ℎ,𝐮^𝑘_𝑎,𝐮^𝑘_𝑏,𝐯^𝑘,𝑞^𝑘,𝐮_ℎ^𝑘+1=𝐮^𝑘_ℎ+𝜇₁𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘+1,𝐮_𝑎^𝑘+1=𝐮^𝑘_𝑎+𝜇₂𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐚^𝑘+1,𝐮_𝑏^𝑘+1=𝐮^𝑘_𝑏+𝜇₃𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐛^𝑘+1,𝐯^𝑘+1=𝐯^𝑘+𝜇₄𝛽^𝑘+1−𝑡^𝑘+1,𝑞^𝑘+1=𝑞^𝑘+𝜇₅1𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽^𝑘+1+𝛽₀^𝑘+1. (12)

The efficiency of the iterative algorithm (12) lies on whether the first equation of (12) can be solved quickly. According to the theory of ADMM, these variables (𝛽,𝛽₀), 𝐡,𝐚,𝐛, and 𝐭 are updated in an alternating or sequential fashion, which accounts for the term alternating direction. So we can get𝛽^𝑘+1,𝛽₀^𝑘+1=argmin_𝛽,𝛽0𝐿𝛽,𝛽₀,𝐡^𝑘,𝐚^𝑘,𝐛^𝑘,𝐭^𝑘,𝐮^𝑘_ℎ,𝐮^𝑘_𝑎,𝐮^𝑘_𝑏,𝐯^𝑘,𝑞^𝑘,𝐡^𝑘+1=argmin_𝐡𝐿𝛽^𝑘+1,𝛽₀^𝑘+1,𝐡,𝐚^𝑘,𝐛^𝑘,𝐭^𝑘,𝐮^𝑘_ℎ,𝐮^𝑘_𝑎,𝐮^𝑘_𝑏,𝐯^𝑘,𝑞^𝑘,𝐚^𝑘+1=argmin_𝐚𝐿𝛽^𝑘+1,𝛽₀^𝑘+1,𝐡^𝑘+1,𝐚,𝐛^𝑘,𝐭^𝑘,𝐮^𝑘_ℎ,𝐮^𝑘_𝑎,𝐮^𝑘_𝑏,𝐯^𝑘,𝑞^𝑘,𝐛^𝑘+1=argmin_𝐛𝐿𝛽^𝑘+1,𝛽₀^𝑘+1,𝐡^𝑘+1,𝐚^𝑘+1,𝐛,𝐭^𝑘,𝐮^𝑘_ℎ,𝐮^𝑘_𝑎,𝐮^𝑘_𝑏,𝐯^𝑘,𝑞^𝑘,𝐭^𝑘+1=argmin_𝐭𝐿𝛽^𝑘+1,𝛽₀^𝑘+1,𝐡^𝑘+1,𝐚^𝑘+1,𝐛^𝑘+1,𝐭,𝐮^𝑘_ℎ,𝐮^𝑘_𝑎,𝐮^𝑘_𝑏,𝐯^𝑘,𝑞^𝑘. (13)

For the first equation in (13), it is equivalent to the following convex optimization: 𝛽^𝑘+1,𝛽₀^𝑘+1=argmin_𝛽,𝛽0𝜆₂2𝛽²₂+𝐮^𝑇_ℎ𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙−𝐡+𝐮^𝑇_𝑎𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢−𝐚+𝐮^𝑇_𝑏𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢−𝐛+𝐯^𝑇1𝛽−𝐭+𝑞𝑟−𝑢^𝑢_𝑖=1𝑥_𝑖𝛽+𝛽₀+𝜇₁2𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙−𝐡²₂+𝜇₂2𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢−𝐚²₂+𝜇₃2𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽+𝛽₀𝟏_𝑢−𝐛²₂+𝜇₄2𝛽−𝐭²₂+𝜇₅21𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽+𝛽₀²₂. (14)

The objective function in the above minimization problem is quadratic and differentiable, and since (𝛽^𝑘+1,𝛽₀^𝑘+1) minimizes this function by definition, the optimal solution can be found by solving a set of linear equations𝜆₂+𝜇₄+𝜇₅𝑢^2𝑛_𝑖=∑_𝑙+1𝑥∑^𝑇_𝑖^𝑛_𝑖=𝑙+1𝑥_𝑖𝐈+𝜇₁𝑋^𝑇_𝑙𝑋_𝑙+𝜇₂+𝜇₃𝑋^𝑇_𝑢𝑋_𝑢𝜇₁𝑋^𝑇_𝑙𝟏_𝑙+𝜇₂+𝜇₃𝑋^𝑇_𝑢𝟏_𝑢+𝜇₅𝑢^𝑛_𝑖=∑_𝑙+1𝑥^𝑇_𝑖𝟏^𝑇_𝑙𝑋_𝑙+𝜇₂+𝜇₃𝟏^𝑇_𝑢𝑋_𝑢+𝜇₅𝑢^𝑛_𝑖=∑_𝑙+1𝑥_𝑖𝟏_𝑢𝑙𝜇₁+2𝑢𝜇₂+𝜇₃+𝜇₅𝛽^𝑘+1𝛽₀^𝑘+1=𝑋^𝑇_𝑙𝑌^𝑇_𝑙𝐮^𝑘_ℎ+𝜇₁1−𝐡^𝑘+𝑋^𝑇_𝑢𝑌^𝑇_𝑢𝐮^𝑘_𝑎+𝐮^𝑘_𝑏+𝜇₂1−𝐚^𝑘+𝜇₃𝑠−𝐛^𝑘−𝐯^𝑘+𝑞^𝑘+𝜇₅𝑟∑^𝑛_𝑖=𝑙+1𝑥^𝑇_𝑖𝑢+𝜇₄𝐭^𝑘𝟏^𝑇_𝑙𝑌^𝑇_𝑙𝐮^𝑘_ℎ+𝜇₁1−𝐡^𝑘+𝟏^𝑇_𝑢𝑌^𝑇_𝑢𝐮^𝑘_𝑎+𝐮^𝑘_𝑏+𝜇₂1−𝐚^𝑘+𝜇₃𝑠−𝐛^𝑘+𝑞^𝑘+𝜇₅𝑟. (15)

In (15), 𝐈 is a 𝑝×𝑝 unit matrix and the coefficient matrix is a (𝑝+1)×(𝑝+1) matrix, independent of the optimization variables. For large 𝑝, small 𝑛 setting, the term 𝑋^𝑇_𝑙𝑋_𝑙 in the coefficient matrix will be a positive low rank matrix with rank at most 𝑙 while the term 𝑋^𝑇_𝑢𝑋_𝑢 in the coefficient matrix will be a positive low rank matrix with rank at most 2𝑢. Therefore, the coefficient matrix is also low rank matrix with rank at most (2𝑢+1). And if we use CG to solve the problem (15), it will converge in less than (2𝑢+1) steps [21].

For the second equation in (13), it is equivalent to solving 𝐡^𝑘+1=argmin_𝐡𝐶^𝑙_𝑖=1ℎ_𝑖++𝐮^𝑘_ℎ^𝑇𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘+1+𝜇₁2𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘+12₂. (16)

In order to solve (16), we need the following Proposition [22].

Proposition 2. Let 𝑠_𝜆(𝑤)=argmin_𝑥∈ℝ𝜆𝑥₊+0.5‖𝑥−𝑤‖²₂ where 𝜆>0. Then 𝑠_𝜆(𝑤)=𝑤−𝜆𝑤>𝜆,00≤𝑤≤𝜆,𝑤𝑤<0. (17)

Combined with Proposition 2 and 𝐮^𝑘_ℎ²₂2𝜇₁+𝐮^𝑘_ℎ^𝑇𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘+1+0.5𝜇₁𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘+12₂=0.5𝜇₁𝐮^𝑘_ℎ𝜇₁+𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘+12₂, (18)

we can update 𝐡^𝑘+1 according to Corollary 3.

Corollary 3. The update of 𝐡^𝑘+1 in (16) is equivalent to 𝐡^𝑘+1=𝑆_𝐶/𝜇1𝟏_𝑙+𝐮^𝑘_ℎ𝜇₁−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙, (19) where 𝑆_𝜆(𝑤)=(𝑠_𝜆(𝑤₁),…,𝑠_𝜆(𝑤_𝑙)).

For the third equation in (13), it is equivalent to solving (20): 𝐚^𝑘+1=argmin_𝐚𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝑎_𝑖++𝐮^𝑘_𝑎^𝑇𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐚^𝑘+1+𝜇₂2𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐚^𝑘+12₂. (20)

Combined with Proposition 2 and 𝐮^𝑘_𝑎²₂2𝜇₁+𝐮^𝑘_𝑎^𝑇𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐚^𝑘+1+0.5𝜇₂𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐚^𝑘+12₂=0.5𝜇₂𝐮^𝑘_𝑎𝜇₂+𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐚^𝑘+12₂, (21)

we can update 𝐚^𝑘+1 according to Corollary 4.

Corollary 4. The update of 𝐚^𝑘+1 in (20) is equivalent to 𝐚^𝑘+1=𝑆_𝐶^∗_/𝜇2𝟏_𝑢+𝐮^𝑘_𝑎𝜇₂−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢, (22) where 𝑆_𝜆(𝑤)=(𝑠_𝜆(𝑤₁),…,𝑠_𝜆(𝑤_2𝑢)).

For the fourth equation in (13), it is equivalent to solving 𝐛^𝑘+1=argmin_𝐛𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝑏_𝑖++𝐮^𝑘_𝑏^𝑇𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐛^𝑘+1+𝜇₃2𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐛^𝑘+12₂. (23) In order to solve (23), we need the following proposition.

Proposition 5. Let 𝑔_𝜆(𝑤)=argmin_𝑥∈ℝ𝜆𝑥₊+0.5‖𝑥−𝑤‖²₂, where 𝜆<0. Then 𝑔_𝜆(𝑤)=𝑤−𝜆𝑤>0,𝑤𝑤<0. (24)

Proof. The function 𝜆𝑥₊+0.5‖𝑥−𝑤‖²₂ is strongly convex, so it has a unique solution. Therefore, by the subdifferential calculus [23], 𝑔_𝜆 is the unique minimizer of the following equation: 0∈𝜆𝜕𝑥₊+𝑥−𝑤, (25) where 𝜕(𝑥₊) is the subdifferential of the function 𝑥₊. According to [23], 𝜕(𝑥₊) can be expressed as 𝜕𝑥₊=[1𝑥>0,𝑝∶𝑝∈0,1]𝑥=0,0𝑥<0. (26)
With (25) and (26), we can get the desired result.
Then, combined with Proposition 5 and 𝐮^𝑘_𝑏²₂2𝜇₃+𝐮^𝑘_𝑏^𝑇𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐛^𝑘+1+0.5𝜇₃𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐛^𝑘+12₂=0.5𝜇₃𝐮^𝑘_𝑏𝜇₃+𝑠𝟏_𝑢−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢−𝐛^𝑘+12₂, (27) we can update 𝐛^𝑘+1 according to Corollary 6.

Corollary 6. The update of 𝐛^𝑘+1 in (23) is equivalent to 𝐛^𝑘+1=𝐺_−𝐶^∗_/𝜇3𝑠𝟏_𝑢+𝐮^𝑘_𝑏𝜇₃−𝑌_𝑢𝑋_𝑢𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑢, (28) where 𝐺_𝜆(𝑤)=(𝑔_𝜆(𝑤₁),…,𝑔_𝜆(𝑤_2𝑢)).

For the fifth equation in (13), it is equivalent to solving 𝐭^𝑘+1=argmin_𝐭𝜆₁𝐭^𝑘+1₁+𝐯^𝑘𝑇𝛽^𝑘+1−𝐭^𝑘+1+0.5𝜇₄𝛽^𝑘+1−𝐭^𝑘+12₂. (29)

Solving (29) can be done efficiently using soft threshold, and we can update 𝐭^𝑘+1 according to Corollary 7.

Corollary 7. The update of 𝐭^𝑘+1 in (29) is 𝐭^𝑘+1=𝑇_𝜆1/𝜇4𝐯^𝑘𝜇₄+𝛽^𝑘+1, (30) where 𝑇_𝜆(𝑤)=(𝑡_𝜆(𝑤₁),…,𝑡_𝜆(𝑤_𝑝)), and 𝑡_𝜆(𝑤)=sgn(𝑤)max{0,|𝑤|−𝜆}.

According to the theory of ADMM, ‖𝑠^𝑘₁‖₂, ‖𝑠^𝑘₂‖₂, ‖𝑠^𝑘₃‖₂, and ‖𝑠^𝑘₄‖₂ must be small below some certain threshold 𝛿. Considering the auxiliary variables in (11), we expect that ‖𝟏_𝑙−𝑌_𝑙(𝑋_𝑙𝛽+𝛽₀𝟏_𝑙)−𝐡‖₂,‖𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢)−𝐚‖₂, ‖𝑠𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢)−𝐛‖₂, and ‖𝛽−𝐭‖₂ also must be small. Therefore, in our experiment, algorithm for SENFS stops whenever 𝑠^𝑘₁₂<𝛿,𝑠^𝑘₂₂<𝛿,𝑠^𝑘₃₂<𝛿,𝑠^𝑘₄₂𝟏<𝛿,_𝑙−𝑌_𝑙(𝑋_𝑙𝛽+𝛽₀𝟏_𝑙)−𝐡₂<√𝟏𝑙𝛿,_𝑢−𝑌_𝑢(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢)−𝐚₂<√2𝑢𝛿,𝑠𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢)−𝐛₂<√2𝑢𝛿,𝛽−𝐭₂<√𝑝𝛿. (31)

Finally, we can get the algorithm for SENFS. The detailed procedure of the algorithm ADMM for SENFS is summarized in Algorithm 8 as follows.

Algorithm 8. ADMM algorithm for SENFS. Input. Labeled data set {𝑥_𝑖,𝑦_𝑖}^𝑙_𝑖=1; unlabeled data set {𝑥_𝑖}^𝑛_𝑖=𝑙+1; tuning parameters 𝜆₁ and 𝜆₂; regularization parameter 𝐶. Output. Selected feature set. Step  1. Initialize 𝛽⁰,𝛽⁰₀,𝐡⁰,𝐚⁰,𝐛⁰,𝐭⁰,𝐮⁰_ℎ,𝐮⁰_𝑎,𝐮⁰_𝑏,𝐯⁰, and 𝑞⁰. Step  2. If (31) is satisfied, go to Step  3; otherwise,(1)update 𝛽^𝑘+1 and 𝛽₀^𝑘+1 according to (15);(2)update 𝐡^𝑘+1, 𝐚^𝑘+1, 𝐛^𝑘+1, and 𝐭^𝑘+1 according to Corollaries 3–7, respectively;(3)update 𝐮_ℎ^𝑘+1,𝐮_𝑎^𝑘+1,𝐮_𝑏^𝑘+1,𝐯^𝑘+1, and 𝑞^𝑘+1 according to (12). Step  3. Get the best feature subset according to 𝛽_𝑗^𝑘+1,𝑗=1,…,𝑝. If 𝛽_𝑗^𝑘+1=0, the corresponding 𝑗th feature is abandoned; otherwise, it is selected as an important feature.

3.2. Convergence Analysis and Computational Cost

The convergence property of Algorithm 8 can be derived from the theory of the alternating direction method of multipliers. According to the standard convergence theory of ADMM, Algorithm 8 satisfies the dual variable convergence [24]. So Theorem 9 holds.

Theorem 9. Suppose that (𝛽^∗,𝛽₀^∗) is one of solution of (5). Then the following property holds:𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^𝑘+𝛽^𝑘₀+𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^𝑘+𝛽^𝑘₀−𝐻_𝑠𝑦_𝑖𝑥_𝑖𝛽^𝑘+𝛽^𝑘₀+𝜆₁𝛽^𝑘₁+𝜆₂𝛽^𝑘2₂2=𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀+𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀−𝐻_𝑠𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀+𝜆₁𝛽^∗₁+𝜆₂𝛽^∗2₂2. (32)

As for the computational issue, it is hard to predict the computational cost because it depends on the all the penalty parameters. According to our experience, we only need to iterate a few hundred iterations to get a reasonable result. On the other hand, the efficiency of Algorithm 8 lies mainly on whether we can quickly solve the linear equations (15). And the computational cost for solving (15) is 𝑂(𝑙²𝑝+4𝑢²𝑝).

3.3. Varying Penalty Parameter

In order to make performance less dependent on the initial choice of the penalty parameter, it is necessary to use different penalty parameters. According to our experiment experience, the penalty parameters 𝜇₁,𝜇₂,𝜇₃, and 𝜇₄ have a huge influence on the performance and the number of iterations involved, so adaptive selections of them are performed.

For 𝜇₁, let 𝑓(𝛽,𝛽₀)=𝐶∑^𝑙_𝑖=1(ℎ_𝑖)₊+𝐶^∗∑^𝑙+2𝑢_𝑖=𝑙+1(𝐻₁(𝑦_𝑖(𝑥_𝑖𝛽+𝛽₀))−𝐻_𝑠(𝑦_𝑖(𝑥_𝑖𝛽+𝛽₀)))+𝜆₁‖𝛽‖₁+0.5𝜆₂‖𝛽‖²₂, and the constraint conditions are 𝐡=𝟏_𝑙−𝑌_𝑙(𝑋_𝑙𝛽+𝛽₀𝟏_𝑙) and the constraint of (3). The optimization task (6) is equivalent to 𝑊𝛽,𝛽₀=𝑓𝛽,𝛽₀+𝐮^𝑇_ℎ𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙1−𝐡+𝑞𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽+𝛽₀+𝜇₁2𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽+𝛽₀𝟏_𝑙−𝐡²₂+𝜇₅21𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽+𝛽₀²₂. (33)

The necessary optimality conditions for the problem (6) are dual feasibility as 𝛽0∈𝜕𝑓^∗,𝛽^∗₀𝜕𝛽^∗−𝑋^𝑇_𝑙𝑌^𝑇_𝑙𝐮^∗_ℎ−𝑞^∗𝑢^𝑛_𝑖=𝑙+1𝑥^𝑇_𝑖. (34) Since 𝛽^𝑘+1 minimizes 𝑊(𝛽^𝑘+1,𝛽₀^𝑘+1,𝐮^𝑘_ℎ,𝑞^𝑘) by definition, we have that 𝛽0∈𝜕𝑓^𝑘+1,𝛽₀^𝑘+1𝜕𝛽^𝑘+1−𝑋^𝑇_𝑙𝑌^𝑇_𝑙𝐮^𝑘_ℎ−𝜇₁𝟏_𝑙−𝑌_𝑙𝑋_𝑙𝛽^𝑘+1+𝛽₀^𝑘+1𝟏_𝑙−𝐡^𝑘−𝑞^𝑘𝑢^𝑛_𝑖=𝑙+1𝑥^𝑇_𝑖−𝜇₅𝑢×^𝑛_𝑖=𝑙+1𝑥^𝑇_𝑖1𝑟−𝑢^𝑛_𝑖=𝑙+1𝑥_𝑖𝛽^𝑘+1+𝛽₀^𝑘+1=𝜕𝑓𝛽^𝑘+1,𝛽₀^𝑘+1𝜕𝛽^𝑘+1−𝑋^𝑇_𝑙𝑌^𝑇_𝑙𝐮_ℎ^𝑘+1−𝑞^𝑘+1𝑢×^𝑛_𝑖=𝑙+1𝑥^𝑇_𝑖+𝜇₁𝑋^𝑇_𝑙𝑌^𝑇_𝑙𝐡^𝑘−𝐡^𝑘+1. (35)

Compared with (34), (35) means that the quantity 𝜇₁𝑋^𝑇_𝑙𝑌^𝑇_𝑙(𝐡^𝑘−𝐡^𝑘+1) can be viewed as a residual for (34), and let 𝑠₁^𝑘+1=𝜇₁𝑋^𝑇_𝑙𝑌^𝑇_𝑙(𝐡^𝑘−𝐡^𝑘+1).

For 𝜇₂, let 𝑓(𝛽,𝛽₀) = 𝐶∑^𝑙_𝑖=1𝐻₁(𝑦_𝑖(𝑥_𝑖𝛽+𝛽₀)) + 𝐶^∗∑^𝑙+2𝑢_𝑖=𝑙+1((𝑎_𝑖)₊−𝐻_𝑠(𝑦_𝑖(𝑥_𝑖𝛽+𝛽₀)))+𝜆₁‖𝛽‖₁+0.5𝜆₂‖𝛽‖²₂, and the constraint conditions are 𝐚=𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽+𝛽₀𝟏_𝑢) and the constraint of (3). Through the same solving process as parameter 𝜇₁, we can get the residual 𝑠₂^𝑘+1=𝜇₂^𝑘+1𝑋^𝑇_𝑢𝑌^𝑇_𝑢(𝐚^𝑘−𝐚^𝑘+1). Similarly, we can get the residual 𝑠₃^𝑘+1=𝜇₃^𝑘+1𝑋^𝑇_𝑢𝑌^𝑇_𝑢(𝐛^𝑘−𝐛^𝑘+1) for parameter 𝜇₃ and 𝑠₄^𝑘+1=𝜇₄^𝑘+1𝑋^𝑇_𝑢𝑌^𝑇_𝑢(𝐭^𝑘+1−𝐭^𝑘) for parameter 𝜇₄. With these residuals, we can get a simple scheme to update 𝜇₁,𝜇₂,𝜇₃, and 𝜇₄, respectively, according to Corollary 10.

Corollary 10. The update of 𝜇₁,𝜇₂,𝜇₃,𝜇₄ is𝜇₁^𝑘+1=𝜏^{𝑖𝑛𝑐𝑟}𝜇^𝑘₁𝑖𝑓𝟏_𝑙−𝑌_𝑙(𝑋_𝑙𝛽^𝑘+𝛽^𝑘₀𝟏_𝑙)−𝐡^𝑘₂>𝜃𝑠^𝑘₁₂𝜇^𝑘₁𝜏^{𝑑𝑒𝑐𝑟}𝑖𝑓𝑠^𝑘₁₂>𝜃𝟏_𝑙−𝑌_𝑙(𝑋_𝑙𝛽^𝑘+𝛽^𝑘₀𝟏_𝑙)−𝐡^𝑘₂𝜇^𝑘₁𝜇𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒,_2𝑘+1=𝜏^{𝑖𝑛𝑐𝑟}𝜇^𝑘₂𝑖𝑓𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽^𝑘+𝛽^𝑘₀𝟏_𝑢)−𝐚^𝑘₂>𝜃𝑠^𝑘₂₂𝜇^𝑘₂𝜏^{𝑑𝑒𝑐𝑟}𝑖𝑓𝑠^𝑘₂₂>𝜃𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽^𝑘+𝛽^𝑘₀𝟏_𝑢)−𝐚^𝑘₂𝜇^𝑘₂𝜇𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒,_3𝑘+1=𝜏^{𝑖𝑛𝑐𝑟}𝜇^𝑘₃𝑖𝑓𝑠𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽^𝑘+𝛽^𝑘₀𝟏_𝑢)−𝐛^𝑘₂>𝜃𝑠^𝑘₃₂𝜇^𝑘₃𝜏^{𝑑𝑒𝑐𝑟}𝑖𝑓𝑠^𝑘₃₂>𝜃𝑠𝟏_𝑢−𝑌_𝑢(𝑋_𝑢𝛽^𝑘+𝛽^𝑘₀𝟏_𝑢)−𝐛^𝑘₂𝜇^𝑘₃𝜇𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒,_4𝑘+1=𝜏^{𝑖𝑛𝑐𝑟}𝜇^𝑘₄𝑖𝑓𝛽−𝐭₂>𝜃𝑠^𝑘₄₂𝜇^𝑘₄𝜏^{𝑑𝑒𝑐𝑟}𝑖𝑓𝑠^𝑘₄₂>𝜃𝛽−𝐭₂𝜇^𝑘₄𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒, (36)where 𝜏^{𝑖𝑛𝑐𝑟}>1,𝜏^{𝑑𝑒𝑐𝑟}>,  and 𝜃>1 are parameters. Typical choices might be 𝜏^{𝑖𝑛𝑐𝑟}=𝜏^{𝑑𝑒𝑐𝑟}=2 and 𝜃=10.

4. Experimental Evaluation

This section examines the performance of SENFS with respect to its feature selection and test error on simulated data and six benchmark data sets. In order to evaluate the effectiveness of SENFS, we compare SENFS with an existing semisupervised feature selection algorithm: Spectral [10], and a supervised feature selection algorithm: DrSVM [11], which also has the characteristics of grouping effect. On the other hand, in order to evaluate the quality of selected features, SVM was executed on these selected features. The experiments are run on a desktop with Pentium(R) 2.0 G CPU, 1.99 G main memory. The programs are compiled in Windows system with Matlab in version R2009a.

The limited number of samples prohibits having enough and independent training and testing data for performance evaluation. It is very common to apply accross-validation (CV) in this scenario. We used 5-fold CV: we partitioned the data set into five complementary subsets of equal size. Four subsets were used as training data; the remaining subset served as test data. We repeated this process five times such that each of the five subsets was used exactly once as test data. To get more reliable estimate, we performed the 5-fold CV for 10 times and the experimental results are average results over test data sets. Moreover, finding the appropriate value of the tuning parameter pair 𝜆₁ and 𝜆₂ is essential for the performance of SENFS. We employed 10-fold CV over a large grid.

4.1. Simulation

We evaluate the performance of SENFS using two parameters: the correlation between relevant features denoted by 𝜌, the number of labeled samples, and the degree of overlapping among classes denoted by 𝑑. Consider 2-class problem in which the samples are lying in a 𝑝 dimensional space with the first 10 dimensional being relevant to classification and the remaining features being noise, where the correlation between the first 10 features is 𝜌. The number of samples is 300 with 𝑝=500. For the samples from +1 class, they are sampled from a normal distribution with mean and covariance as follows: 𝛼₊=𝑑,…,𝑑₁₀,0,…,0_𝑝−10,=^∗₁₀∑_×100_{10×(𝑝−10)}0_{(𝑝−10)×10}𝐈_{(𝑝−10)×(𝑝−10)}, (37) where the diagonal elements of ∑^∗ are 1 and the off-diagonal elements are all equal to 𝜌. The −1 class has a similar distribution expect that its mean is: 𝛼₋=−𝑑,…,−𝑑₁₀,0,…,0_𝑝−10. (38)

To evaluate the effect of the correlation between relevant features, SENFS is compared with Spectral and DrSVM, measured by the number of selected features with two labeled samples and 𝑑=1. The results are summarized in Table 1. As shown in Table 1, on this simulated data, when the relevant features are highly correlated (e.g., 𝜌=0.9), Spectral and DrSVM tend to keep only a small subset of the relevant correlated variables and overlook the others, while the SENFS tends to identify all of them, due to the grouping effect. These three methods seem to work well in removing irrelevant features.

Table 1 

Comparisons of selected features and their standard errors (in parenthesis). 𝐹_relevant is the number of selected relevant features and 𝐹_noise is the number of selected noise features.

The effects of the number of labeled samples on test error over the top 10 selected features are summarized in Figure 1 with 𝑑=1 and 𝜌=0.9. As can be seen, the test errors of SENFS, Spectral, and DrSVM decrease with the increase of the number of labeled samples, but SENFS seems to achieve the best classification performance when the number of labeled samples is varying, which may imply that SENFS can make better use of the labeled samples than spectral and DrSVM. The supervised feature selection method DrSVM achieves the worst results because it only relies on the few labeled samples and discards the large amount of unlabeled samples.

Figure 1 

Test error versus the number of labeled samples.

In Table 2, the effect of the degree of overlapping among classes on test error over the top 5 selected features is evaluated with two labeled samples and 𝑑=1 and 𝜌=0.9, also reporting the typical computational time of our experimental campaign. As we can see, SENFS seems to have the best prediction performance. When 𝑑 is small, the two classes overlap largely and in this case, other methods achieved worse performance compared with SENFS. However, SENFS, solving the programming problem (6) needing an iterative procedure, requires more computational time than the other methods as you can see in Table 2. It is noted that the absolute values are not as important as the relative differences between the individual methods.

Table 2 

Comparisons of test errors, computational time needed, and their standard errors (in parenthesis).

4.2. Application to Benchmark Data Sets

Several benchmark data sets are selected to test the performance of SENFS, which are used as benchmark data sets in [7, 8] to test the performances of semisupervised algorithms. These benchmark data sets consist of 9 semisupervised learning data sets. We did not test the SSL6, SSL8, and SSL9 data sets since the SSL6 data set includes six classes, the SSL8 data set contains too many samples (𝑛 is over one million) and the SSL9 data set has too many dimensions (𝑝 is over ten thousand). The names and characteristics of the left six data sets are given in Table 3.

Table 3 

Data sets used in the experiments. 𝑛, 𝑝, and 𝑙 are the number of samples, features, and labeled samples, respectively.

In this study, we examine performance evaluation through 5-fold cross-validation that is, we randomly select four fifths of the unlabeled samples, plus all the labeled samples, for SENFS, Spectral, and DrSVM to select optimal feature subsets, while leaving the remaining one fifth for testing test error on the selected features using SVM, where all the labeled samples are used for training SVM. The results measured by test error are reported in Table 4. As can be seen, SENFS outperforms the semisupervised and supervised feature selection methods on all the six data sets when 𝑙=10 and 𝑙=100. When 𝑙=10, Spectral performs the second best, on USPS, COIL2, and BCI data sets, while DrSVM performs the second best on Digit1, BCI, g241c and g241n data sets when 𝑙=100.

Table 4 

Comparisons of test errors and their standard errors (in parenthesis) on benchmark data sets.

5. Conclusion

This paper has proposed a novel semisupervised feature selection algorithm based on SVM and the elastic net penalty. The whole methodology of SENFS and the solution path based on ADMM have been described in detail in this paper. The experimental results illustrate that SENFS can identify the relevant features and encourage highly correlated features to be selected or removed together.

Future work will address how these selected features interpret their semantic relationship with the data they are selected from, which can be used for unknown data analysis, and extend SENFS to be suitable for multiclass case.

Appendix

Proof of Theorem 1. Consider another set of coefficients 𝛽^∗₀=𝛽₀,𝛽^∗_𝑗′=𝛽_𝑗+𝛽_𝑘2,if𝑗^′=𝑗or𝑗^′𝛽=𝑘,_𝑗′otherwise. (A.1) Then we have 𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀+𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀−𝐻_𝑠𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀+𝜆₁𝛽^∗₁+𝜆₂𝛽^∗2₂2≥𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽+𝛽₀+𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑥_𝑖𝛽+𝛽₀−𝐻_𝑠𝑦_𝑖𝑥_𝑖𝛽+𝛽₀+𝜆₁𝛽₁+𝜆₂𝛽²₂2. (A.2)
It is simple to verify that both the loss function 𝐻₁(𝑧) and 𝐻_𝑠(𝑧) are Lipschitz continuous, so we have 𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀−𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽+𝛽₀≤𝐶^𝑙_𝑖=1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀−𝐻₁𝑦_𝑖𝑥_𝑖𝛽+𝛽₀≤𝑀^𝑙_𝑖=1𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀−𝑦_𝑖𝑥_𝑖𝛽+𝛽₀=𝑀2𝛽_𝑗−𝛽_𝑙^𝑙_𝑖=1𝑥_𝑖𝑗−𝑥_𝑖𝑘. (A.3)
Similarly, we can get 𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀−𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻₁𝑦_𝑖𝑥_𝑖𝛽+𝛽₀≤𝑀₁2𝛽_𝑗−𝛽_𝑘^𝑙+2𝑢_𝑖=𝑙+1𝑥_𝑖𝑗−𝑥_𝑖𝑘,𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻_𝑠𝑦_𝑖𝑥_𝑖𝛽+𝛽₀−𝐶^{∗𝑙+2𝑢}_𝑖=𝑙+1𝐻_𝑠𝑦_𝑖𝑥_𝑖𝛽^∗+𝛽^∗₀≤𝑀₂2𝛽_𝑗−𝛽_𝑘^𝑙+2𝑢_𝑖=𝑙+1𝑥_𝑖𝑗−𝑥_𝑖𝑘, (A.4) where 𝑀₁ and 𝑀₂ are positive constants. As described in [11], we get ‖𝛽^∗‖₁−‖𝛽‖₁≤0 and ‖𝛽^∗‖²₂−‖𝛽‖²₂=−(1/2)|𝛽_𝑗−𝛽_𝑘|. Then combining (A.3), (A.4), and (A.2) implies that 𝑀2𝛽_𝑗−𝛽_𝑘^𝑙_𝑖=1𝑥_𝑖𝑗−𝑥_𝑖𝑘+𝑀₁+𝑀₂2𝛽_𝑗−𝛽_𝑘×^𝑙+2𝑢_𝑖=𝑙+1𝑥_𝑖𝑗−𝑥_𝑖𝑘−𝜆₂4𝛽_𝑗−𝛽_𝑘²≥0. (A.5)
Let 𝑀^′=𝑀₁+𝑀₂, and (7) is obtained. For (8), we simply use the two inequalities ∑^𝑙_𝑖=1|𝑥_𝑖𝑗−𝑥_𝑖𝑘|≤√𝑙2(1−𝜌_𝑙), and ∑^𝑙+2𝑢_𝑖=𝑙+1|𝑥_𝑖𝑗−𝑥_𝑖𝑘|≤√2𝑢2(1−𝜌_𝑢).

Acknowledgments

This work was supported by the National Key Basic Research Program of China (973 Program) under Grant no. 613148 and the National Science and Technology Major Project of the Ministry of Science and Technology under Grant no. 2010ZX03006-002, 2011ZX03005-003-03. The authors are grateful to the anonymous reviewers for their helpful comments.