Robust multi-class feature selection via l 2,0 -norm regularization minimization

3.1 Problem formulation

In this work, we use the least squares regression combined with a l 2 , 0 -norm regularization for multi-class feature selection, the optimization function is formulated as

φ λ ⁢ ( 𝐖 ) = min 𝐖 ⁢ 1 2 ⁢ || 𝐖 T ⁢ 𝐗 + 𝐛 𝟏 T - 𝐘 || F 2 + λ ⁢ || 𝐖 || 2 , 0 ,(4)

where the notations are defined as in 2.2.

According to the Karush-Kuhn-Tucker (KKT) theorem [35], set the derivative of Eq. (4) with respect to b to be 0, then the optimal solution of b can be obtained as follows:

𝐛 = 1 N ⁢ ( 𝐘 - 𝐖 T ⁢ 𝐗 ) ⁢ 𝟏 .(5)

Replacing b with Eq. (5), we can rewrite Eq. (4) as

φ λ ⁢ ( 𝐖 ) = min 𝐖 ⁢ 1 2 ⁢ || ( 𝐖 T ⁢ 𝐗 - 𝐘 ) ⁢ 𝐇 || F 2 + λ ⁢ || 𝐖 || 2, ⁢ 0 ,(6)

where 𝐇 = 𝐈 - ( 1 / N ) ⁢ 𝟏𝟏 T is referred as the centering matrix, which has following properties

𝐇 = 𝐇 T = 𝐇𝐇 T = 𝐇 T ⁢ 𝐇 .(7)

Now, the objective function of l 2 , 0 -FS is presented totally, next we will show how to optimize it. Iterative hard thresholding (IHT) algorithm [36, 37, 38] is usually used to solve l 0 -norm regularization minimization for the vector sparsity problem. Inspired by it, we extended IHT to solve the matrix sparsity Eq. (6).

Denoting 𝐗 ∼ = 𝐗𝐇 and 𝐘 ∼ = 𝐘𝐇 , Eq. (6) can be reformulated as

φ λ ⁢ ( 𝐖 ) = min 𝐖 ⁢ 1 2 ⁢ || 𝐖 T ⁢ 𝐗 ∼ - 𝐘 ∼ || F 2 + λ ⁢ || 𝐖 || 2 , 0 .(8)

IHT algorithms use the proximal gradient technique to update the solution iteratively. Firstly, we rewrite Eq. (8) as

φ λ ⁢ ( 𝐖 ) = min 𝐖 ⁢ f ⁢ ( 𝐖 ) + λ ⁢ || 𝐖 || 2 , 0 ,(9)

where f ⁢ ( 𝐖 ) = 1 2 ⁢ || 𝐖 T ⁢ 𝐗 ∼ - 𝐘 ∼ || F 2 is a differentiable convex function, and its gradient ∇ ⁡ f ⁢ ( 𝐖 ) satisfies the Lipschitz continuous condition (denote its Lipschitz constant as L f ). Therefore, f ⁢ ( 𝐖 ) can be approximately iteratively updated by the projected gradient method

𝐖 t + 1 = arg ⁢ min 𝐖 f ⁢ ( 𝐖 t ) + t ⁢ r ⁢ ( ∇ ⁡ f ⁢ ( 𝐖 t ) T ⁢ ( 𝐖 - 𝐖 t ) ) + L 2 ⁢ || 𝐖 - 𝐖 t || F 2 ,(10)

where L ⩾ 0 is a constant, which should satisfy the condition of L ⩾ L f .

Adding λ ⁢ || 𝐖 || 2 , 0 into both sides of Eq. (10), Eq. (8) can be solved by iteratively optimizing the subproblem

𝐖 t + 1 = arg ⁢ min 𝐖 f ⁢ ( 𝐖 t ) + t ⁢ r ⁢ ( ∇ ⁡ f ⁢ ( 𝐖 t ) T ⁢ ( 𝐖 - 𝐖 t ) ) + L 2 ⁢ || 𝐖 - 𝐖 t || F 2 + λ ⁢ || 𝐖 || 2 , 0 .(11)

