Feature Selection using Stochastic Gates
Authors: Yutaro Yamada
Ofir Lindenbaum
Sahand Negahban
Yuval Kluger
Venue:  (2020)
Keywords: Feature selection
classification
regression
survival analysis

Abstract: Feature selection problems have been extensively studied in the setting of
linear estimation, for instance LASSO, but less emphasis has been placed on
feature selection for non-linear functions. In this study, we propose a method
for feature selection in high-dimensional non-linear function estimation
problems. The new procedure is based on directly penalizing the $\ell_0$ norm of
features, or the count of the number of selected features. Our $\ell_0$ based regularization relies on a continuous relaxation of the Bernoulli distribution, which
allows our model to learn the parameters of the approximate Bernoulli
distributions via gradient descent. The proposed framework simultaneously learns
a non-linear regression or classification function while selecting a small
subset of features. We provide an information-theoretic justification for
incorporating Bernoulli distribution into our approach. Furthermore, we evaluate
our method using synthetic and real-life data and demonstrate that our approach
outperforms other embedded methods in terms of predictive performance and feature selection.

Under review as a conference paper at ICLR 2020
FEATURE SELECTION USING STOCHASTIC GATES
Anonymous authors
Paper under double-blind review
ABSTRACT
Feature selection problems have been extensively studied in the setting of linear estimation, for instance LASSO, but less emphasis has been placed on feature selection for non-linear functions. In this study, we propose a method for feature selection in non-linear function estimation problems. The new procedure is based on directly penalizing the to norm of features, or the count of the number of se- lected features. Our &o based regularization relies on a continuous relaxation of the Bernoulli distribution, which allows our model to learn the parameters of the approximate Bernoulli distributions via gradient descent. The proposed framework simultaneously learns a non-linear regression or classification function while select- ing small subset of features. We provide an information-theoree incorporating Bernoulli distribution into our approach. Furthermore, evaluate our method using synthetic and real-life data and demonstrate that our approach outperforms other embedded methods in terms of predictive performance and feature selection.
justification for
INTRODUCTION
Feature selection is fundamental task machine learning and statistics. Feature selection leads to number of potential benefits: reducing experimental costs Min et al. (2144) enhancing interpretability (Ribeiro et al. 2) computational speed up and even improving model generalization on unseen data (Chandrashekar & Sahin 2014). In biomedicine, scientists collect multitude datasets comprising of many biomarkers (e.g., genes or proteins) that require development of effective diagnosties or prognostics models. For instance, in Genome wide association studies (GWAS), feature selection can help identify such models and lead to improved risk assessment and reduced cost.
Feature selection methods may be classified into three major categories: filter methods, wrapper methods, and embedded methods. Filter methods attempt to remove irrelevant features prior to learning model. These methods filter features using per-feature relevance score that is created
based on some statistical measure (Battiti 1994, Peng et al. 2005) Estévez et al. 2009; Song et al..
2007 2012) Chen et al.) 2017). Wrapper methods (Kohavi & John 1997b Stein et al. 2005 Zhu et al. 2007 Reunanen 2003] Allen 2013) use the outcome of a classifier to determine the relevance of each feature, which requires recomputing the classifier for each subset of features. This becomes computationally expensive for neural network based wrapper methods (Verikas & Bacauskiene/ 2002
Kabir et al. 2010) Roy et al. 2015).
Embedded methods aim to remove this burden by learning the model while simultancously selecting the subset of relevant features. The Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani 1996) is a well-known embedded method, whose objective is to minimize the loss while enforcing an &, constraint on the weights of the features. Although LASSO is scalable and widely used (Hans, 2009 Li et al. 2011 2006), it is restricted to the domain of linear functions. To allow the model to capture nonlinear interaction, appealing to consider the non-convex extension of the LASSO formulation using neural networks.
We develop a fully embedded feature selection method for nonlinear models. To the best of our knowledge it is the first to embedded feature selection method. Our method improves upon the LASSO formulation in two aspects: a) it captures nonlinear interactions between the features via neural network modeling and b) it employs an like regularization using gates whose weights parametrized by a smooth variant of a Bermoulli distribution. Altogether these twofold improvements formulated as fully differentiable neural network. Specifically, our contributions are follows:

---

