Robust Flexible Feature Selection via Exclusive L21 Regularization

1 Introduction

Feature selection plays an important role in many machine learning tasks. The main purpose is to remove irrelevant and redundant noise features in high-dimensional data space. The selected features will help to reduce the computation cost and improve the performance on real-world applications.

There are many research works on feature selection over the years. Generally, feature selection methods can be divided into three main categories [Guyon and Elisseeff, 2003]: wrapper method, filter method, and sparse coding based method (also known as embedded method). The most representative wrapper method is support vector machine recursive feature elimination (SVM-RFE) [Guyon et al., 2002], but the computation cost is extremely high. Contrarily, filter method is very efficient such as F-statistic [Ding and Peng, 2003], ReliefF [Robnik-ˇSikonja and Kononenko, 2003], minimum redundancy maximum relevance (mRMR) [Peng et al., 2005].

Recently, sparse coding based methods have been widely investigated, and applied to the study of feature selections.

Least absolute shrinkage and selection operator (LASSO) [Tibshirani, 1996] is a regression based analysis method that incurs the sparsity on weights via ℓ1-norm. ℓ1-SVM [Zhu et al., 2003] and hybrid huberized SVM (HHSVM) [Wang et al., 2007] are introduced to further improve performance on two-class problem. LASSO can be derived from probabilistic selection on ridge regression [Ming et al., 2019].

To solve multi-class problem, researchers search a subset of features shared by all the classes, also known as multi-task feature learning (MTFL). In this area, ℓ2,1-norm is the most widely used regularization developed in [Liu et al., 2009; Nie et al., 2010; Gui et al., 2017]. In [Quattoni et al., 2009], authors propose ℓ1,∞-norm regularization, which shares same property of row-sparsity as ℓ2,1-norm. As compared to classshared feature selection, exclusive lasso (eLASSO) [Zhou et al., 2010; Campbell and Allen, 2017] proposes to capture the negative correlation among different classes via ℓ1,2-norm, which is first introduced in [Zhao et al., 2009] called composite absolute penalties (CAP). In exclusive feature learning, discriminative features are selected for each class to provide certain flexibility. Based on this, Kong et al. of [Kong et al., 2014] propose to solve the mix of ℓ1-norm and ℓ1,2-norm, towards minimizing the feature correlation.

Motivated by previous works, in this paper, we introduce a novel regularization called "exclusive ℓ2,1", which is short for "ℓ2,1 with exclusive lasso". The exclusive ℓ2,1 regularization brings out joint sparsity at inter-group level and exclusive sparsity at intra-group level simultaneously. Thus, the proposed regularization can combine the advantages from different sparsity-induced terms, which not only removes irrelevant noise features (i.e. increase the robustness via ℓ2,1-norm) but also selects discriminative features for each class (i.e. provide the flexibility via ℓ1,2-norm).

The main contribution of this paper includes: (i) a novel "exclusive ℓ2,1" regularization is proposed to conduct robust flexible feature selection; (ii) we point out some interesting properties of ‖w‖21 regularization as compared to ‖w‖1 regularization; (iii) a sorting based explicit approach is introduced to directly solve the ℓ1,2-norm regularization; (iv) an efficient augmented Lagrange multipliers (ALM) based optimization algorithm is proposed to iteratively solve the "exclusive ℓ2,1" regularization in a row-wise fashion; (v) experimental results on twelve benchmark datasets demonstrate that the proposed regularization outperforms state-of-the-arts.

2 Notations and Definitions

Throughout this paper, scalars, vectors, and matrices are denoted as lower-case/capital letters, boldface lower-case letters, and boldface capital letters, respectively.

The i-th element of vector w is represented by wi. Given a matrix W = (Wij) ∈ Rd×k , the i-th row is represented by wi (i.e. W = [w1 ; · · ·;wd ]), and the j-th column is represented by wj (i.e. W = [w1, · · ·, wk]). The Frobenius norm

X = [x1, · · ·, xn] ∈ Rd×n represents n data points, where xi ∈ Rd , and corresponding class labels are defined as Y = [y1 ; · · ·; yn ] ∈ Rn×k , where yi ∈ Rk is one-hot vector and yij = 1 or Yij = 1 means i-th sample belonging to j-th class.

3 Exclusive ℓ2,1 Regularization

Generally, sparse coding based methods can be formulated as minW {f(W)+λΩ(W)}, where f(W) is the loss function, Ω(W) is the regularization, and λ is the hyperparameter.