[b] Homotopy iterative hard thresholding method to solve Eq. (4)(Input:) Training data 𝐗 ∈ R d × N , training labels 𝐘 ∈ R C × N , centering matrix 𝐇 ∈ R N × N , parameters L 0 , λ 0 , λ t ⁢ g ⁢ t ; (Output:) 𝐖 * ; 1: initialize ${\text{𝐖}}^0=\text{𝟎},\rho \in (0,1),ϵ>0,\eta >0,\gamma >1,k\leftarrow 0$ ; 2: do3: t ← 0 ;4: 𝐖 k , 0 = 𝐖 k ;5: L k , 0 ← L k ;6: do An L-tuning iteration indexed by t 7: update 𝐖 k , t + 1 by Eq. (14);8: while do9: L k , t ← γ ⁢ L k , t ;10: update 𝐖 k , t + 1 by Eq. (14);11: end while12: L k , t + 1 ← L k , t ;13: t ← t + 1 ;14: until || 𝐖 k , t - 𝐖 k , t + 1 || F 2 ⩽ ϵ 15: 𝐖 k + 1 ← 𝐖 k , t ;16: L k + 1 ← L k , t .17: λ k + 1 ← ρ ⁢ λ k ;18: k ← k + 1 ;19: until λ k + 1 ⩽ λ 𝑡𝑔𝑡 ;20: 𝐖 ⋆ ← 𝐖 k .

By removing the item f ⁢ ( 𝐖 t ) and adding an item 1 2 ⁢ L ⁢ || ∇ ⁡ f ⁢ ( 𝐖 t ) || F 2 , both of which are independent on W) and thus can be seen as constant items, the right hand of Eq. (11) can be rewritten as

𝐖 t + 1 = arg ⁢ min 𝐖 L 2 ⁢ { || 𝐖 - ( 𝐖 t - 1 L ⁢ ∇ ⁡ f ⁢ ( 𝐖 t ) ) || F 2 + 2 ⁢ λ L ⁢ || 𝐖 || 2 , 0 } ,(12)

Since the Frobenius norm and l 2 , 0 -norm are all separable function (both of these two functions can be represented as a sum of single-variable functions), each row of W can be updated individually

( 𝐰 i ) t + 1 = arg ⁢ min 𝐰 i L 2 ⁢ { || 𝐰 i - ( ( 𝐰 i ) t - 1 L ⁢ ∇ ⁡ f ⁢ ( ( 𝐰 i ) t ) ) || 2 2 + 2 ⁢ λ L ⁢ 1 || 𝐰 i || 2 ≠ 0 } ,(13)

this subproblem has a closed-form solution as