Under review as a conference paper at ICLR 2020
1. By utilizing recent development of continuous and differentiable approximation discrete distributions (Maddison et al. 2016). (Jang et al. 2017). (Louizos et al. 2017). we introduce solution to the long standing problem of feature selection with an &o regularization.
2. We present novel simple relaxation of Bernoulli distribution to sparsify the imput layer (the feature space) which we call stochastic gate (STG) and show its advantage over the distribution presented by (Louizos et al., 2017) both in performance and convergence time.
3. By applying these two relaxations to input layer of neural network, perform embed- ded feature selection classification, regression or survival analysis tasks and demonstrate its capabilitiiss artificial and real data sets.
4. We justify our probabilistic approach by analyzing the constrained Mutual Information maximization objective of feature selection. We demonstrate the applicability of our method using numerous examples (see Section and Appendix).
Notation: We refer to vectors as bold lowercase and random vectors bold uppercase letters X. Scalars are non-bold case y, while random variables are capital case Y. A set represented by script fonts X, y. S. For example the ntth vector-valued observation is denoted as nn, whereas Xa represents the t feature of the vector-valued random variable X. Let [n] 2, 2. n,, For a set {0, 1) i be the characteristic function for the set. That is 8; 1 if S and otherwise. For two vectors æ and z we denote a @ z to be the element-wise product between æ and z. Thus, if we let E [0, be the characteristic vector of S, then we may define æş
S SC [D] let the vector
The & norm of vector denoted as ||>|||| Σp1 i-1 | Finally, the Eo norm of vector denoted as ||æ]|o and counts the total number of non-zero entries in the vector
2 PROBLEM SETUP AND BACKGROUND
Let X' RR be the input domain with corresponding response domain Y. Given realizations from some unknown data distribution Px,y, the goal of embedded feature selection methods is to simultaneously find subset of indices S (} and construct model that predicts Y based on the selected features Xs
2.1 RISK MINIMIZATION OBJECTIVE
We assume that we are given a family of functions F such that any function So E F is indexed by set of parameters Given some loss L, a selection of features C [D]. and a choice of parameters 0, we denote the risk of our model as
R(0 s) Ex ()),
where we recall that = {0. 1) is vector indicator variables for the set S and denotes the point-wise product. Thus, the goal of the feature selection problem is to find the parameters A and that minimize R(O,s) a) such that ||s|] ||||oo small compared to D.
2.2 FEATURE SELECTION FOR LINEAR MODELS
Before proceeding with our proposed method, we review the feature selection problem in the linear regression setting with the least squares loss. Thus, we restrict F to be the space of linear functions and the loss function to be the quadratic loss.
Given observations {æa, Da NN problem
we may consider the constrained empirical risk minimization
N Σ
min 0
L(0"
wn, Xn.ya) s.t. 10llo= S
Since the above problem intractable, a number authors replace the lo constraint with a surrogate function 22(0): E + R. designed to penalize the number of selected features in 9. A natural choice for S2 is the E, norm, which yields convex problem and more precisely the LASSO optimization problem. The & is known to be the closest convex relaxation to the fg. In fact, in certain settings the
2

---

Under review as a conference paper at ICLR 2020
LASSO optimization and the &o based objective have the same solution. While the original LASSO problem focuses on the constrained optimization problem, the regularized least squares problem, which is often used in practice, yields the following minimization objective:
N
1
min
E0"',, x, - yn)) A||0|11.
(3)
1=1
The hyperparameter À trades off the amount of regularization versus the fit of the objective The & regularized method is very effective for feature selection and prediction; however, it achieves this through shrinkage of the coefficients. As a result, (Fan Li 2001) have considered non-convex choices for f that perform well both theoretically and empirically for prediction and feature selection.
Our goal to apply such regularization techniques to perform feature selection while learning non-linear function. Kernel methods have been considered (Yamada et al.) (014), but scale quadratically in the number of observations. To alleviate this burden, (Gregorová et al. (2018) use random Fourier features approximate the kernel. (Li et al./ (016) and (Scardapane et al. 2017)) take an alternative approach by modeling fg using a neural network with & regularization on the input weights. However, practice, introducing an & penalty into gradient descent does not provide sufficient sparsification. Below, we discuss our method that works to directly use an Ep penalty.
3 PROPOSED METHOD
We take a probabilistic approach to approximate the lo norm, which can extend non-linear models while remaining computationally efficient. To motivate such probabilistic approach, we provide theoretical support (see Section 4) based on a Mutual Information perspective of the feature selection problem.
To view the la regularized version of the risk (Fq. I from a probabilistic perspective, one can introduce Bernoulli random vector Š whose entries are independent and the ath entry satisfies Ty = P(35 = 1) for d E [D]. If we denote the empirical expectation over our observations as Ex,y, then, the empirical regularized risk (Eq. becomes
min (0, π = min Ex Eg 0,
L(fa(X O 5), Y) + X||
(4)
0,
where we have Eg 5/% = Σt Ta and we constrain na E {0, 1]. Clearly, this formulation is equivalent to Eq. [ with a regularized penalty on cardinality rather than an explicit constraint. We may then relax the discrete constraint on Fy to Fy E (0, 1].
Now, our goal is to find the model parameters 0" and Bernoulli parameters " that minimize the empirical risk R((0, π)) via gradient descent. However, an optimization of loss function which includes discrete random variables suffers from high variance Fsn the Appendix for more details). Therefore, inspired by recently developed continuous approximation for discrete random variables, suggested by (Jang et al. 2017 Maddison et al.) 2016), we develop and use a novel and simple continuous distribution that is fully differentiable and suited for the task of feature selection.
3.1 CONTINUOUS RELAXATION
Our continuous relaxation for the Bernoulli variables 5. for E [D] is termed stochastic gate (STG). The STG relies on the reparametrization trick, which is widely used for reducing the variance of gradient estimators (Miller et al., [017]: Figurnov et al. 2018). To construct a con- tinuous approximation to Bernoulli random variable via the reparametrization trick, we define zd g(ss + ea) = max(0, min(1, eu + pa 0.5)) where ed is drawn from Gaussian distribution (0 oa), where σa fixed throughout training. This approximation can be viewed clipped, mean-shifted, Gaussian random vector. Furthermore, the gradient of the objective with respect to pa can be computed via the chain rule.
We can now rewrite the objective in Eq. 4 as
min A(e,, 11) = min E Ex E, (X O Z), Y) + A Z101 0, μ
(5)
0, μ
'A has one-to-one correspondence to à in the convex setting via Lagrangian duality.
3

---

Under review as a conference paper at ICLR 2020
where Z is random vector with D independent variables 2d for [D]. To optimize the empirical surrogate of the objective (Eq. [5), we first differentiate it with respect to μ. Then, Monte Carlo sampling leads us to the following gradient estimator
1 Σ L' R(O, p1) - K
Ozh
O
Σ
Pr Pr(za 0}.
dpd
Opd
K Opd
where K is the number of Monte Carlo samples. Thus, we can update the parameters jaa for [D] 0= via gradient descent. We note that if replace Bug with the above gradient estimator for L' is reduced to the Straight-Through estimator (Bengio al. 2013).
Under the continuous relaxation, the expected regularization term in the objective R((O, (t) is simply
the sum of the probability that the gates {=d (2 active, which is equal to E + Ma Ha+2 du
where the standard Gaussian CDF To conclude, we can now optimize the objective in Eq. using gradient descent over the model parameters and the parameters μ representing the Gaussian' mean (instead of the Bernoulli parameters πT).
After training, to remove the stochasticity from the learned gates, we set بي max (0, min(1, pd 0.5)), which informs what features are selected. Note that when |pa| is less than 1. 2 returns the value between (0, 1). such a case, we can treat the value of & as feature importance or employ an additional thresholding (i.e. 1 if &a >> 0.5 and 0 otherwise) depending on application-specfff needs. In the Appendix, we provide the pseudo-code of our algorithm as well as the discussion of the choice of d.
4 CONNECTION TO MUTUAL INFORMATION
In this section we show an equivalence between the Bernoulli formulation of the feature selection problem and the to regularized approach.
4.1 MUTUAL INFORMATION BASED ORJECTIVE
From an information theoretic perspective, the goal of feature selection is to find subset of features S that has the highest Mutual Information (MT) with the target variable Y. Recall that the MI between two random variables can be defined as T(X: Y) = H(Y) H where HY)), H are the entropy of pYY and the conditional entropy of Py[x (Y) respectively (Cover & Thomas 2006). Then we can formulate the task as selecting S such that the mutual informatron between X and Y are maximized:
max [(X Y)
s.t. |S| = k,
where & the hypothesised number of relevant features.
4.2 INTRODUCING RANDOMNESS
We first demonstrate that under mild assumptions we can replace the deterministic search over the set S (or corresponding indicator vector 8), by a search over the parameters of the distributions that model Our proposition is based on the following two assumptions:
There exists a subset indices S* with cardinality equal to such that for any i S* we have I(X;; x|X),)) (1)) 0.
Assumption 1:
Assumption 2: T(X x 0.
Discussion of assumptions: Assumption that including element from S improves prediction accuracy. This assumption is equivalent to stating that feature is strongly relevant (Kohavi & John 1997a Brown et al. 2012). Assumption 2 simply states that S* is a superset of the Markov Blanket of the variable Y (Brown et al. 2012). The assumptions are quite benign. For instance they are satisfied if X is drawn from a Gaussian with a non-degenerate covariance matrix and Y f(Xs-) w, where w is noise independent of and f is not degenerate. With these assumptions in hand, we may present our result.
4

---

Under review as a conference paper at ICLR 2020
Proposition 1. Suppose that the above assumptions hold for r the model. Then, solving the optimization in Eq. is equivalent to solving the optimization
[(X ④ S.Y) s.f. ΣEIS, S k,
max OTT=
where the coordinates St are drawn independently at random according to Bernoulli distribution with parameter T.
Due to length constraints, we leave the proof of this proposition and how bridges the MI maximiza- tion (Eq. 6 and risk minimization (Eq. to the Appendix (see Section C and D
5 RELATED WORK
The two most related works to this study are (Louizos et al. 2017) and (Chen et al. 2018). In (Louizos et al. 2017). they introduce the Hard-Concrete (HC) distribution as a continuous surrogate for Bernoulli distributions in the context of model compression. The HC distribution is induced by applying a hard sigmoid 2 = min(1, (1, max(0, $)) to 8, where we construct s by applying Sigmoid function to the logistic distribution follows:
U(,,1), log(u) log(1 (1- u), 8 + exp(
(lug (luga+))) s(5 - +
The interval (T, (). with 7 0 and Ç > allows the distribution have more probability mass on the edge of the support. The authors demonstrate how incorporating the HC in deep neural networks leads to fast convergence and improved generalization due to the sparsification effect. this study, we introduce this type of sparsification for the long standing problem of feature selection with to regularization. We demonstrate that a simple relaxation of Bernoulli distributions sufficient and works better than the HC distribution for feature selection. We find empirically that the HC yields high variance gradient estimates compared to the STG. This high-variance does not affect model sparsification as the sparsity pattern does not matter. Contrarily, for feature selection it is crucial as practitioners require stable procedure. Furthermore, higher variance gradients will result in slower convergence of the model (this has been demonstrated empirically in the supp. material). In (Chen) ct al. 2018), the Gumbel-softmall trick is used to develop a framework for interpreting pre-trained models. Their method is focused finding subset of features given particular instance, and therefore is not appropriate for general feature selection.
Some studies tackle embedded feature selection problems by extending LASSO and group LASSO to neural network models. The authors (Li et al. 2016 Scardapane et al. 2017) and (Feng & Simon 2017) have similar goal ours in performing feature selection, but instead rely the & relaxation to the fo. Our approach, which utilizes stochastic gates coupled with the tu norm achieves substantially better empirical performance compared against other & based baselines. This point is demonstrated in the next section.
EXPERIMENTS
Here we provide empirical evaluation of our method in wide range of settings. We begin with simple linear regression evaluating the capabilities of our approach to select a sparse set of relevant features (section 6. We then evaluate our method in nonlinear regression tasks (section 6.2). 1.2..1m Next, we demonstrate the applicability of the method in a highly nonlinear classification task (section 6.3 Finally, we present two potential applications (section 7)
We implement our method using both the STG and HC distributions and compare their applicability for feature selection. The hyperparamters of all the methods are optimized using validation sets via Optuna (Takuya Akiba & Koyama| 2019). In the Appendix, we provide all details the Optuna based parameter tuning procedure. Additional experiments also appear in the Appendix. including high dimensional datasets and an extensive comparisons between the STG and HC (Louizos et al. 2017) distributions.
5

---

Under review as a conference paper at ICLR 2020
6.1 SUPPORT RECOVERY IN LINEAR REGRESSION
In the setting of noisy linear regression, (Wainwright 2009) have analyzed the probability of LASSO correctly identify sparse subset of active variables. The problem is know as support recovery and is formulated as follows; let f E RD be fixed sparse vector, such that 8 E {-0.5, 0.5) (with equal probability) if S, and 3 = 0 otherwise. Suppose the cardinality of the support |S] is known. Given matrix of measurements X E RxN with values drawn independently from N (0, 1), the response y is defined as
y X 3 w,
where the values of the noise w,, 4 1,..., N are drawn independently from (0,00.5)
Here, reproduce this setting evaluate the probability of the proposed approach perfectly recover the support of 3* As suggested by (Wainwright, 2009), we use sparsity that scales with D such that & 0..4 For each number of samples N in the range [10, 500], run 200 simulations and count the portion of correctly recovered supports. We repeat this process for 2 different values of and compare our performance LASSO. For LASSO, the regularization
parameter was set to its optimal value aN = 20 - log( log(k)
Wainwright 2009). For STG and
HC we set ' = CaN, such that C is constant, which is selected using grid search in the range [0.1,10]. As evident from Fig. [ even when restricting to linear functions our method has clear advantage over LASSO. This implies that the to based penalty, although is non convex in nature, allows perfect recovery of the support using less samples. Furthermore, the proposed STG requires less samples than the HC distribution and suffers from smaller variance.
Prob. of sucess
Prob. of sucess
0.9
0.0
0.4
STG LASSO
STG
02
IC
LASSO
0.
200
300
samples N
# of samoles N
(a)
(b)
Figure Probability of successfully recovering the support sparse signal 3 E RDD vs. the number of observations N. Comparison between the proposed method, HC and LASSO under
linear additive noise model y = 3 + w using (a) D = 64 and (b) D = 128. Central lines are the
means while shaded are represent the standard deviation.
6.2 NONLINEAR REGRESSION USING SYNTHETIC AND REAL DATASETS
In this section, we evaluate our method for regression tasks against other embedded feature selection methods: LASSO and Sparse Random Fourier Features (Gregorová et al. 2018). Following the same format as (Gregorová al. 2018). the following functions are used generate synthetic data: (SE1: 100/5) y = sin (ax x² sin (#7 + (0, 0.1). (SE2: 18/5) y = S (/ 0.1). (SE3: 1000/10) y =0(=] + 2e - + NO0 0.01), where each z; is drawn from the standard Gaussian (o 1). For (SE3), the five consecutive coordinates are generated by #s(j-)++ 2; + (O 0.1), where 2j (0, for = 1, ..., 5 and j 1, ..., 200. The numbers next to the experiment code indicate total dimensions/relevant dimensions the feature space. We also evaluate our method using three real datasets: (RCP: 21/-), of computer systems activity, (REL 17/-) of F16 elevators and (RAI 39/-) of F16 ailernos all taken from the LIACC repository ? For each dataset, we generate 30 different replications and randomly split the data into train, validation, and test set (see Appendix for more details). The root mean squared error on the test set averaged over 30 random replicated datasets reported in Table 1. Our method outperforms all alternative methods for most cases. Note that (SE1) is generated using sine function, which is in favor of the random Fourier feature based method (SRFF).
xx)?)

---

Under review as a conference paper at ICLR 2020
Table 1: Regression performance comparison in terms of root mean squared error. The mean and
standard deviation are computed across 30 resamples. N is 1000 for SE1,2,3 and 6000 for RCP. The values for LASSO and SRFF are borrowed from (Gregorová al. 2018). DNN represents a deep neural network without feature selection.
EXP
LASSO
RF
SG-L11NN
SRFF
DNN
HC
STG
SEI
0.29 (0.01)
0.30 (0.01)
0.29 (0.01)
0.27 (0.01)
0.29 (0.01)
0.29 (0.01)
0.29 (0.01)
SE2
2.22 (0.10)
2.34 (0.17)
2.35 (0.18)
1.60 (0.10)
2.05 (0.11)
1.19 (0.31)
0.87 (0.15)
SE3
0.68 (0.002)
0.50(0.01)
0.68 (0.01)
0,48 (0.03)
0.73 (0.02)
0.33 (0.05)
0.14 (0.10)
RCP
9.69 (0.71)
3.52 (0.12)
9.64 (0.65)
(0.180
2.89 (0.25)
2.74 (0.75)
2.44 (0.08)
REL
0.47 (0.01)
0.58 (0.01)
0..00)
0..03)
0.61 (0.01)
0.33 (0.04)
0.27 (0.002)
RAI
0.43 (0.02)
0.48 (0.003)
0.47(0.002)
0.41(0.02)
0.44 (0.01)
0.41 (0.01)
0.39 (0.01)
6.3
NOISY BINARY XOR CLASSIFICATION
In the following evaluation, we consider the problem of learning a binary XOR function for classifica- tion task. The first two coordinates 21, I2 are drawn from a binary "fair" Bernoulli distribution. The response variable is set an XOR of the first coordinates. such that y = æ1 I2. The coordinates Iz = 3, ..... D are nuisance features, also drawn from a binary "fair" Bernoulli distribution. The number of points we generate N 1, 500, of which 70 % are reserved for test and 10% of the re- maining training set was reserved for validation. We compare the proposed method to four embedded feature selection methods (LASSO (Tibshirani] 1996). C-support vectors (SVC) (Chang & Lin 2011), deep feature selection (DFS) (Li et 2016), sparse group regularized NN (SG-LI-NN) (Scardapanc) ctt al. (0177). To provide more benchmarks, we also compare our embedded method against three wrapper methods (Extremely Randomized Trees (Tree) (Rastogi & Shim 2000), Random Forests (RF) (Strobl et al.) 2008)) and Extreme Gradient Boosting (XGBOOST) (Chen & Guestrin 2016).
To evaluate the feature selection performance, calculate the Informative Features Weight Ratio (IFWR). IFWR is defined as the sum of weights Wa over the informative features divided by the sum over all weights. In the case of binary weights the IFWR is in fact recall measure for the relevant features (See the Appendix for more details.)
The experiment is repeated 20 times for different values of D, and the average test classification accuracy and standard deviation are presented in Fig. 2(a) followed by the IFWR in Fig. 2(b) The number of selected features affects the accuracy. Therefore, to treat all the methods in fair manner,
we tune the hyperparameter that controls the sparsity level using Optuna (Takuya Akiba & Koyama 2019) which optimizes the overall accuracy across different Ds. For instance, the wrapper methods (Tree, RF and XGBOOST) has threshold value to retain features. We retrain them using only such features whose weight higher than the threshold. terms of feature ranking (see Fig. [[c] only the tree based methods and the proposed (STG and HC based) provide the optimal median rank (which 1.5) for the two relevant features. Nonetheless, the ranking provided by STG is the most stable comparing to all the alternative methods.
Median rank
Accuracy
features
of featurds
(c)
Figure 2: (a) Classification accuracy (mean and standard deviation) the number of irrelevant noisy dimension (D) for the XOR problem. (b) The Informative Features Weight Ratio (IFWR). central lines the means while shaded are represent the standard deviation. IFWR is the sum of weights Wa over the informative features divided by the sum over all weights. (c) Box plots for the median rank of the two informative features. Black line and dashed red represent the median and mean of each method. Optimal median rank in this experiment is 1.5.
APPLICATION
7.1 PURIFIED POPULATIONS OF PERIPHERAL BLOOD MONOCYTES (PBMCS)
In (Zheng et al. 2017), the authors have subjected more than 90, 000 purified cell populations of peripheral blood monocytes (PBMCs) to Single-cell RNA sequencing (scRNA-seq) analysis. Here
7

---

Under review as a conference paper at ICLR 2020
Table 2: Performance comparison of survival analysis on METABRIC. We run the same experiment 5 times with different train/test split and report the mean and the standard deviation on the test set. In (Katzman et al., 2018). it is reported that DeepSurv outperforms other existing survival analysis methods such as Random Survival Forest (RSF) (Ishwaran et al. 2008) and the original Cox Propotional Hazard Model.
DEEPSURV
RSF
CoX-LASSO
Cox-IIC
CoX-STG
C-INDEX # FEATURES
0.612 (0.009)
0.626 (0.006)
0.580 (0.003)
0.636 (0.007)
0.633 (0.005)
221 (ALL)
221 (ALL)
44 (0)
8 (0.89)
2 (0)
we use this data and focus on classifying two subpopulations of T-cells, namely Naive and regulatory T-cells. We use the proposed method to select subset of genes for which the network discriminates between Naive and regulatory T-cells. We first filter out the genes that are lowly expressed in the cells, which leaves us with D = 2538 genes (features). The total number of cells in these two classes is N 20742, of which we only use 10% of the data for training. We apply the proposed method for different values of X and report the number of selected features and classification accuracy on the test set. Here compare our performance (STG and HC) to RF and LASSO. A least squares polynomial fit plot of the accuracy vs. number of selected features is presented in Fig. 3 The accuracy obtained by a NN without feature selection is 91.09%. which is comparable what achieve with a small fraction of the features. We have also evaluated the performance of the layers (HC-Full), following the procedure in (Louizos al., (017). Empirically we observed that using this type of regularization across all layers provides inferior capabilities in terms of feature selection. Moreover, when the HC-Full converges to a larger subset of features (d > 50), it does not generalize at all and the test accuracy is around 0.5.
applied to all
Accuracy
0.52
0.55
0.57
HC
LASSO
0.6
RF
015
# of features
Figure 3: Classification of T-cells sub-populations. Aceuracy vs. number of features selected by each method. Comparison of the proposed method (STG and HC) to Random Forests (RF) and LASSO. 7.2 Cox PROPORTIONAL HAZARD MODELS FOR SURVIVAL ANALYSIS
We incorporate our STG into DeepSurv Katzman et al. 2018) to examine how this procedure improves survival analysis. The input to this analysis are gene expression profiles (Curtis et al. 2012) along with additional commonly used clinical variables, the goal is to predict survival time. this setting, identifying small subset of predictive covariates is important it may lead to cheaper and more stable medical assays. See the Appendix for additional description of the data and the experimental setup.
We compared our method (Cox-STG and COX-HC) against three other methods: Cox model with { regularization (Cox-LASSO), Random Survival Forest (RSF) (Ishwaran et al. (008) and the original SeepSurv. We evaluate the predictive ability of the learned models based on the concordance index (CI), standard performance metric for model assessment in survival analysis, which measures the agreement between the rankings of the predicted and observed survival times. The performance in terms of the CI reported in Table We see that Cox-HC and Cox-STG outperform the alternative methods, indicating that our method shrinks the feature size while learning predictive model.
8 CONCLUSION
In this paper, we propose a novel embedded feature selection method for neural networks and linear models based on stochastic gates. It has advantage over previous E based regularization methods in terms of achieving a high level of sparsity while learning effective non linear models. We justify our probabilistic feature selection framework from the information theoretic perspective. In experiments, we demonstrate that our method consistently outperforms existing embedded feature selection methods in both synthetic and real datasets.
8

---

Under review as a conference paper at ICLR 2020
REFERENCES
Genevera I Allen. Automatic feature selection via weighted kernels and regularization. Journal of Computational and Graphical Statistics, 22(2):284-299, 2013.
Roberto Battiti. Using mutual information for selecting features supervised neural net learning. IEEE Transactions neural networks, 5(4):537-550, 1994.
Yoshua Bengio, Nicholas Léonard, and Aaron C. Courville. Estimating propagating gradients through stochastic neurons for conditional computation. CoRR, abs/1 308.3432, 2013.
Gavin Brown, Adam Pocock, Ming-Jie Zhao, and Mikel Luján. Conditional likelihood maximisation: unifying framework for information theoretic feature selection. Journal of machine learning research, 13(aan):27-66, 2012.
Sohail Chand. On tuning parameter selection of lasso-type methods-a monte carlo study. In Proceedings of 2012 9th International Bhurban Conference on Applied Sciences & Technology (IBCAST), pp. 120-129, IEEE, 2012.
Girish Chandrashekar and Ferat Sahin. survey on feature selection methods. Computers & Electrical Engineering, 40(1):16-25 2014.
Chih Chih-Chung Chang and Chih-Jen Lin. Libsvm: A library for support vector machines. ACM
transactions on intelligent systems and technology (TIST), 2(3):27, 2011.
Jianbo Chen, Mitchell Stern, Martin J Wainwright, and Michael I Jordan. Kernel feature selection via conditional covariance minimization. In Advances Neural Information Processing Systems, pp. 6946-66 6955, 2017.
Jianbo Chen, Le Song, Martin J Wainwright, and Michael Jordan. Learning explain: An cntormatiiontheoretti perspective on model interpretation. arXiv preprint arXiv:1802.0771 08880008.14,,. 2018.
Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pp. 785-794.
ACM, 2016.
Tom Coen, Wouter Saeys, Herman Ramon, and Josse De Baerdemaeker. Optimizing the tuning param- eters of least squares support vector machines regression for nir spectra. Journal of Chemometrics: A Journal of the Chemometries Society, 20(5):184- 192, 2006.
Thomas M. Cover and Joy A. Thomas, Elements of Information Theory (Wiley Series in Telecom- munications and Signal Processing). Wiley-Interscience, New York, NY, USA, 2006. ISBN 0471241954.
Christina Curtis, Sohrab P Shah, Suet-Feung Chin, Gulisa Turashvili, Oscar M Rueda, Mark Dunning, Doug Speed, Andy G Lynch, Shamith Samarajiwa, Yinyin Yuan, et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature, 486(7403): 346, 2012.
Pablo A Estévez, Michel Tesmer, Claudio A Perez, and Jacek M Zurada. Normalized mutual information feature selection. IEEE Transactions on Neural Networks, 10(2):189-201, 2009.
Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American statistical Association, 96(456): 1348-1360, 2001.
Yingying Fan and Cheng Yong Tang. Tuning parameter selection in high dimensional penalized likelihood. Journal of the Royal Statistical Society: Series (Statistical Methodology). 75(3): 531-552, 2013.
J. Feng and N. Simon. Sparse-Input Neural Networks for High-dimensional Nonparametric Regres- sion and Classification. ArXiv e-prints, November 2017.
Mikhail Figurnv,, Shakir Mohamed, and Andriy Mnih. Implicit reparameterization gradients. In
Advances in Neural Information Processing Systems, pp. 441-452, 2018.

---

Under review as a conference paper at ICLR 2020
Magda Gregorová, Jason Ramapuram, Alexandros Kalousis, and Stéphane Marchand-Mailler Large- scale nonlinear variable selection via kernel random features. arXiv preprint 1844.0:080 2018.
Chris Hans. Bayesian lasso regression. Biometrika, 96(4):835 845, 2009.
Hemant Ishwaran and Udaya B. Kogalur. Random survival forests for r, 2007.
Hemant Ishwaran, Udaya B. Kogalur, Eugene H. Blackstone, and Michael S. Lauer. Random
survival forests. Annals of Applied Statistics, 2(3):811-600, 2008. ISSN 1932-6157, doi:
11214/08-AAAS66
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax, 2017. URL htttps://axxii.....// www.aaiwwi.....aaaaooooooooooooe abs/1611. 01144
Md Monirul Kabir, Md Monirul Islam, and Kazuyuki Murase. A new wrapper feature selection approach using neural network. Neurocomputing, (((1611))3273-3283, 2010.
Jared L. Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval personalized treatment recommender system using a cox proportional hazards deep neural network. BMC Medical Research Methodology, 18, 2018.
Kluger. Deepsurv:
Ron Kohavi and George H John. Wrappers for feature subset selection. Artificial intelligence, 97
(1-2):273-324, 1997a.
Ron Kohavi and George John. Wrappers for feature subset selection. Artificial intelligence, 97 (1-2):273-324, 1997b.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278-22 1998.
Johannes Lederer and Christian Müller. Don't fall for tuning parameters: tuiig--free variable selection in high dimensions with the trex. Twenty-Ninth AAAI Conference on Artificial Intelligence,
2015.
Fan Li, Yiming Yang, and Erie P Xing. From lasso regression feature vector machine. In Advances in Neural Information Processing Systems, pp. 779-786, 2006.
Wei Li, Jianxing Feng, and Tao Jiang. Isolasso: a lasso regression approach to rm- based transcriptome assembly. In International Conference on Research in Computational Molecular Biology, pp. 168-188. Springer, 2011.
Yifeng Li, Chih- Yu Chen, and Wyeth W Wasserman. Deep feature selection: theory and application to identify enhancers and promoters. Journal of Computational Biology, 23(5):322-336, 2016.
Christos Louizos, Max Welling, and Diederik P. Kingma. Learning sparse neural networks through 10 regularization. CoRR, abs/1712.01312, 2017.
Chris Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. CoRR, abs/1611.007// 2016. http:///aa/
Stg/abs/Ine 00712
Andrew Miller, Nick Foti, Alexander D' Amour, and Ryan P Adams. Reducing reparameterization gradient variance. In Advances Neural Information Processing Systems, pp. 3708-3718, 2017.
Fan Min, Qinghua Hu, and William Zhu. Feature selection with test cost constraint. International Journal of Approximate Reasoning, 55(1):167-179, 2014.
Hanchuan Peng, Fuhui Long, and Chris Ding. Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy. IEEE Transactions pattern analysis and machine intelligence, 27(8): 1222-1238. 2005.
Rajeev Rastogi and Kyuseok Shim. Public: A decision tree classifier that integrates building and pruning. Data Mining and Knowledge Discovery, 4(4):315-344, 2000.
10

---

Under review as a conference paper at ICLR 2020
Juha Reunanen. Overfitting in making comparisons between variable selection methods. Journal of Machine Learning Research, 33Mar):1371-1382) 2003.
Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. "why should trust you?": Explaining the predictions of any classifier. In Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '16, pp. 1 1135-1144, 1144, New York, NY, USA, 2016.
ACM. ISBN 978-1-4503--232-22 doi: 10. 10.1455/29997722933978 URL battp://doot acm.
org/10. 000//00.119/299967722
Debaditya Roy, K Sri Rama Murty, and C Krishna Mohan. Feature selection using deep neural networks. In Neural Networks (IJCNN), 2015 International Joint Conference on, pp. 1-6. IEEE, 2015.
Simone Scardapane, Danilo Comminiello, Amir Hussain, and Aurelio Uncini. Group sparse reg- ularization for deep neural networks. Neurocomput, 2C1(C):8 June 2017. ISSN 0925- 2312. doi: 0..1016/j.eucco..c www.::.....1........ 7.02.029. URL httpt::/t/ooooo http://oooorrg////// 101 6/3 neucom 2017.02. 029
Le Song, Alex Smola, Arthur Gretton, Karsten Borgwardt, and Justin Bedo. Supervised feature selection via dependence estimation. In Proceedings of the 24th international conference Machine learning, pp. 823-830. ACM, 2007.
Le Song, Alex Smola, Arthur Gretton, Justin Bedo, and Karsten Borgwardt. Feature selection via dependence maximization. Journal of Machine Learning Research, (3(May):1393-1443) 2012.
Gary Stein, Bing Chen, Annie S Wu, and Kien A Hua. Decision tree classifier for network intrusion detection with ga-based feature selection. In Proceedings of the 43rd annual Southeast regional conference-Vee Volume 2, pp. 136-1411 ACM, 2005.
Carolin Strobl, Anne-Laure Boulesteix, Thomas Kneib, Thomas Augustin, and Achim Zeileis. Conditional variable importance for random forests. BMC bioinformatics, 9(1):307, 2008.
Toshihiko Yanase Takeru Ohta Takuya Akiba, Shotaro Sano and Masanori Koyama. Optuna: A next-generation hyperparameter optimization framework. KDD, 2019.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological)) pp. 267-288, 1996.
Antanas Verikas and Marija Bacauskiene. Feature selection with neural networks. Pattern Recognition Letters, 1333133333331335 2002.
Martin J. Wainwright. Sharp thresholds for high- sigh-dimensional and noisy sparsity recovery using
en-constrained quadratic programming (lasso). IEEE Transactions on Information Theory, 55
(5):2183 2202, May 2009. ISSN 0018-9448,
doi: 10. 211100///..2002200000000
URL http:
//dx. doi org /10 1.....1199tiie. 2009. 2016018
Ronald Williams. Simple statistical foolowng algorithms for connectionist reinforcement learning. Machine Learning. 8:229-256, 1992.
Makoto Yamada, Wittawat Jitkrittum, Lconid Sigal, Eric P Xing, and Masashi Sugiyama. High- dimensional feature selection by fcature-wise kernelized lasso. Neural computation, 26(1):185-207, 2014.
Grace XY Zheng, Jessica M Terry, Phillip Belgrader, Paul Ryvkin, Zachary W Bent, Ryan Wilson, Solongo B Ziraldo, Tobias D Wheeler, Geoff McDermott, Junjie Zhu, et al. Massively parallel digital transcriptional profiling of single cells. Nature communications, 14049, 2017.
Zexuan Zhu, Yew-Soon Ong, and Manoranjan Dash. Wrapper-filter feature selection algorithm using a memetic framework. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics). 37(1):70-76, 2007.
11

---

Under review as a conference paper at ICLR 2020
A APPENDIX
B ALGORITHM
Here present the pseudo code of our method, presented in Algorithm The loss L is typically negative-log likelihood for classification and squared loss for regression. N the sample size, D is the number of features, K is the number of Monte Carlo samples.
Input: X RNDD target variables RV regularization parameter number of epochs M, learning rate
Output: Trained model fg and parameter μ RDD Initialize the model parameter 0. Set μ 0.
2: for i 1, .,.., M do
3: for n = 1,.., N do
4: for d = 1 ..... D do
5: for 1,... K do
6:
Sample dd
(k)
(0,
7:
Compute the gate 22 = max(0, min(1, Ha + cd - -0.5))
8:
end for
9:
Set 2d
K Σ== Rd
end for
11:
Set = = [Z1,,.., =D!²
end for
Compute the loss = N EN ((f @ 2))
13:
Compute the regularization term R > $ mo+0..5
Update 0 = 0 VV and μ == μ Vu R)
16:
end for
Algorithm STG: Feature selection using stochastic gates
We empirically observe that setting the number of Monte Carlo samples K K= 1 and the standard deviation of the Gaussian distribution ay 0.5 for d I,..., suffices for feature selection in our experiments. See Section E in the Appendix for more details about the specific choice of aa. After training. the set of indices for selected features is: {d: min(1, max(0, pa + 0.5)) 0}.
C PROOF OF PROPOSITION
We now provide proof for Proposition 1, showing the equivalence between the stochastic optimiza- tion (Eq 7 and the deterministic one (Eq. 6p. Let S be subset such that S* \ S * 0. That is there exists some element in S* that is not in S. For any such set S we have that T(X gi Y) IX Indeed, if let E S* n Se then we have
T(X Y) S T(XM)):
= =((XYY) -
< <X())),
where the final inequality follows by Assumption 1. Assumption 2 also yields that for any set Š such that S* c: S, we have I(X 3: Y) YY Now, when we consider the Bernoulli optimization problem we have
T(X @ S:Y) s.t. Σ T k and 0 क < 1.
The mutual information can be expanded as
1(X o S;YY= = Σ 11
12

---

Under review as a conference paper at ICLR 2020
where we have used the fact that S is independent of everything else. Recall that in optimization (Eq. 7 the coordinates of are sampled at random. Therefore, the distribution that is being optimized over Prr is product distribution. Our goal is to understand the form of this distribution. To that end, we will consider problem dropping the independence constraint. If we can show that the distribution found by solving this optimization problem with less constraints is still product distribution, then we obtain a solution to the original optimization (Eq. 6
Now, from above we know that the optimal value of the optimization [(X @ S; Y) for any set S* C S. Hence, any unconstrained distribution should place all of its mass on such subsets in order to maximize the mutual information. As a result Eies- p(S 1) k. However, there is an optimization constraint that E Si] k. Therefore, E, [t] solution is to select the distribution so that all of the mass is placed on the subset 5 and no mass elsewhere. As this is also product distribution, this complete the proof of the claim.
for any & S Hence, the optimal
D BRIDGING THE TWO PERSPECTIVES
To motivate the introduction of randomness into the risk, we have looked at the feature selection problem from MI perspective. Based the MI objective, have observed that introducing randomness into the constrained maximization procedure, does not change the objective (Proposition
1). Here we provide relation between the MI objective (Eq. 6 the empirical risk (Eq. D which supports our proposed procedure.
We first note that the MI maximization over the S can be reformulated the minimization of the conditional entropy HY since HYY does not depend on S:
max s: Y) max H(Y) (Y S
min HY))ss).
S
Recall that Xs X S. By Proposition 1, rewrite the deterministic search over the set S by a search over the Bernoulli parameters π:
min H X 5 = min E x,Y - log Pg- (Y 5 π
min min E x log Pg (Y O))
where the expectation is over X, Y ~ Po.: which is the true data distribution, and S Bern(S]m Put our model distribution Po, then we can rewrite the right hand side as:
PY X o 5) Po X X
Ex log Po. 6Y| @ 5) E 3 log Po- (Y X C S)
PY X S)
Ex, log
E xx log Po 5).
P O S)
Since KL(P-- [Y X C PY + $) is non-negative, ELKL (Y @ S) X + also non-negative because is a weighted sum of non-negative terms. Noting that
PY o S)
(yK @ X @ 0)) = E Ex log
Po X o $)
we can conclude that
E Ex Po- (Y| @ 5) 5 Ex - log Pe 3 5).
If we consider the negative log likelihood of the target given the observations (i.e. log Po ( as a loss function L (which encodes the classification or regression function fo), then we see that minimizing the risk approximately maximizes the MI objective in Eq. 6
13

---

Under review as a conference paper at ICLR 2020
E DETAILS OF REGULARIZATION TERM
Here we provide the detail description of the regularization term. For the vector of stochastic gates I. the regularization term is expressed follows:
Ez
12
Piza > 0J = Σ
Plus + oaed -
d=1
dm
Σ
D
(1 Plud + oded S - 1
-
d=1
D
3 Md }}
(1
d=1
+
$
dd
d=1
The derivative of the regularization term with respect to the distribution parameter pa is simply the Gaussian PDF:
Ez 1 0pa
i
a
Md
2)
1
(
2n2
Oya
0d
2
The effect of can be understood by looking at the value of opa 1 In the first iteration during the training, l4a is O. Therefore, during the initial phase of training, it is close V2naz In order remove irrelevant features, this term has to be greater than the derivative of the loss with respect μa because otherwise μ4 is updated in the incorrect direction. To encourage such behavior, we set a 0.5, which around the maximum
Bo
gradient during the initial phase as shown in Fig. 4 Although the point that attains the maximum moves as changes, we empirically observe that setting o = 0.5 performs well our experiments when the regularization parameter is appropriately tuned (see Section
16 2/ at p=0
0.2
0000
Figure 4: The plot of EZ
1olp=0
for [0.001, 2].
V2
F ISSUES IN GRADIENT ESTIMATION OF DISCRETE RANDOM VARIABLES
In Section 3.1 we have introduced Bernoulli random 3., d = 1,.., D variables with corresponding parameters тu into the risk objective (Eq. 4 Taking the expectation over the & norm of $ boils down to the sum the Bernoulli parameters тa. However, the optimization of the resulting objective suffers from high variance due to the discrete nature of S. Here, we attempt to convey this problem by analyzing the risk term in the objective in Eq. 4 Using the Bernoulii paramterization the empirical risk R(, ) is expressed as
14

---

Under review as a conference paper at ICLR 2020
E 2 a n==1
T ynl d=1
T (1 nT2"" "d II
In practice, as the outer sum involves enumerating 2n possibilities of the indicator variables, one can replace the outer sum with Monte Carlo samples from the product of Bernoulli distributions (=) However, Monte Carlo estimate of Ōπa R(O, π suffers from high variance. To see this, consider the following exact gradient of the empirical risk with respect to TT. which
Σ
[2))
E (: .5
[[())Pscaa]
where p(Ziza) = 1eipa TH -1- -.) (1-x;)) by absorbing the model fo(-) and data into L(-). Due the discrete nature of z, see that even the sign of the gradient estimate becomes inaccurate we can only access small number of Monte Carlo samples. While score-function estimator such as REINFORCE (Williams) T992) can used, is known that the reparametrization trick reduces the variance more in practice.
G COMPARISON TO HARD-CONCRETE DISTRIBUTION
To evaluate the strength of the proposed continuous relaxation of the Bernoulli distribution, described in Subsection 3.1, we compare with the Hard-Concrete distribution (Louizos et al.) (207).. another continuous surrogate for Bernoulli distributions, which was originally developed for neural network model compression. The details of the Hard-Concrete distribution is described below.
G.1 HARD-CONCRETE DISTRIBUTION
The authors in (Louizos et al. 2017) introduce a modification of Binary Concrete, whose sampling procedure is follows:
u ~ U(0, 1), L= =Uo((U) -log(1 U)
exp( - 1log a+)))
8 = s(5 - 7)
min(1, m0, 8))
where (T, is interval, with and <> 1. This induces new distribution, whose support is [0, 1] instead of (0, 1). With < 3 < 1, the probability density concentrates its mass near the end points, since values larger than (-r are rounded one, whereas values smaller than rounded to zero.
The CDF of is
Qs(3, Sigmoid( log(1 s))B log a))
(9)
and the CDF of is
5 5-7 (s(5))) Sigmoid (log(
( )- log(1 --1))² 1),8 log a
(10)
C-7 7
where ø = (3, log a, T)
Now, the probability of being the gate 2; being active is - Q=(0lØ) and can be written
Q. Sigmoid(log Blog -) )
(11)
15

---

We define a deterministic gate using žα = g(pa) = max(0, min pa + 0.5)), where pat a define the Under review as a conference paper at ICLR 2020
G.2
THE EFFECT OF THE HEAVY TAIL DISTRIBUTION
The main difference between our proposed distribution (STG) and the Hard-Concrete (Louizos et al. 2017) distribution is that the latter based on the logistic distribution, which has heavier tail than the Gaussian distribution we have employed. As shown in Fig the heavy-tailness results in instability during training. Furthermore, our method converges much faster and more reliably than the feature selection method using the Hard-Concrete distribution on a the two-moons, XOR and MADELON datasets (see Subsection and J.2p. A similar phenomenon also demonstrated in the MNIST experiment (see Subsection
JFW
STG
STG Conc-Concidte
epachs
(a)
(b)
(c)
Figure 5: Comparison between STG and Hard- Concrete the Two-Moon dataset XOR (b) and MADELON The shaded area represents the standard deviation, calculated by running the same experiment 10 times with different random initializations.
H NON-CONVEX DETERMINISTIC REGULARIZATION
Our proposed objective for feature selection (summarized in Eq. consists of three major differences compared to the LASSO formulation.
Non-linearity, which is obtained by learning parameters of a multi-layer network with non linear activation's s.
Feature gates, these are specific form on non convex regularization.
Stochasticity, which is achieved by injecting noise via the reparametrization trick.
In section 6.1 have demonstrated that even in the linear regime the proposed formulation outperforms the LASSO in the task of identifying active features (see Fig. [(a)).
The stochastic property of the gates proposed in this paper is justified using a MI prospective (see Section 4. Nonetheless, this section we evaluate whether the practical performance gain of our method stems from the non-convex regularization from it's combination with the injected noise.
parameter learned for each feature 1,... D. In the linear regression setting non-convex deterministic (NCD) as
N Σ min 0,p N
D Σ d=
[dd + 0.5
10" xr O )2 +
(12)
0.5
where the standard Gaussian CDF. Combined with & this non-convex regularized objective is deterministic and differentiable, and its solution can be approximated via gradient decent. We evaluate the applicability of such alternative objective for the task of feature selection in linear regression. The experimental setting follows the description in Section 6.1 using D = 64.
Fig. 6(b) we compare the values the deteministic gates (DNC) and stochastic gates (STG) using N 60. As depicted from this figure, the STG correctly identifies the active features (first 10) and sparsifies the non-active features. The deterministic version DNC does induce partial sparsity, but fails to remove some of the irrelevant features. Furthemore, it does not push the gates values of active features up to 1. Clearly, ; in the deterministic setting, shrinking the value of μa while compensating with large coefficient in 0 is preferable for the objective in Eq. 12
16

---

Under review as a conference paper at ICLR 2020
To evaluate the probability of support recovery using DNC, we apply additional thresholding to the gates. Specifically, we use a threshold of 0.5 to define if feature is active or not. A comparison with STG, HC and LASSO in this setting appears in Fig. 6(a] The DNC is superior LASSO, but inferior to the HC and STG. Which means that there is true advantage of stochacity in our formulation.
One aspect that is crucially difference between the STG and DNC what we call "second chance". For the deterministic formulation, if feature is zeroed out in an early training phase, the gradient of the feature vanishes for further training and will never be pulled back. For the STG, even the feature has been pulled zero, the gate may still be pushed back up one at a later training phase. This is due to the injected noise which allows to reevaluate the gradient of each gate. In Fig. [6() we demonstrate this "second chance" effect using N 60 samples and presenting the gate's values (throughout training) for an active feature.
To conclude, this form of deterministic gates (Eq. 12 is not suitable for an embedded feature selection method. Using thresholding seems somewhat effective; however such thresholding was not required for the STG in any our experiments. While the shrinkage phenomenon might be avoided using further regularization on 0 this is out of the scope of the current study.
1.0
Prob. of sucess
0.8
HC
DNC
0.0
50
100
150
200
of samples
(a)
DNC STG
Gate value
0.8
0.0
62
0.4
%
Feature index
of epoch
(b)
(c)
Figure 6: Support recovery in the linear regression (see Section 6. the dimension 64. a
Phase transition comparison between the STG, HC, LASSO and the deterministic non convex (DNC). An example of the output gates, STG (z.) and DNC (2.) for N 60. Each marker represents the gate value of the corresponding feature, of which the first 10 features are active. The deterministic
gates are not pushed towards 1. (c) Demonstration of the "second chance" inherit to the STG. In the deterministic setting if feature is eliminated its gradient vanishes for the rest of the training. The STG reevaluates this features and "recovers" after 4000 epochs.
17

---

Under review as a conference paper at ICLR 2020
I TUNING THE REGULARIZATION PARAMETER À
There have been several studies (Lederer & Müller 2015 Coen et al. 2006) on the problem of tuning regularization parameter regression and classihccation. Perhaps the most related regularized
objective is the LASSO formulation. Studies such as (Chand, 2012 Fan & Tang 2013) develop methods for tuning the regularization parameter LASSO. In Eq. 5 as the LASSO formulation, increasing the value of X effectively forces sparser solution with less features. For simplicity, we focus on one popular way for tuning the LASSO. The idea is to minimize the generalization error using cross validation procedure over range of À values.
In this section, we demonstrate that this procedure is also applicable for the proposed approach. We apply the STG the RAI dataset and evaluate the MSE using 5-folds cross validation procedure. We compute the MSE on the validation data using 100 values logarithmically scaled [10-3, 10+]. As demonstrated in Fig. using cross validation procedure for tuning provides clear optimal value for this example.
0.70
0.65
0.60
0.55
MSE
0.50
0.45
0.40
0.35
0.30
10
10
10
10"
10A
X
Demonstrating a cross validation procedure on the RAI dataset (see section ! for details). We perform 5-folds cross validation and evaluate the MSE on the test test. The optimal X 0.135 which seems stable in this example based on the low standard deviation.
Figure 7:
J ADDITIONAL EXPERIMENTS
J.1 MADELON DATASET
The MADELON dataset, first suggested for the NIPS 2003 feature selection problem, is a multivariate highly nonlinear binary classification problem. The MADELON dataset is generated using 32 groups of data points placed on a 5 dimensional hyper-cube and randomly labeling them by one of the two class labels. The first 5 informative features are then used create 15 additional coordinates that are formed based on random linear transformation of the first 5. A Gaussian noise (0,, 1) is added to each feature. Next, additional 480 nuisance coordinates are added in the same manner. These features have no effect on the class label. Finally, 1% of the labels are flipped. To be consistent with the XOR and two moons experiments, we use 1, 500 points from this dataset, and evaluate our proposed method in terms of its predictive power and ability to detect the informative features. We vary the regularizaton parameter À in the range [0.01, 10 and evaluate the classification accuracy using folds cross validation. In this example, focus on Random Forest as was the strongest competitor to our method all of our experiments, while LASSO is evaluated because it is widely used embedded feature selection method.
restrict comparison to Random Forest and LASSO. We
18

---

Under review as a conference paper at ICLR 2020
8(a) method achieves the highest accuracy while using less features. Moreover, as depicted from this figure, peak performence occurs when selecting 5 features, thus, our method provides clear indication to the true number of informative features. Both LASSO and RF on the other hand, do not provide clear indication of the true number of relevant features. Fig. 8(b) we evaluate the effect of À on the number and quality of selected features. As shown in this plot, there is wide range of A's such that our method only selects relevant features. Finally, as evident from the plato the right hand side of the red plot, there range A'S such that exactly 5 features
As evident from Fig.
are selected.
12
selected features
0.95
20
Accuracy
0.40
0-
0.5
20.4.00
IFWR
0.6
0.75
SIG
0.70
RF
IM LASSO
065
100
200
10
100
# of selected features (a)
(b)
5.00
4.75
Median rank
4.50
4.25
3.75
3.50
3.25
LASSO
RF
STG
(c)
Figure 8: (a) Classification accuracy on the MADELON data sets. We evaluate performance using 5-fold cross validation for different number of selected features. In this dataset, only the first 20 coordinates informative. that regime the proposed method outperforms RF and LASSO. (b) An empirical evaluation of the effect the regularization parameter X. The IFWR and the number of
selected features are presented on both sides of the y-axis of this plot. For both plots, the mean presented as solid/dashed line, while the standard deviation is marked mean. [c) Box plots for the median rank the 5 original informative features. Black line and dashed red represent the median and mean of each method. Optimal median rank in this experiment is 3.
shaded color around the
J.2 Two MOONS CLASSIFICATION WITH NUISANCE FEATURES
In this experiment, we construct a dataset based on "two moons" shape classes, concatenated with noisy features. The first two coordinates Ii, 42 are generated by adding Gaussian noise with zero mean and the variance of o² 0.1 onto two nested half circles, presented in Fig. 9()) Nuisance features T,, = 3, .... D, are drawn from a Gaussian distribution with zero mean and variance of o = 1. We reserve the 70% as test set, and use 10% of the remaining training set as a validation set. We follow the same hyperparameter tuning procedure as in the XOR experiment. The classification accuracy is in Fig. [b)) Based the classification accuracies, it is evident that for small number of nuisance dimensions all methods correctly identify the most relevant features. The proposed method (STG) and Random Forest (RF) are the only methods that achieve near perfect classification accuracy
19

---

Under review as a conference paper at ICLR 2020
for wide range of nuisance dimensions. The other NN based method (DFS) seem to converge to sub-optimal solutions. We note that the median rank for all the methods 1.5.
10
Accuracy
0.81
0.0
0.7
LASSO
0.0 10 20 30
60
70
100
X1
# of features
(a)
(b)
Figure 9: Realization from the "Two moons" shaped binary classification class. X1 and X₂ are the relevant features, X,, i ..., D are noisy features drawn from a Gaussian with zero mean and variance of 1. b Classification accuracy (mean and standard deviation based on 20 runs) vs. the number irrelevant noisy dimension.
J.3 SPARSE HANDWRITTEN DIGITS CLASSIFICATION
In the following toy example, we attempt to distinguish between images of handwritten digits 3's and 8's using samples from MNIST (LeCun et al. 1998). The orientation and location of the digits is more or less the same throughout this dataset, therefore for these two classes (3's and 8's), we expect that of the left side features (pixels) would be sufficient for the separation. The experiment is performed as followed. We reserve 90% of the data the test set, and train on the remaining 10%. We then apply STG and evaluate the classification accuracy and the number of selected features. We use the architecture [200, 50, 10] with tanh activations. The experiment was repeated 10 times, the extracted features and accuracies were consistent over 20 trials. We noticed relatively small number of selected features, which are positioned southwest and close the center of the images, achieve very high classification accuracy. An example of randomly selected samples overlaid with also evaluate the effect of À on the the number of selected features and accuracy of the method. We apply our method distribution
the weights of the selected features is presented in Fig. T0(a) Furthermore,
randomly sampled training set of size V=1500 and vary À in the range of (0.001, 0.01]. ]. In Fig. [0(b) we present the accuracy and sparsity level vs. the À parameter. This experiment demonstrates the improved performance of the proposed distribution compared to the Hard-Concres neural model compression. Not only that the overall accuracy is superior, but also it seems that the transition as function of À is smoother, which suggests that the method is less sensitive to the choice of X.
(STG) and its variant using the
(HC (Louizos et al. 2017)). which was designed for
J.4 GISETTE DATASET
The Gisette is a handwritten dataset also appears in the NIPS 2003 feature selection challenge. The data consists 13, 500 handwritten digits of and '9'. The original digits have x 28 pixels, which were modified for the feature selection challenge. Specifically, random subset of the features are embedded in a 2500 dimensional space and additional 2500 irrelevant probes are added. The train, validation and test dimensions are 6000/ 0000/1100/ /5000 respectively. We apply the proposed STG based feature selection method and compare its performance to Random Forests (RF) and Extremely Randomized Trees (Tree). As evident from Fig. [I] obtain high accuracy for dramatic reduction in the feature size. LASSO not presented in this experiment as its performance is dramatically inferior to the alternatives (accuracy < 0.6).
20

---

Under review as a conference paper at ICLR 2020
3
3
3
Accuracy-HC Sparsity-HC Accuracy-STG
8
3
1
3
Sparsity STG
3
8
3
(a)
(b)
Figure 10: (a) Nine samples from MNIST (white) overlaid with the subset of 13 features (black) selected by STG. Based on only these features, the binary classification accuracy reaches 92.2%. For these nine randomly selected samples, all the 8's have values within the support of the selected features, whereas for the there no intersection. The comparison of accuracy and sparsity level performance for À in the range of [10-3 10 between using our proposed method (STG) and its variant using the Hard-Concrete (HC) distribution.
0.8
Accuracy
0.00
STG Tree
# of features
Figure 11: Classification of the binary noisy digits from the Gisette dataset. The total number of feature is 5000 of which 2500 irrelevant probes. Here compare the performance of the proposed approach to RF and Tree based classifier. The lines represent a least squares polynomial fit plot of the accuracy number of selected features.
J.5 REUTERS CORPUS VOLUME
The Reuters Corpus Volume (RCV1) consists of 800, 000 newswire stories manually labeled by 103 categories. This is multilable regime, i.e. each story is assigned to multiple categories. Here focus on a binary subset of this corpus, with 23, 203 stories. The total number of feature is 47, 236 and the train, validation and test portions are 10%, 8.5%, and 81.5% respectively. We evaluate the performance of our method using 5-fold cross validation. A comparison the LASSO and RF appears in Fig. 12 This example demonstrates that our method is also effective in an extremely high dimensional regime of non linear function estimation.
J.6 DETAILS OF COX PROPORTIONAL HAZARD MODEL
Survival times are assumed to follow a distribution, which is characterized by the survival function S()) P(T 1). A hazard function, which measures the instantaneous rate of death, is defined by
h(t) = lim P(tTT (HAT))))))))
p(t) We can relate the two functions in the following way: 5(t)
At
S(t) e k(t)dt
Proportional hazard models assume multiplicative effect of the covariates æ on the hazard function such that
(t)e 0"
h(((=)))
21

---

Under review as a conference paper at ICLR 2020
Accuracy
0:00
0.70
STG LASSO
# of features
Figure 12: Classification of the version of the RCVI textual dataset. The total number of feature is 47, 236. Here compare the performance of the proposed approach to RF and LASSO. The lines represent a least squares polynomial fit plot of the accuracy vs. number of selected features.
where (t) is a baseline hazard function, which is often the exponential or Weibull distribution, and O is the parameter of interests.
One of the difficulties estimating survival analysis is that a large portion of the available data
is censored. However, in order obtain estimates, Cox observed that it is sufficient to maximize the
0 x
partial-likelihod, which is defined as follows: [(()) = Il uncensored
Σrye, x
In (Katzman et al. 2018), the authors propose DeepSurv, which uses a deep neural network model to replace the linear relations between the covariate and 6, demonstrating improvements of survival time prediction over existing models such as CPH and the random survival forest (Ishwaran & Kogalur 2007). (Ishwaran al. 2008).
The Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) dataset consists of gene expression data and clinical features for 980 patients, and 57.72% have an observed death due to breast cancer with a median survival time of 116 months.
The METABRIC dataset involves 24,368 features (genes). Most genes irrelevant for outcome prediction. To demonstrate the advantage of STG in the context of survival analysis, we selected 16 well-known genes relevant for survival (out of the 24,368 genes) that correspond the Oncotype DX test, a gene panel used for treatment decision making. We also include five additional clinical features (hormone treatment indicator, radiotherapy indicator, chemotherapy indicator, ER-positive indicator, and age at diagnosis). We then added 200 additional irrelevant gene variables that we selected randomly from the remaining list genes.
After we omit the null values, the number of samples is 1969, We reserve the 20% for test, and use 20% of the remaining training set as validation. (That is, Train: 1260, Valid: 315, Test: 394 samples.)
Experimental Detail For DeepSurv, we manually select the architecture using the validation set so that we obtain a similar performance reported in (Katzman et al., 2018). The learning rate decay is set to 1. The learning rate and the regularization parameter are optimized via Optuna using the validation set, where the search range set LR [le - 3, 1] and x: [1e- 3, 1]. The hyperparameters the following: architecture :[60, 20, 3], activation: selu (as suggested by
used in the experiment (Katzman et al. 2018)), learning rate 0.152, learning rate decay: 1.0, a : 0.5, 0.023, training epoch: 2000. Note that to see the effect of feature selection, used the hyperparameters optimized for DeepSurv to test our method (STG-DeepSurv).
J.7 ADDITIONAL EXPERIMENTAL DETAILS
Here we provide a full description of the procedures we have performed the experimental parts of the paper.
For synthetic datasets are first split into train, validation and test. Validation is always 10% of the train, while the exact ratios between train and test is detailed for each experiment separately. All the neural network weights are initialized by drawing from (0, 0.1) and bias terms are set to 0. All the batch sizes equal to the number of training samples. In table we detail the search range
22

---

Under review as a conference paper at ICLR 2020
of hyperparamters well as the exact values used in experiments. We set n-trials 1000 for Optuna, which is a bu-rne based hyperparameter optimization software equipped with efficient search and pruning strategies. We use SGD for all the experiments, except for the Cox model where we use Adam. All the experiments are conducted using Intel(R) Xeon(R) CPU E5-2620 v3 @2,4Ghh x2 (12 cores total).
Table 3: List of the search range for the hyperparameters used in our expirements for XOR and Two-Moon
Param
Search range
# dense layers
[1.3]
# hidden units activation LR n-epoch (DFS, SG-LI-NN)
[10, 500]
[tanh, relu, sigmoid]
[le-4, le-1]
[50, 20000]
(SG-LI-NN)
[le-3, 1]
À (SG-L1-NN)
[1e-7, 1]
À (STG, DFS)
[1e-3, 1]
A (LASSO)
[0.01, 1]
n-est (RF, XGBoost, Tree)
[5,100]
n- n-boost-round (XGBoost)
[1,100]
Thresh (RF, XGBoost, Tree)
[0.01,0.5]
max-depth (XGBoost) (SVC)
[0.01,0.5]
[le-7, 1]
For the Phase Transition experiment, we use 0.1 as a learning rate. For all the experiments when we use Tree and RF, we use the default value for max-depth, so that nodes are expanded based on the purity. For the XOR problem, the exact architectures used for the NN based methods are: (STG/HC): [476, 490, 14] with Tanh, (DFS): [100, 10) with Tanh, (SG-LI-NN): [100, 10, 5] with Tanh. For the two moons we use (STG): [490, 406, 18] with Tanh, (DFS): [158, 27, 224] with Tanh, (SG-LI-NN): [88, 28, 27] with Tanh. For the XOR problem, we attempted to use Optúna to optimize parameters of DFS and SG-LI-NN, but we ended up using the architecture suggested by the authors (Li et al.] 2016 Scardapane et al., 2017) as they outperform the values suggested by Optuna. The number of epochs used for the XOR problem is 20K, 14K, 800 for STG/HC, DFS and SG-LI-NN respectively. Regularization parameters are 0.17,33 3.3e 5 and 3e 5 respectively. The number of epochs used for the two-moons problem is 20K, 1570, 708 for STG/HC, DFS and SG-LI-NN respectively. Regularization parameters are 0.48, 9e - and le - 3 respectively. We note that the regularization parameters and learning procedure is different in nature, as we use an to type penalty. For the PBMC experiment, the architecture was hand-tuned using the validation set and set [200, 100, 50, 10) with Tanh. Learning rate was 0.2 and the number of epochs 4000. The hyperparameter X varies in the range [0.001, 0.11] to achieve different levels of sparsity. For MADELON,
the architecture optimized for the binary XOR classification. The number of epochs used is 20K, the learning rate is 0.06 and À varies in the range [0.01, 10). For regression, the learning rate and the are optimized via Optuna using the search range LR [le 4, 1] and A: [Le- - 4, 10] based on validation. The parameters used are the following: (SE1) architecture [600, 200, 100, 50, 1] with ReLu activations, num epochs: 5, : 5, LR 0.0001 (SE2) architecture [600, 300, 150, 60, 20) with ReLu activations, num epochs=2000, À 5, LR 0.001 (SE3) architecture [600, 300, 150, 60, 20] with ReLu activations, num epochs : 1000, À 1, LR 0.005. (RCP) architecture [1000, 300, 150, 60, 20] with ReLu activation. num-epochs: 2000, X: 5.00, LR: 0.001. (REL) architecture [26,91.63] with ReLu activation, num-epochs: 1600, A: 0.031, LR 0.007. (RAI) architecture [10.177] with ReLu activation, num-epochs: 1800, À 0.019, LR: 0.002. Architectures for SE1-SE3 and RCP where tuned manually.
The ratio of train/test/valid split is 1:1:1 for synthetic data. For the real data (RCP and REL), the train size is 6000, the test and valid size is 1000 samples. For RAI the train size is 5000, the test and valid size is 1000 samples.
In order to define the IFWR, for STG, the th feature weight is set to max(0, min(1, pa + 0.5)). For other neural net based methods, it is given by E, Waj. where W is the weight matrix of the first.
23

---

Under review as a conference paper at ICLR 2020
layer. For other methods just used the feature relevance returned by the trained model. Finally, the LASSO's IFWR in the XOR experiment was omitted from the manuscript as suffered from high variance.
Regarding the comparison performed in the two-moons and XOR problem, we believe that adding IFWR along with classification accuracy versus number of feature selected provides complementary perspective in demonstrating the efficacy of feature selection techniques. We emphasize that our goal just rank features but select features by assigning the weight of 0 irrelevant features while simultaneously obtaining good predictive accuracy.
24


Reviews:
[{"review_id": "ryeWRH4TFr", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The author rebuttal sufficiently addresses my concerns, so I am upgrading my score.\n\n***\n\nThe paper considers the problem of embedded feature selection for supervised learning with nonlinear functions. A feature subset is evaluated via the loss function in a \"soft\" manner: a fraction of an individual feature can be \"selected\". Sparsity in the feature selection is enforced via a relaxation of l0 regularization. The resulting objective function is differentiable both in the feature selection and learned function making (simultaneous) gradient-based optimization possible. A variety of experiments in several supervised learning tasks demonstrates that the proposed method has superior performance to other embedded and wrapper methods.\n\nMy decision is to reject, but I'm on the fence regarding this paper. I'm not clearly seeing the motivation for an embedded feature selection method for neural network models: for the datasets considered in the paper, it would seem that training a nonlinear model that used all the features would result in performance at least as good as training the nonlinear model with a prepended STG layer. Perhaps there is evidence that filtering features, e.g., irrelevant features, results in higher accuracy and that the prepended STG layer achieves this accuracy, but that evidence is missing from the paper. Also, there could be downstream computational savings, e.g., at prediction time, if the dimension was very large, but this is not the setting tested in the experiments. I suppose interpretability could be considered motivation, but, even so, isn't there at least one simpler, deterministic approach (described below) that also \"solves\" the problem? Finally, it isn't clear how the method scales with increasing sample size and dimension as all the datasets tested are relatively small in these respects.\n\n***\n\nQuestions and suggestions related to decision:\n\n* The performance values using all features should be included in the experimental results so that the value added by STG can be assessed.\n\n* Why not use the simpler deterministic and differentiable relaxation z = \\sigma(\\mu), where \\sigma() is a \"squashing\" function from the real numbers to [0,1] applied element-by-element to the vector \\mu? What specifically is/are the advantage(s) that the randomness in the definition of z at the bottom of pg. 3 provide over this deterministic alternative?\n\n* Though well-described and methodologically rigorous, the experimental comparison is none-the-less a little disappointing: one dataset for classification and half the datasets for regression are synthetic and low-dimensional. The remaining regression datasets are real but also low-dimensional. The survival analysis dataset is also low-dimensional (as described in the supplementary material). This leaves one real classification dataset which was on the order of 20,000 examples and 2500 features. Why were larger sample-size and dimensionality datasets not tested? These should be readily available. For example, the gisette dataset from the NIPS 2003 feature selection challenge has 5000 features. See \"MISSION: Ultra Large-Scale Feature Selection using Count-Sketches\" by Aghazadeh & Spring et al. (2018) for other high-dimensional datasets. Even a single run for each large dataset would have provided some evidence of scalability.\n\n***\n\nOther minor comments not related to decision:\n\n* \"Concrete Autoencoders for Differentiable Feature Selection and Reconstruction\" by Abid et al. (2019) targets unsupervised feature selection but has enough similarities in the approach that it should be considered related work.\n\n* [Typo?] The unnumbered equation after (5) should not have a sum over d in the second term. Perhaps a sum over k was intended? Also, in this equation, the gradient of the loss wrt/ z samples, average of gradients over z samples times..., does not seem to match what the gradient would be given the algorithmic description in the supplementary material, a gradient of the (sample) average z times...\n\n* The abstract states the paper is proposing a method for high-dimensional feature selection, but all of the experiments have datasets with max. dimensionality 2538.\n\n* Some discussion of how the regularization parameter can be selected by a user of the proposed method would be good to include.", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "rye3z0ytcr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper is concerned with embedding a supervised feature selection within a classification setting. \nThe originality is to use an L_0 regularization (counting the number of retained features), besides the classification loss; the authors leverage the ability to include boolean variables in a neural network and to optimize their value using gradient descent through the reparameterization trick.\n\nI am mildly convinced by the paper:\n* Out of the four contributions listed p. 2, STG is the most convincing one; still, the description thereof is not cristal clear: the reparametrization trick is not due to the authors. The discussion (section 5) needs be more detailed, adding the HC details (presently in appendix); could you comment upon the difference between the proposed STG and the Gumbel-Softmax due to Jang et al, cited ?\n* Likewise the authors delve into details regarding the early state of the art, while omitting some key points. For instance, p. 3, the fact that many authors replaced an L_0 penalization with an L_1 one is rooted on the fact that, provided that the optimal L_0 solution is sparse enough, the L_0 and L_1 problems have same solutions. This section can be summarized;\n* the sought sparsity is assumed to be known, which is bold; \n* Assumption 2 is debatable; one would like to find at most the Markov blanket of the label variable. See Markov Blanket Feature Selection for Support Vector Machines, AAAI 08.\n* There are digressions in the paper which make it harder to follow the argumentation (section 6.1); section 6.2 is not at the state of the art; in Guyon et al's Feature Selection Challenge (2003), the Arcene artificial problem involves a XOR with 5 key features, and 15 additional features are functions of the key features.\n\nSuggestion, you might compare with the L_0 inspired regularization setting used for unsupervised feature selection in Agnostic Feature Selection, Doquet et al, 2019.\n\nDetails: check the citation style: use \\citep instead of \\cite.", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}, {"review_id": "H1e8061p9H", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors propose a feature selection method for high-dimensional datasets that attempts to fit a model while selecting relevant features. \nThe strategy they follow is below:\n\n1. They formulate feature selection as an optimization problem by augmenting standard empirical risk minimization with zero-one variables associated with each feature representing the absence-presence, and adding a penalty proportional to the number of included features. They relax the discrete variables using a continuous relaxation and provide a simple unbiased estimator for the gradient of the relaxation. After training the relaxation is rounded to a zero-one solution by a simple scheme. \n2. They provide an information theoretic motivation for their formulation of feature selection\n3. They exhibit the performance of their method on a number of synthetic and real data scenarios: (i) linear models with a true underlying sparse parameter, (ii) binary classification with a small number of true determining features, (iii) regression performance post-feature selection with synthetic non-linear models (with a few determining features) and two real datasets. They also use the method for a classification problem with RNA-seq data on T-cells and a survival analysis based on a breast-cancer dataset called METABRIC. \n\nDespite my recommendation, there are a number of things that I like about the paper that I list below, along with directions where I believe the article can be improved. \n1. At a certain abstraction, the main idea of the paper is to do feature selection at the same time as model fitting (as the LASSO for e.g. does) while ignoring constraints of convexity raised in the optimization, and simply using stochastic gradient with a reasonable unbiased estimate of the gradient. This is a reasonable idea, particularly if under some reasonable assumptions, the non-convex formulation that is obtained is expected to be computationally 'benign'. \n2. In a number of the experiments, and particularly 6.1 (sparse linear model) 6.2 (noisy XOR classification) I suspect the non-convex formulation is what is providing a lot of the improvement. This has been observed empirically in a number of other settings, for e.g. in matrix completion/factorization problems. Verifying this hypothesis in a simple, synthetic (and therefore controlled) dataset would be a good contribution for a future version.   \n3. The authors have done a fairly good job of validating the method in a number of different settings, even if some of the presentation of their results can possibly be somewhat improved. For e.g. the median rank is better shown with box plots (as in the Chen et al 2018 paper cited by the authors).\n4. There are a number of relaxations of discrete variables used in optimization and theoretical computer science literature. For instance, the approach of the authors is reminiscent to 'mean field' methods, or standard linear programming relaxation of combinatorial optimization problems (i.e. the first level of the Sherali-Adams LP hierarchy). On the other hand, naive versions of this are not likely to work well on (say) sparse linear regression. The current methods do which suggests that the continuous relaxation is useful. \n\nAt an expository level, I also think the paper could do with quite a bit of improvement:\n1. The introduction is sparse and hurried, and without providing sufficient motivation and intuition for the contributions of the article. \n2. In 6.4, 6.5, the introduction about RNA-seq or Cox models can be removed and relevant work cited instead. \n3.  Organizing the experiments as real data, and synthetic data might be semantically better, though that would necessitate splitting Table 1. I am also unclear on why the authors show  performance in Tables 1, 2 independent of the number of features selected, while for the experiment on RNA-seq data the full accuracy/#features tradeoff is given. The sparse explanation about using the Optuna paper is certainly not enough. \n\nMinor comments not related to decision:\n1. The value for \\alpha_N in synthetic sparse linear model experiment of 6.1 likely has an extraneous \\sqrt \\log k \n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors propose a method for feature selection in non linear models by using an appropriate continuous relaxation of binary feature selection variables. The reviewers found that the paper contains several interesting methodological contributions. However, they thought that the foundations of the methodology make very strong assumptions. Moreover the experimental evaluation is lacking comparison with other methods for non linear feature selection such as that of Doquet et al and Chang et al.
Decision: Reject