Our work is motivated from the following observations. The ℓ2,1 norm based feature selection (i.e. f(W)+λ‖W‖2,1) incurs joint sparsity on rows. A selected non-zero row could still have some elements with small (in magnitude) numerical values. Suppose one of them is Wij. This implies i-feature is not highly correlated with j-th class. Thus ℓ2,1 alone is too rigid for feature selection.

On the other end, exclusive lasso (i.e. f(W) +λ‖W‖21,2 ) selects discriminative features for each class. Here, as λ increases, different elements in squared ℓ1-norm of i-th row wi are competing with each other to survive. Thus, at least one element in row wi survive (remaining non-zero). The problem with exclusive lasso in this context is: all features/rows will be selected, because for each feature/row i, there will be some non-zero elements even at large regularization strength.

Towards resolving above main concerns for using ℓ2,1 regularization alone or using exclusive lasso alone, we propose to combine them together as a new regularization defined as Ω(W) = α ‖W‖2,1 + β ‖W‖21,2 , which will be called "exclusive ℓ2,1" short for "ℓ2,1 with exclusive lasso". As a result, ℓ2,1-normwill increase the robustness to help ℓ1,2-norm, and ℓ1,2-norm will provide the flexibility to help ℓ2,1-norm.

3.1 An Illustration

The synthetic data X, Y is given in Eq. (1), where d = 7, n = 8, k = 3. The loss function is f(W) = ‖XT W- Y‖2F . The learned matrices are given in Eq. (2), where the number of non-zero elements in W is enforced to 12 for each regularization. The difference is explained as follows:

• (i) W21 (ℓ2,1): a feature can be selected by all the classes (e.g. 3rd row is selected for 1st, 2nd, 3rd classes), or can be discarded (e.g. 4th row is a zero vector).
• (ii) W12 (exclusive lasso): a feature can be selected by some classes (e.g. 5th row is selected for 1st, 3rd classes; 6th row is selected only for 3rd class), but can not be discarded since the matrix has no zero rows.
• (iii) Wex21 (the proposed "exclusive ℓ2,1"): a feature can be selected by all the classes (e.g. 2nd row is selected for 1st, 2nd, 3rd classes), or can be selected by some classes (e.g. 3rd row is selected for 2nd, 3rd classes), or can be discarded (e.g. 4th row is a zero vector).

4 Understanding the Exclusive Sparsity

4.1 Interesting Property of ‖w‖21 Regularization

In this paper we use ‖w‖21 regularization for flexible feature selection. Here, we point out some interesting properties of this regularization.

Consider ‖w‖21 regularization first. We investigate the following simple proximal operator-type problem:

This is very similar to the standard ℓ1-norm regularization problem

which has been thoroughly studied in connection to lasso [Tibshirani, 1996].

There exists a widely held belief that optimization problems Eq. (3) and Eq. (4) behave very similarly and their solutions have identical sparsity pattern.

This belief comes from the following reasoning. Problem (3) is equivalent to

for some parameter t. And problem (4) is equivalent to

for some parameter t.

However, this widely held belief is incorrect.

Let w∗ ℓ12 be the optimal solution for problem (3). Let w∗ ℓ1 be the optimal solution for problem (4). We illustrate their significant differences in two simple cases.

Case 1 is a simple problem in 2-dim. a = (2, 1). Optimal solutions are (computed using algorithm explained later 1 )

Clearly as λ increases above 1, w∗ ℓ1 is all zeros, but w∗ ℓ12 has non-zero component.

Case 2. Consider the dimension is one with a = 1. These problems can be solved analytically. The solutions are

Clearly, when λ > 2, w∗ ℓ1 = 0, but w∗ ℓ12 is never zero not matter how large λ is.

These two cases show that as λ increases to large values, w∗ ℓ1 will become exact zero for all components, while w∗ ℓ12 will become zero for d - 1 components and one component approaches 1 1+λ asymptotically.

4.2 Solving ℓ1,2-NormRegularization

In [Zhou et al., 2010], authors illustrate the sparsity of ℓ1,2norm from a projection point of view, then solve a min-max optimization problem. Kong et al of [Kong et al., 2014] use a re-weight strategy to solve ℓ1,2-norm regularization.

However, both methods are inefficient in high-dimensional data space. Inspired by non-negative shrinkage thresholding operator [Cavazza et al., 2018], we introduce a sorting based explicit approach to solve ℓ1,2-norm regularization. Here, we focus on its simplified formulation in Eq. (3), which then can be applied to solve multi-class problem in section 5.

Lemma 1. The optimal solution w∗ of Eq. (3) has the following property of its sign: for i = 1, · · ·, d, (i) if ai = 0, w∗ i = 0; (ii) if ai ̸= 0, sign(w∗ i ) = sign(ai).