${({\text{𝐰}}^i)}^{t+1}=\{\begin{array}{c}{({\text{𝐰}}^i)}^t-{\displaystyle \frac{1}{L}}\nabla f({({\text{𝐰}}^i)}^t),if{||{({\text{𝐰}}^i)}^t-{\displaystyle \frac{\nabla f({({\text{𝐰}}^i)}^t)}{L}}||}_2^2>{\displaystyle \frac{2\lambda }{L}}\hfill \\ \text{𝟎},\text{otherwise}\hfill \end{array}$(14)

in which the parameter L needs to be tuned. The upper bound on Lipschitz constant L f is not be easily calculated, thus we use the line search method to tune L as suggested in [37] until the objective value descent.

Homotopy Strategy: many works [39, 40] have verified that the sparse coding approaches benefit from a good starting point. Thus we can use the solution of Eq. (8), for a given value of λ , to initialize IHT in a nearby value of λ . Usually, the next loop with warm starting will require fewer iterations than the current loop. Using this warm-starting technique, we can efficiently solve Eq. (8) for a sequence of values of λ , which is called homotopy strategy. The complete routine of the HIHT algorithm for solving Eq. (8) is described in Algorithm 3.2.

For a fixed λ t , since ▽ f is Lipschitz continuous with constant L f , we have

f ⁢ ( 𝐖 k + 1 ) ⩽ f ⁢ ( 𝐖 k ) + t ⁢ r ⁢ ( ▽ f ⁢ ( 𝐖 k ) T ⁢ ( 𝐖 k + 1 - 𝐖 k ) ) + L f 2 ⁢ || 𝐖 k + 1 - 𝐖 k || 2 .(15)

Using this inequality, the fact that $L>L_f$ , and Eq. (11), we obtain that

φ λ t ⁢ ( 𝐖 k + 1 ) = f ⁢ ( 𝐖 k + 1 ) + λ t ⁢ || 𝐖 k + 1 || 0 ⩽ f ⁢ ( 𝐖 k ) + t ⁢ r ⁢ ( ▽ f ⁢ ( 𝐖 k ) T ⁢ ( 𝐖 k + 1 - 𝐖 k ) ) + L f 2 ⁢ || 𝐖 k + 1 - 𝐖 k || 2 + λ t ⁢ || 𝐖 k + 1 || 0 (16) ⩽ f ⁢ ( 𝐖 k ) + t ⁢ r ⁢ ( ▽ f ⁢ ( 𝐖 k ) T ⁢ ( 𝐖 k + 1 - 𝐖 k ) ) + L 2 ⁢ || 𝐖 k + 1 - 𝐖 k || 2 + λ t ⁢ || 𝐖 k + 1 || 0 ⩽ f ⁢ ( 𝐖 k ) + λ t ⁢ || 𝐖 k || 0 = φ λ t ⁢ ( 𝐖 k ) ,

where the last inequality follows from Eq. (11).1

φ λ t ⁢ ( 𝐖 k ) - φ λ t ⁢ ( 𝐖 k + 1 ) ⩾ L - L f 2 ⁢ || 𝐖 k + 1 - 𝐖 k || 2 .(17)

Since f ⁢ ( 𝐖 ) is bounded below, it then follows that φ λ t ⁢ { 𝐖 k } is bounded below. Hence, φ λ t ⁢ { 𝐖 k } converges to a finite value as k → ∞ and a local optimal solution 𝐖 λ t * can be achieved.

Since the λ is monotonically decreased, and 𝐖 λ t * is set as the initial solution for HIHT in λ t + 1 , we obtain that:

${\phi }_{{\lambda }_t}({\text{𝐖}}_{{\lambda }_t}^{*})>{\phi }_{{\lambda }_{t+1}}({\text{𝐖}}_{{\lambda }_t}^{*})={\phi }_{{\lambda }_{t+1}}({\text{𝐖}}_{{\lambda }_{t+1}}^0)⩾{\phi }_{{\lambda }_{t+1}}({\text{𝐖}}_{{\lambda }_{t+1}}^{*}),$(18)

it implies that the objective value is monotonically decreasing and a local optimal solution can be achieved by the proposed algorithm. We show an example of the objective value at each iteration in Fig. 1, it can be seen that the objective function values monotonically decrease at each iteration until convergence, which verifies the convergence of Algorithm 3.2 experimentally.

OPEN IN VIEWER

Inspired by [39], we derived an accelerated version for HIHT to reduce the computational time. For each fixed value of λ , Algorithm 3.2 repeats steps 6–14 to get a exact solution of Eq. (8). In the following, for acceleration of Algorithm 3.2, just one outer loop is called to replace steps 6–14. The complete routine of the the AHIHT method is described in Algorithm 3.4.

[t] Acceleration of HIHT(Input:) Training data 𝐗 ∈ R d × N , training labels 𝐘 ∈ R C × N , centering matrix 𝐇 ∈ R N × N , parameters L 0 , λ 0 , λ 𝑡𝑔𝑡 ; (Output:) 𝐖 * ; 1: initialize ${\text{𝐖}}^0=\text{𝟎},\rho \in (0,1),ϵ>0,\eta >0,\gamma >1,k\leftarrow 0$ ; 2: do An L-tuning iteration3: update 𝐖 k + 1 by Eq. (14);4: while do5: L k ← γ ⁢ L k ;6: update 𝐖 k + 1 by Eq. (14);7: end while8: L k + 1 ← L k ;9: λ k + 1 ← ρ ⁢ λ k ;10: k ← k + 1 ;11: until λ k + 1 ⩽ λ 𝑡𝑔𝑡 12: 𝐖 ⋆ ← 𝐖 k .

4.1 Data sets description

Eight biological benchmark data sets are used to validate the performance of our approach in the experiments: Brain [41], Breast3 [42], Leukemia [43], Lung [44], Lymphoma [45], NCI [46], Prostate [47], and Srbct [48]. The 8 benchmark data sets were downloaded from Feiping Nie’s homepage,2

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

OPEN IN VIEWER

In the experiments, we evaluate the feature selection performance on classification accuracy, where two popular classifiers KNN and Softmax are used. For each data set, 2 3 of samples per class are randomly selected as training samples, and the rest samples are responsible for testing. Ten repeated trials are carried out and the average results are recorded for comparison. We compare our feature selection method with Baseline (without feature selection) and six state-of-the-art feature selection methods: (1) two basic filter methods: Relief [10], and mRMR [12]; (2) two l 2 , 1 -norm based methods: RFS [18], and RLSR [21]; (3) two l 2 , 0 -norm constrained methods: RPMFS [24], and EFSF [25]. The codes of Relief and mRMR are provided by the FEAST package [49], and others can be downloaded from the authors’ homepages.

For our method, the parameter λ is searched in the grid of { 10 - 5 , 10 - 4 , … , 10 0 } , and for RFS and RLSR it is searched in { 10 - 5 , 10 - 3 , … , 10 5 } . For all methods, the number of selected features is tuned from { 20 , 40 , … , 400 } in this paper. The k nearest neighbors of KNN is set as k = 5 . All other parameters take the default values as suggested by the authors. For RLSR, all training data are used as labeled data, thus it can be seen as a supervised method here. The number of selected features with highest classification accuracy are recorded.

4.3 Classification performance

The classification performance of our approach is evaluated in this part. Table 2 shows the best Acc of each method with corresponding No.fea. From this table it can be seen that, our method can achieve highest Acc with fewest No.fea on most data sets, or comparable results to the best ones. On NCI dataset, the Acc of our approach obtained by Softmax is 78.33%, which is 3.33% higher than the second one (expect Baseline). RPMFS, ESFS and our algorithm both solve the original l 2 , 0 -norm problem, but RPMFS and ESFS try to solve an equality constraint problem so that the number of selected features needs to be tuned carefully to obtain a satisfactory classification result. The numerical comparison shows that our method outperforms RPMFS and ESFS most of the time, which indicates that our approach is more efficient than RPMFS and ESFS. The classification results demonstrate that our method can remove more redundant features while maintains the discriminative performance.

Figures 2 and 3 show Acc V.S. No.fea obtained by KNN and Softmax, respectively. It is obvious that the proposed method distinctly outperforms other approaches in the most experimental data sets. What’s more, compared with other lines, ours are not fluctuated too much as the No.fea changed, which indicates the performance of the proposed approach is much more robust to the number of selected features than other methods.

4.4 Parameter sensitivity

We evaluate the sensitivity of λ in this part, the Acc is employed to evaluate the performance of classification with λ searched in the grid of { 10 - 5 , 10 - 4 , 10 - 3 , 10 - 2 , 10 - 1 } . The number of selected features varies in { 30 , 60 , … , 150 , 180 } . The data sets of Brain, Leukemia, Lung and Prostate are used for testing, and the experimental results are shown in Fig. 4. It can be seen from this figure that, the performance of l 2 , 0 -FS is not very sensitive to λ . When λ is in the range of { 10 - 5 , … , 10 - 2 } , there is almost no fluctuation on Acc. Thus the proposed approach is robust to λ and there is no need to spend much time to tune the value of λ . It also can be seen that the performance gets worst when λ = 10 - 1 , the reason is that a large value of regularization factor will hurt data approximation performance.

L 2 , 0 -norm is a non-convex problem, which may make the solution sensitive to the initialization. Three kinds of initialization are used to explore the effect on classification performance: zero-valued initialization, random Gaussian distributed initialization, and random uniform distributed initialization. For Gaussian and uniform distributed initialization, ten repeated trials are carried out and the average results are recorded. Figure 5 shows the results (in this figure, the No.fea is equal to the number of non-zero rows of W), in which data sets Breast3, Lung, and Srbct are used. It can be seen that these three different initialization methods can get the same results, indicating that AHIHT is not sensitive to initialization when applied to feature selection.

4.5 Comparison of convergence speed and time consumption

Figure 6 plots the objective function value for each iteration of the three l 2 , 0 -norm based methods. As can be observed, our method can decrease the objective value quickly in early iterations, which indicates that our method converges faster than RPMFS and ESFS. Figure 7 shows the average computational time of each sparsity regularization method on the eight data sets, which also compares AHIHT with HIHT. It can be seen that the computational time of RFS, RPMFS and AHIHT are much less than the other three methods, the reason is that the computational complexity of RFS, RPMFS and AHIHT is in propotion to d while for other methods it is in propotion to d 2 or d 3 . Comparing AHIHT with HIHT, it can be found that AHIHT can reduce the computational time of HIHT effectively, thus is more suitable for practical applications.

4.6 Comparison of AHIHT and HIHT

In this part, we compare the performance of HIHT and its accelerated version in terms of accuracy and the row-sparsity of the produced solutions with different values of λ . Figure 8 shows the results, in which data sets Breast3, Lung, and NCI are used. From this figure, it can be seen that, when λ is less than 0.1, HIHT can get a more sparse solution than AHIHT, while AHIHT can get higher classification accuracy than HIHT, it means that the features selected by AHIHT are more useful for classification than those of HIHT. What’s more, the results obtained by AHIHT demonstrate AHIHT is more robust to λ than HIHT. For HIHT, it should tune a small value of λ to get satisfactory results, while a small value of λ will spend more computational time, which has been shown in Fig. 7. In conclusion, AHIHT is more practical than HIHT for feature selection.