Proof of Lemma 1. If ai = 0, w∗ i = 0 can be easily verified. If ai ̸= 0, suppose w∗ i = c and sign(c) ̸= sign(ai). However, w∗ i = -c gives the lower objective value, since |c| = | - c| and (c - ai)2 > (-c - ai)2 . Thus, sign(w∗ i ) = sign(ai). □

1 For standard ‖w‖1 regularization, Eq. (4) has the closed-form solution as w∗ ℓ1 = sign(a)⊙[|a|-λ/2]+. For ‖w‖21 regularization, we propose a sorting based explicit approach (see Theorem 5) to solve Eq. (3), and the optimal solution w∗ ℓ12 is given in Eq. (10).

Lemma 2. The optimal solution w∗ of Eq. (3) has the following property of its magnitude: for i = 1, · · ·, d,

where ξi is the subgradient of f(x) = |x|, x ≥ 0 at x = 0.

Proof of Lemma 2. Eq. (3) can be rewritten equivalently as

since [sign(wi)]2 = [sign(ai)]2 = 1, according to Lemma 1. Taking derivative of J(w) in Eq. (9) w.r.t |wi| and setting ∂J(w) ∂|wi| = 0, we will have the same first-order optimality conditions defined in Eq. (7) and Eq. (8). □

Proposition 3. As λ increases to large values, at least one element wi in w will survive (i.e. |wi| > 0), given a ̸= 0. Otherwise, w = 0 will lead to a = 0 according to Eq. (8).

Definition 4. Given a = (a1, · · ·, ad) ∈ Rd , S denotes a ddimensional vector with Si ̸=Sj(i ̸=j), ⋃di=1Si ={1, · · ·, d}, and each Si represents the indexes of a descending order with respect to a, such as |aS1 | ≥ |aS2 | ≥ · · · ≥ |aSd |.

Theorem 5. The optimal solution of Eq. (3) is given by

where ⊙ is the Hadamard product, i.e. [x ⊙ y]i = xiyi, [·]+ = max(·, 0), µτ = 1 τ ∑τi=1 |aSi |, and τ is the largest coordinate of S satisfying |aSτ | - λτ 1+λτ µτ > 0.

Proof of Theorem 5. Suppose that w∗S1 , w∗S2 , · · · , w∗Sτ are non-zeros. By adding Eq. (7) for S1, S2, · · · , Sτ (i.e. the first τ indexes saved in S), we have

which can be equivalently rewritten as ‖w∗ ‖1 = τ 1+λτ µτ , where µτ = 1 τ ∑τi=1 |aSi |. Thus, Lemma 2 and Eq. (11) give the optimal solution w∗ w.r.t its magnitude as follows

which is equivalent to the definition of w∗ in Eq. (10), since w∗ j = sign(w∗ j )|w∗ j | = sign(aj)|w∗ j | for j = S1, · · ·, Sτ , and w∗ j = 0, |aj| - λτ 1+λτ µτ < 0 for j = Sτ+1, · · ·, Sd. □

Theorem 6. When τ is the largest coordinate of S satisfying |aSτ | - λτ 1+λτ µτ > 0, the solution w∗ defined in Eq. (10) achieves the global minimum of J(w).

Proof of Theorem 6. If τ = d, we have |w∗Si | > 0 for i = 1, · · ·, d, and each w∗Si given by Eq. (12) satisfies optimal condition Eq. (7). Thus, w∗ is the global minimizer of J(w).

If τ < d, we have |w∗Si | > 0 for i = 1, · · ·, τ , and |w∗Si | = 0 for i = τ+1, · · ·, d. Since τ is the largest coordinate of S

Algorithm 1 Search the largest coordinate τ of S.

Input: a ∈ Rd , λ, S. Output: τ , µτ .

• 1: Initialize: τ = d, µτ = 1d ∑di=1 |aSi |. 2: while τ > 1 and |aSτ | - λτ 1+λτ µτ < 0 do 3: µτ = τ τ-1 µτ - 1 τ-1 |aSτ |. 4: τ = τ -1.
• 5: end while
• 6: return τ, µτ .

satisfying |aSτ |- λτ 1+λτ µτ > 0, for (τ+1)-th coordinate of S, we have |aSτ+1 |- λ(τ+1) 1+λ(τ+1) µτ+1 < 0, which can be rewritten equivalently as |aSτ+1 |- λτ 1+λτ µτ < 0, i.e. |aSτ+1 |<λ‖w∗ ‖1. This implies w∗Sτ+1 satisfies Eq. (8). For i = τ+2, · · ·, d, w∗Si satisfies Eq. (8), since |aSi | ≤ |aSτ+1 | < λ‖w∗ ‖1. Besides, w∗Si satisfies Eq. (7) for i = 1, · · ·, τ. Thus, w∗ is the global minimizer of J(w), which completes the proof. □

Since w∗ depends on τ , µτ , here we introduce an efficient algorithm (given in Algorithm 1) to search the largest coordinate τ of S satisfying |aSτ | - λτ 1+λτ µτ > 0 in linear time.

5 Optimization Algorithm

To select robust and flexible features, we are interested in the following optimization problem

where the least square loss is penalized by the proposed "exclusive ℓ2,1" regularization, and α, β are hyperparameters.

First, we add an auxiliary variable Z to make the optimization separable between ℓ2,1-norm and ℓ1,2-norm. Thus, original problem (14) becomes

Then, augmented Lagrange multipliers (ALM) method is applied to enforce the constraint in problem (15) explicitly

where 〈·, ·〉 is the inner product, i.e. 〈A,B〉 = ∑ij AijBij, Λ is the Lagrange multiplier, and ν is the penalty parameter. Problem (16) can be rewritten equivalently as

Thus, our task is to solve the variables Z, W and update the parameters Λ, ν.

5.1 Solving for Z

Firstly, we solve Z while fixing W. Then, problem (17) w.r.t Z becomes

Since the optimizations of each row zi in Z are separable, we can minimize problem (18) in a row-wise fashion. Thus, the optimization of Eq. (18) w.r.t zi becomes

where i = 1, · · ·, d is the feature/row index, e = wit -λit /νt, wit is the i-th row ofWt, and λit is the i-th row of Λt. Using Theorem 5, the optimal solution of Eq. (19) is

where τ , µτ are computed using Algorithm 1, given the input (e, 2β/νt, S), and S is a k-dimensional vector representing the indexes of descending order |eS1 | ≥ |eS2 | ≥ · · · ≥ |eSk |.

5.2 Solving forW

Secondly, we solve W while fixing Z. Then, problem (17) w.r.t W becomes

Since ℓ2,1-norm is defined on each row wi inW, here we can solve W in the similar way as Z.

To solve W in a row-wise fashion, we decompose the least square loss w.r.t wi as follows

where Y-i = Y - ∑j̸=i (xj )T wj , a = ‖xi ‖22 , b = xiY-i , and c = Tr((Y-i )T Y-i ).

Thus, the optimization of Eq. (21) w.r.t wi becomes

where i = 1, · · ·, d is the feature/row index, d = 1 2a+νt (2b+ νtzit+1 +λit ), zit+1 is the i-th row of Zt+1, and λit is the i-th row of Λt.

The optimal solution of Eq. (23) is given by 2

where [·]+ = max(·, 0).

5.3 Updating Parameters

Finally, we update parameters Λ, ν at the end of t-th iteration as the following

where ρ > 1 is a constant.

2 For standard ‖w‖2 regularization, the optimization problem minw {‖w - a‖22 +λ‖w‖2 } has the closed-form solution as w∗ ℓ2 = max(1 - λ 2‖a‖2 , 0)a.

Algorithm 2 ALM based optimization algorithm for solving the "exclusive ℓ2,1" regularization in problem (14).

Input: data matrix X ∈ Rd×n , class labels Y ∈ Rn×k , hyperparameters α, β. Output: weight matrix W ∈ Rd×k . 1: Initialize: t = 0, νt = 1/‖X‖F , ρ = 1.1, ϵ1 = 1e-8, ϵ2 = 1e-5, Λt = 0, random initialization weights Wt. 2: repeat 3: for i ∈ {1, · · · , d} do 4: Compute e via Eq. (19). 5: Compute the descending order S of |e1|, · · ·, |ek|. 6: Compute τ, µτ via Algorithm 1 given (e,2β/νt,S). 7: Compute zit+1 via Eq. (20). 8: end for 9: for i ∈ {1, · · · , d} do 10: Compute Y-i , a, b via Eq. (22). 11: Compute d via Eq. (23). 12: Compute wit+1 via Eq. (24). 13: end for 14: Update Λt+1 via Eq. (25). 15: Update νt+1 via Eq. (26). 16: Set t = t +1. 17: until convergence condition is satisfied: ∣∣Jex21(Wt+1 ) - Jex21(Wt )∣∣/Jex21(Wt ) ≤ ϵ1, ∥∥Zt+1 -Wt+1∥∥∞ ≤ ϵ2. 18: return the optimal solution: W∗ .

5.4 The Summary of Optimization Algorithm

The complete framework of the proposed augmented Lagrange multipliers (ALM) based optimization algorithm is summarized in Algorithm 2.

6 Experiments

6.1 Benchmark Datasets

Experiments on twelve benchmark datasets are conducted to evaluate the performance of feature selection methods on classification. Among those benchmarks, there are 4 image datasets: MNIST 3 [Lecun et al., 1998], Yale 4 , YaleB 5 , PIE [Sim et al., 2002]; 1 spoken letter recognition dataset: ISOLET 6 ; 5 bio-microarray datasets: Carcinomas [Yang et al., 2006], Lung [Bhattacharjee et al., 2001], Glioma [Nutt et al., 2003], TOX 6 , Tumor-14 [Ramaswamy et al., 2001]; and 2 text datasets: CNAE-9 [Ciarelli and Oliveira, 2009], 20Newsgroups 7 . The details of all the benchmark datasets is summarized in Table 1.

6.2 EvaluationMetrics

In subsequent experiments, the proposed exclusive ℓ2,1 regularization is compared to five state-of-the-arts, including three filter methods: F-statistic [Ding and Peng, 2003], ReliefF [Robnik-ˇSikonja and Kononenko, 2003], minimum redundancy maximum relevance (mRMR) [Peng et al., 2005], and two sparse coding based methods: multi-task feature learning via ℓ2,1-norm (ℓ2,1) [Liu et al., 2009; Nie et al., 2010; Gui et al., 2017], exclusive Lasso (eLASSO) [Zhou et al., 2010; Campbell and Allen, 2017].

3 In MNIST, 100 images are randomly selected out of each digit.

4 http://vision.ucsd.edu/content/yale-face-database

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

6 http://featureselection.asu.edu/datasets.php

7 http://qwone.com/∼jason/20Newsgroups/ In 20-Newsgroups,

100 documents are randomly selected out of each newsgroup, and F-statistic method is used to prescreen 5,000 keywords.

To evaluate the performance on classification, 5-fold crossvalidation accuracy with SVM as classifier are computed on average. LIBSVM [Chang and Lin, 2011] is used as the practical implementation of SVM, where kernel is set as linear and parameter C is set as 1 for all the experiments.

When training different models, hyperparameters are adjusted to enforce the same level of sparsity on learned weight matrices. Then we select the top features with largest weights for each class. For testing, an SVM classifier is built for each class separately, by using the selected features. The final classification result is obtained via majority voting.

6.3 Analysis of the Results

Convergence Study

Convergence of our proposed ALM based optimization algorithm is shown in Fig. 1, where x-axis and y-axis denote the number of iterations and the objective value respectively.

We use the same hyperparameter setting, i.e. α=1, β=1, for four benchmark datasets. As it can be seen in Fig. 1, the proposed optimization algorithm takes around 100∼150 iterations to converge. This shows our ALM based optimization algorithm is efficient and converge fast in real applications.

Classification Results Comparison

Experimental results of our proposed exclusive ℓ2,1 regularization versus five state-of-the-arts are shown in Fig. 2, where x-axis denotes the number of selected features ranging from 10 to 80, and y-axis denotes the average of 5-fold crossvalidation classification accuracy.

In general, sparse coding based methods (ℓ2,1, eLASSO, Ours) achieve better performances than filter methods (FStatistic, ReliefF, mRMR). Among filter methods, mRMR has relatively higher classification accuracy, since it takes consideration of minimizing the correlation between features.

eLASSO performs well in image and spoken letter recognition datasets. However, its performance has a great degradation in bio-microarray and text datasets, since ℓ1,2-norm cannot remove a large amount of irrelevant noise features in highdimensional data space. ℓ2,1 has a very stable performance in all datasets via selecting class-shared features. In some cases, ℓ2,1 performs even close to our method around top 60∼80 features. Overall, our method obtains the best classification result on twelve benchmark datasets. Additionally, in the small number of selected features setting, e.g. top 10∼20, our method has an overwhelming advantage over other methods, with around 5%∼10% improvement on accuracy.

7 Conclusion

In this paper, we introduce a novel "exclusive ℓ2,1" regularization for robust flexible feature selection. Besides, we point out some interesting property of ‖w‖21 regularization, which can be solved directly by a sorting based explicit approach. Then, an efficient augmented Lagrange multipliers based optimization algorithm is proposed to iteratively solve the "exclusive ℓ2,1" regularization in a row-wise fashion. Extensive experiments validate the effectiveness of the proposed robust flexible feature selection, which outperforms state-of-the-arts on twelve benchmark datasets.