Title: INVASE: INSTANCE-WISE VARIABLE SELECTION US-ING NEURAL NETWORKS

Abstract: The advent of big data brings with it data with more and more dimensions and thus a growing need to be able to efficiently select which features to use for a variety of problems. While global feature selection has been a well-studied problem for quite some time, only recently has the paradigm of instance-wise feature selection been developed. In this paper, we propose a new instance-wise feature selection method, which we term INVASE. INVASE consists of 3 neural networks, a selector network, a predictor network and a baseline network which are used to train the selector network using the actor-critic methodology. Using this methodology, INVASE is capable of flexibly discovering feature subsets of a different size for each instance, which is a key limitation of existing state-of-the-art methods. We demonstrate through a mixture of synthetic and real data experiments that IN-VASE significantly outperforms state-of-the-art benchmarks.
Let X = X 1 × ... × X d be a d-dimensional feature space and Y = {1, ..., c} be a discrete label space 1 . Let X = (X 1 , ..., X d ) ∈ X and Y ∈ Y be random variables with joint density (or mass) p and marginal densities (or masses) p X and p Y respectively. We will refer to s ∈ {0, 1} d as the

1. INTRODUCTION:

High-dimensional data is becoming more readily available, and it brings with it a growing need to be able to efficiently select which features to use for a variety of problems. When doing predictions, it is well known that using too many variables with too few samples can lead to overfitting, which can significantly hinder the performance of predictive models. In the realm of interpretability, the large dimensionality of the data is often too much information to present to a human who may be using the machine learning model as a support system. Understanding which features are most relevant to an outcome or to a model output is an important first step in improving predictions and interpretability and many works exist that tackle feature selection on a global level. However, in the heterogeneous data we typically encounter, the prediction made by a model (and indeed the true label) may rely on a different subset of the features for different subgroups within the data [14] . In this paper we propose a novel instance-wise feature selection method, INVASE (INstance-wise VAriable SElection), which attempts to learn which subset of the features is relevant for each sample, allowing us to display the minimal information required to explain each prediction and also to reduce overfitting of predictive models.

Discovering a global subset of relevant features for a particular task is a well-studied problem and there are several existing methods for solving it such as Sequential Correlation Feature Selection [11] , Mutual Information Feature Selection [21] , Knockoff models [3] , and more [10; 16] . However, global feature selection suffers from a key limitation -the features discovered by global feature selection are the same for all samples. In many cases, in particular when populations are highly heterogeneous, the relevant features may differ across samples [33; 32] . For instance, different patient subgroups have different relevant features for predicting heart failure [14] . Instance-wise feature selection methods such as [4; 27] instead try to discover the features that are relevant for each sample. When the goal is to provide an interpretable explanation of the predictions made, a key challenge is in ensuring that we do not over-explain by providing too much information (i.e. choosing too many features). Naturally, by performing feature selection on an individualized level we are able to select features that are more relevant to each sample, rather than having to choose the top k features globally, which may not explain the predictions for some samples very well, but simply perform well on average across all samples.

In this paper, we propose a novel instance-wise feature selection method which we term INVASE. We draw influence from actor-critic models [22] to solve the problem of backpropagating through subset sampling. Our model consists of 3 neural networks: a selector network, a predictor network and a baseline network. During training, each of these are trained iteratively, with the selector network being trained to minimize a Kullback-Leibler (KL) divergence between the full conditional distribution and the selected-features-only conditional distribution of the outcome. Our model is capable of discovering a different number of relevant variables for each sample which is a key limitation in existing instance-wise approaches (such as [4] ). We show significant improvements over the state-of-the-art in both synthetic data and real-word data in terms of true positive rates, false discovery rates, and show better predictive performance with respect to several prediction metrics. Our model can also be easily extended to handle both continuous and discrete outputs and time-series inputs (see the Appendix for details).

1.1. RELATED WORKS:

There are many existing works on global variable selection (see [10] for a good summary paper). [21] and [11] use max-dependency min-redundancy criteria [17] with mutual information and Pearson correlation, respectively. [3] uses multiple hypothesis testing for global variable selection. As noted above, these global selection methods are not capable of learning sample-specific relevance.

Instance-wise variable selection is also closely related to model interpretation methods. Some previous works are based on backpropagation from the output of the predictive model to the input variables [29] . DeepLIFT [27] decomposes the output of the neural network on a reference input to compute the contribution of each input variable. However, both methods need white-box access to the pre-trained predictive models to compute the gradient and decomposition. [2] approximates the predictive models using a Parzen window approximator when there is only black-box access to the predictive models. Some other works are based on input perturbation such as [1] , [15] , [30] and [5] . [18] uses Shapley values to compute the variable importance, and [24] uses locally linear models to explain the linear dependency for each sample. [19] tries to interpret tree ensemble models using Shapley values but cannot generalize to other predictive models such as neural networks.

Our work is most closely related to L2X (Learning to Explain) [4] . However, there are 3 key differences between our work and theirs. In L2X, they try to maximize a lower bound of the mutual information between the target Y and the selected input variables X S . In contrast, we try to minimize the KL divergence between the conditional distributions Y |X and Y |X S . In order to be able to backpropagate through subset sampling, L2X use the Gumbel-softmax trick [13] to approximately discretize the continuous outputs of the neural network. In our work, we use methods from actorcritic models [22] to bypass backpropagation through the sampling and instead use the predictor network to provide a reward to the selector network. Finally, due to the Gumbel-softmax used in L2X, the number of variables to be detected must be fixed in advance and is necessarily the same for every sample. The actor-critic methodology used in our model has no such limitations and so we are able to flexibly select a different number of relevant variables for each sample and instead induce sparsity using an l 0 penalty term. In fact, using the actor-critic methodology allows us to directly use the l 0 penalty term (which is not differentiable and therefore not practical to use in general). A summary table highlighting the key features of all of the related works can be found in the Appendix. selection vector, where s i = 1 will indicate that variable i is selected, and s i = 0 will indicate that variable i is not selected. Let * be any point not in any of the spaces X 1 , ..., X d and define

X * i = X i ∪ { * } and X * = X * 1 × ... × X * d .

Given x ∈ X we will write x (s) to denote the suppressed feature vector defined by

x (s) i = x i if s i = 1

* if s i = 0 so that * represents that a feature is not selected.

In the global feature selection literature, the goal is to find the smallest s (i.e. the one with fewest 1s) such that E(Y |X (s) ) = E(Y |X), or equivalently such that the conditional distribution of Y given X (s) is the same as Y given all of X. Note that this definition is given fully in terms of random variables, rather than realizations of those random variables.

In contrast, our problem necessarily needs to be defined in terms of realizations since we are aiming to select features for a given realization. We will write x to denote realizations of the random variable X. Then we formalize our problem as one of finding a selector function, S : X → {0, 1} d such that for almost every x ∈ X (w.r.t. p X ) we have

EQUATION

where

d.

= denotes equality in distribution and S(x) is minimal (i.e. fewest 1s) such that (1) holds.

We suppose that we have a dataset D = {(x j , y j )} n j=1 consisting of n i.i.d. realizations of the pair (X, Y ). 2 Note that Y can be viewed as having either come from a dataset, in which case the problem is of selecting predictive features, or as having come from a predictive model, in which case the problem is of explaining the model's predictions.

2.1. OPTIMIZATION PROBLEM:

In order to learn a suitable selector function, we transform the constraint (1) into a soft constraint using the Kullback-Leibler (KL) divergence which, for random variables W and V with densities p W and p V is defined as

KL(W ||V ) = E log p W (W ) p V (W ) .

We define the following loss for our selector function S

EQUATION

where || • || simply denotes the number of non-zero entries of a vector (or equivalently in this case, the number of 1s) and λ is a hyper-parameter that trades off between the constraint in (1) and the number of selected features. The KL divergence in (2) can be rewritten as

KL(Y |X = x||Y |X (S(x)) = x (S(x)) ) = E y∼Y |X=x log p Y (y|x) p Y (y|x (S(x)) ) = E y∼Y |X=x log(p Y (y|x)) -log(p Y (y|x (S(x)) )) = Y p Y (y|x) log(p Y (y|x)) -log(p Y (y|x (S(x)) )) dy

where p Y (•|•) denotes the appropriate conditional densities of Y . We will write

EQUATION

so that our final loss can be written as

EQUATION

where || • || denotes the l 0 (pseudo-)norm.

3. PROPOSED MODEL:

There are two main challenges in minimizing the loss in (4). First, the output space of the selector function ({0, 1} d ) is large -its size increases exponentially with the dimension of the feature space; thus a complete search is impractical in high dimensional settings (and it should be noted that it is in high dimensional settings where feature selection is most necessary). Second, we do not have access to the densities p Y (•|x (S(x)) ) and p Y (y|x) required to compute (4).

3.1. LOSS ESTIMATION:

To approximate the densities in (3), we introduce a pair of functions f φ :

X * × {0, 1} d → [0, 1] c parametrized by φ and f γ : X → [0, 1] c parametrized by γ that will estimate p Y (•|x (S(x))

) and p Y (•|x) respectively.

3.1.1. PREDICTOR NETWORK:

We refer to f φ as the predictor network. This will take as input a suppressed3 feature vector x (s) and its corresponding selection vector s and will output a probability distribution (using a softmax layer) over the c-dimensional output space.

f φ is trained to minimize the cross entropy loss given by

l 1 (φ) = -E (x,y)∼p,s∼π θ (x,•) c i=1 y i log(f φ i (x (s) , s))

where y i is the ith component of the one-hot encoding of y and π θ is the distribution induced by our selector network which will be defined in the following section. f φ is implemented as a fully connected neural network4 .

3.1.2. BASELINE NETWORK:

We refer to f γ as the baseline network, which is standard in the actor-critic literature for variance reduction. f γ is implemented as a fully connected neural network and is trained to minimize

l 3 (γ) = -E (x,y)∼p c i=1 y i log(f γ i (x)) .

For fixed φ, γ we define our loss estimator, l, by

EQUATION

3.2. SELECTOR FUNCTION OPTIMIZATION:

We approximate the selector function S : X → {0, 1} d by using a single neural network, Ŝθ : X → [0, 1] d parameterized by weights θ, that outputs a probability for selecting each feature (i.e. the ith component of Ŝθ (x) will denote the probability with which we select the ith feature). The selector network induces a probability distribution over the selection space ({0, 1} d ), with the probability of a given joint selection vector s ∈ {0, 1} d being given by5 Instances are fed into the selector network which outputs a vector of selection probabilities. The selection vector is then sampled according to these probabilities. The predictor network then receives the selected features and makes a prediction and the baseline network is given the entire feature vector and makes a prediction. Each of these networks are trained using backpropagation using the real label. The loss of the baseline network is then subtracted from the prediction network's loss and this is used to update the selector network.

π θ (x, s) = Π d i=1 Ŝθ i (x) si (1 -Ŝθ i (x)) 1-si .

Using this, we define the following loss for our selector network

l 2 (θ) = E (x,y)∼p E s∼π θ (x,•) l(x, s) + λ||s|| 0 = X ×Y p(x, y)   s∈{0,1} d π θ (x, s) l(x, s) + λ||s|| 0   dxdy.

Taking the gradient of this loss with respect to θ gives us

∇ θ l 2 (θ) = X ×Y p(x, y)   s∈{0,1} d ∇ θ π θ (x, s) l(x, s) + λ||s|| 0   dxdy = X ×Y p(x, y)   s∈{0,1} d ∇ θ π θ (x, s) π θ (x, s) π θ (x, s) l(x, s) + λ||s|| 0   dxdy = X ×Y p(x, y)   s∈{0,1} d ∇ θ log π θ (x, s)π θ (x, s) l(x, s) + λ||s|| 0   dxdy = E (x,y)∼p E s∼π θ (x,•) l(x, s) + λ||s|| 0 ∇ θ log π θ (x, s) .

We update each of Ŝθ , f φ and f γ iteratively using stochastic gradient descent. Pseudo-code of INVASE is given in Algorithm 1 and a block representation of INVASE can be found in Fig. 1 .

4. EXPERIMENTS:

In this section, we quantitatively evaluate INVASE against various state-of-the-art benchmarks on both synthetic and real-world datasets. We evaluate our performance both at identifying ground truth relevance and at enhancing predictions. We compare our model with 4 global variable selection models: Knockoffs [3] , Tree Ensembles (Tree) [7] , Sequential Correlation Feature Selection (SCFS) [11] , and LASSO regularized linear model; and 3 instance-wise feature selection methods: L2X Sample a mini-batch from the dataset (x j , y j ) n mb j=1 ∼ D 5:

for j = 1, ..., n mb do 6:

Calculate selection probabilities

(p j 1 , ..., p j d ) ← Ŝθ (x j ) 7:

Sample selection vector 8:

for i = 1, ..., d do s j i ∼ Ber(p j i ) 9:

Calculate loss

lj (x j , s j ) ← - c i=1 y j i log(f φ i (x (sj ) j , s j )) - c i=1 y j i log(f γ i (x j )) 10:

Update the selector network parameters θ

θ ← θ -α 1 n mb n mb j=1 lj (x j , s j ) + λ||s j || ∇ θ log π θ (x j , s j ) 11: Update the predictor network parameters φ φ ← φ -β 1 n mb n mb j=1 c i=1 y j i × ∇ φ log(f φ i (x (sj ) j , s j )) 12:

Update the baseline network parameters γ

γ ← γ -β 1 n mb n mb j=1 c i=1 y j i × ∇ γ log(f γ i (x j ))

[4], LIME [24] , and Shapley [18] . The details of benchmark implementation can be found in the appendix. Implementation of INVASE can be found at https://github.com/jsyoon0823/ INVASE.

4.1. SYNTHETIC DATA EXPERIMENTS:



4.1.1. EXPERIMENTAL SETTINGS:

For our first set of experiments, we use the same synthetic data generation models as in L2X [4] . The input features are generated from an 11-dimensional 67 Gaussian distribution with no correlations across the features (X ∼ N (0, I)). The label Y is sampled as a Bernoulli random variable with

P(Y = 1|X) = 1 1+logit(X)

, where logit(X) is varied to create 3 different synthetic datasets:

• Syn1: exp(X 1 X 2 ) • Syn2: exp( 6 i=3 X 2 i -4) • Syn3: -10 × sin 2X 7 + 2|X 8 | + X 9 + exp(-X 10 )

In each of these datasets, the label depends on the same subset of features for every sample. To highlight the capability of INVASE to detect instance-wise dependence, we generate 3 further synthetic datasets as follows:

• Syn4: If X 11 < 0, logit follows Syn1, otherwise, logit follows Syn2.

• Syn5: If X 11 < 0, logit follows Syn1, otherwise, logit follows Syn3.

• Syn6: If X 11 < 0, logit follows Syn2, otherwise, logit follows Syn3.

Note that in Syn4 and Syn5, the number of relevant features is different for different samples.

For each of Syn1 to Syn6 we draw 20,000 samples from the data generation model and separate each into training (D train = (x i , y i ) 10000 i=1 ) and testing (D test = (x j , y j ) 10000 j=1 ) sets. For each method we try to find the top k relevant features for each sample (we set k = 4 for Syn1, Syn2, Syn3, Syn4, Syn5 and k = 5 for Syn6), note, however, that k is not given as an input to INVASE (but is necessary for other methods). The performance metrics we use are the true positive rate (TPR) (higher is better) and false discovery rate (FDR) 8 (lower is better) to measure the performance of the methods when the focus is on discovery (i.e. discovering which features are relevant) and we use Area Under the Receiver Operating Characteristic Curve (AUROC), Area Under the Precision Recall Curve (AUPRC) and accuracy when the focus is on predictions. The results for Syn4, Syn5 and Syn6 demonstrate that INVASE is capable of detecting a different number of relevant features for each sample when necessary -the performance improvement over L2X is greater in Syn4 and Syn5 than Syn6. In particular, in Syn4, L2X is forced to overselect features when X 11 < 0 and underselect when X 11 ≥ 0 thus resulting in higher FDR and lower TPR, respectively. To highlight this, in Table 2 we report the group specific FDR and TPR on Syn4 and Syn5 when setting k = 3, 4, 5, where Group 1 refers to samples with X 11 < 0 and Group 2 to samples with X 11 ≥ 0.

For k = 3 in Syn4, we see that INVASE and L2X have comparable FDR in Group 1, since the total number of relevant features for each sample is 3 (X 1 , X 2 , X 11 ). However, when we increase k, we see that the FDR increases for L2X as it is forced to select more than 3 features, which necessarily means that the FDR must be at least 40% even if L2X was finding the relevant features perfectly. On the other hand, for Group 2 we see that the TPR is low for k = 3 since necessarily, L2X cannot possibly select all of the 5 relevant features. INVASE, however, is able to select the correct number in both and hence enjoys low FDR and high TPR.

Syn5 reinforces the conclusions we drew for L2X in Syn4. Interestingly, though, for INVASE, we found that X 11 was almost never selected for Group 1 in Syn5. We believe this is because the lack of overlap between the relevant features for each group means that the predictor network can essentially learn two separate networks -one for each group. This is because it is possible to create two subnetworks with non-overlapping weights that each take as input the features of a given group. X 11 is therefore unnecessary for prediction. Note, however, that X 11 is highly relevant for the selector network in deciding which features to pass on and so it is not true that X 11 isn't relevant, but simply that the selector network does not need to "pass on" its relevance to the predictor network.

To investigate this further, results for settings where the features overlap between groups (and so it is not possible to disentangle the networks) can be found in the Appendix. In this experiment we analyze the effect of using feature selection as a pre-processing step for prediction. We first perform feature selection (either instance-wise or global) and then train a 3-layer fully connected network with Batch Normalization [12] in every layer (to avoid overfitting) to perform predictions on top of the (feature-selected) data. In this setting we compare the two global feature selection methods (LASSO and Tree) and one instance-wise feature selection method (L2X). Furthermore, we also compare with the predictive model without any feature selections (w/o FS) and the predictive model with ground truth globally relevant features 9 (with Global). In particular, this allows us to demonstrate that the improvements in prediction performance are not just because the global feature selection performed implicitly by INVASE is better than the other global feature selection methods but are also due to the fact that we select features on an instance-wise level. Experiments here are conducted on synthetic data with 100 features but the same labelling procedures as above.

Datasets:

As can be seen in Table 3 , there is a significant performance improvement when discarding all of the irrelevant features (with Global). However, neither of the global feature selection methods (Tree and Lasso) are capable of achieving this improvement. On the other hand, INVASE is capable of achieving (and beating -in Syn4 and Syn6) this improvement, demonstrating its capability both at selecting features globally better than existing methods but also at improving on global selection with instance-wise selection (where relevant), to provide further improvements. On the other hand, L2X performs worse than the global methods in Syn1-3, demonstrating an inability to perform even global feature selection in this higher dimensional setting (this is supported by the high dimensional discovery results in the Appendix), and in Syn4-6 is performing worse than with Global (which now is not even optimal).

Furthermore, even though we include Batch Normalization to avoid overfitting, with a small number of samples and high number of dimensions, the 3-layer fully connected network still suffers from overfitting as demonstrated by the significant difference in performance between w/o FS and with Global. This demonstrates the necessity of feature selection as a pre-processing step. The first experiment we carried out was to create semi-synthetic datasets by using the labelling procedures Syn1-6 from above but with the features now coming from real data (instead of being i.i.d. Gaussian). The results of this experiment can be found in the Appendix.

4.2.2. THE DISCOVERED FEATURE IMPORTANCE IN MAGGIC DATASET:

In this next experiment, we visualize the ability of INVASE to select features on an individualized level. Fig. 2 (left) shows the selection probability (given by INVASE) of each feature for 20 randomly selected patients in the MAGGIC dataset. Fig. 2 (right) shows the selection probability of each feature averaged over different binary splits of the data (i.e. when split into Male and Female). In Table 4 , we also report the mean and variance of the number of selected features in each subgroup.

As can be seen, INVASE discovers significantly different features for both individuals and for different subgroups of the dataset.

4.2.3. RESULTS: PREDICTION USING REAL DATA VARIABLES WITH REAL LABEL:

Evaluating the performance of feature selection methods on real data is difficult, since ground truth relevance is often not known. We therefore cannot use TPR and FDR to evaluate the performance on real data. In our final experiment, therefore, we instead focus on prediction performance exactly as in 4.1.3 (except now both the features and label come from real data). As can be seen in Table 5 , INVASE consistently improves prediction performance in each of the two settings (different time horizons) in each dataset.

5. FUTURE WORK:

While this paper has focused on discovering relevant features in the static setting, this could also be extended to apply in the temporal setting. One such avenue of exploration for this would be to replace each of the networks with an RNN. Particular care will need to be taken in defining the problem, though; do we treat each stream as a feature or each time point of each stream? We leave this investigation to future work.

HIGH DIMENSIONAL DISCOVERY:

To demonstrate the scalability of our method, we run an experiment in which we increase the total number of features to 100. The features are generated as a 100-dimensional Gaussian with no correlations (N (0, I)) and the relationships between features and label remains as in Table 1 As can be seen in Table 7 , INVASE also works consistently better than all other benchmarks in all 6 synthetic datasets in this setting. In fact, we see a significant reduction in performance (compared to the 11 feature setting) for L2X in Syn1, with the TPR dropping more than 90% leading to an almost complete failure of the method to detect any relevant features. In particular, we see that L2X does not scale as well as INVASE with the dimensionality of the data, which is particularly limiting for a feature selection method.

We also compare the CPU times of the algorithm for training and testing with other instance-wise feature selection benchmarks to show the scalability in terms of computational complexity. As can be seen in 

ADDITIONAL RESULTS ON COMPLEX SYNTHETIC DATASETS:

In the main paper, the relevant subset for a sample in each of our variable synthetic datasets (Syn4-6) depended on X 11 only, which was unused in the rest of the model (i.e. X 11 determined only the relevant subset, and was otherwise unused as a predictive variable). In this set of experiments, we investigate the effect of having the subset relevance depend on a variable that is also used in the model itself (Syn4A, Syn5A, Syn6A). We then investigate the effect of having more than one variable being used to determine subset relevance (Syn4B, Syn5B, Syn6B, Syn7). The results for these are reported in Tables 10 and 11 , respectively.

The input features are generated from a 100-dimensional Gaussian distribution with no correlations across the features (X ∼ N (0, I)). Y is generated according to P(Y = 1|X) = 1 1+logit(X) with the logit value for each synthetic dataset now defined as follows:

• Syn4A: If X 1 < 0, logit = exp(X 1 X 2 ), otherwise, logit =exp( 6 i=3 X 2 i -4). • Syn5A: If X 1 < 0, logit = exp(X 1 X 2 ), otherwise, logit =-10 × sin 2X 7 + 2|X 8 | + X 9 +

exp(-X 10 ).

• Syn6A: If X 7 < 0, logit =exp( 

• Syn4B: If X 1 X 3 < 0, logit = exp(X 1 X 2 ), otherwise, logit =exp( 6 i=3 X 2 i -4). • Syn5B: If X 1 X 7 < 0, logit = exp(X 1 X 2 ), otherwise, logit =-10 × sin 2X 7 + 2|X 8 | + X 9 + exp(-X 10 ). • Syn6B: If X 3 X 7 < 0, logit =exp( 6 i=3 X 2 i -4), otherwise, logit =-10 × sin 2X 7 + 2|X 8 | + X 9 + exp(-X 10 ). • Syn7: -If X 1 < 0, X 2 < 0, logit = exp(X 1 X 2 ) -If X 1 < 0, X 2 ≥ 0, logit =exp( 6 i=3 X 2 i -4). -If X 1 ≥ 0, X 2 < 0, logit =-10 × sin 2X 7 + 2|X 8 | + X 9 + exp(-X 10 ). -If X 1 ≥ 0, X 2 ≥ 0, logit =0.5 × exp(X 1 X 2 ) + 0.5 × exp( 4 i=3 X 2 i -2).

RESULTS ON SEMI-SYNTHETIC DATASETS:

In this experiment, we use real features (which have correlation across features) but generate the labels as in the synthetic experiments from the main paper, using Syn1-Syn6. This allows us to know the ground truth relevance of the features, and calculate TPR and FDR, while using unknown and correlated feature distributions (instead of the unrealistic setting of i.i.d. Gaussian used in the fully synthetic experiment). The results for the MAGGIC and PLCO datasets are given below. We do the same analysis for the MAGGIC dataset; results are given in Table 16 . We see that here the linear correlation with the label is stronger and this is reflected in Tables 12 and 13 , where all of the linear models performed better than in the fully-synthetic settings. However, we note that although they had a better performance, in most cases it was still not comparable with INVASE.

COMPUTER VISION:

Another natural application of INVASE is in computer vision. To briefly demonstrate the applicability and capability of INVASE to computer vision, we conduct two experiments using the Kaggle Dogs vs. Cats dataset (https://www.kaggle.com/c/dogs-vs-cats) [6] and the Oxford Pet dataset (http://www.robots.ox.ac.uk/ ˜vgg/data/pets/) [20] . The is to select a set 16 x 16 patches of each image that maximize the predictive capability of a model. In order to apply INVASE to this problem, we simply treat each 16 x 16 patch as a feature.

We use the U-Net [25] architecture for the selector network and the VGG network [28] architecture for the predictor and baseline networks. Below we give qualitative results of INVASE applied to these datasets, where we see that INVASE successfully identifies patches of each image in which the animal's face is visible. 

:

In this paper we focus on classification; we discuss an extension of our model to regression in the Appendix.

We will occasionally abuse notation and write yi to denote the ith element of the one-hot encoding of y, though the context should make it clear when this is the case.

When implemented we set * = 0 and include the selection vector to differentiate this from the case xi = 0.

f φ , f γ and Ŝθ could also be implemented as CNNs or RNNs, when appropriate.

Note that, when d is large, this becomes vanishingly small, however, π θ appears in our loss only via its log and so in practice this is not a problem.

In L2X they use a 10-dimensional Gaussian, we introduce X11 to act as a "switch" to create instance-wise relevance. Experiments where instead the "switch" variable is one of X1, ..., X10 can be found in the appendix.

We also perform experiments using 100 features in the Appendix to demonstrate the scalability of our method.

Definitions of TPR and FDR can be found in the Appendix.

For example, in Syn1 the predictor network in the with Global setting is trained on only X1 and X2 and in Syn4 it would be trained on X1, X2, X3, X4, X5, X6, X11.

Figures and Tables:

Figure 1: Block diagram of INVASE.Instances are fed into the selector network which outputs a vector of selection probabilities. The selection vector is then sampled according to these probabilities. The predictor network then receives the selected features and makes a prediction and the baseline network is given the entire feature vector and makes a prediction. Each of these networks are trained using backpropagation using the real label. The loss of the baseline network is then subtracted from the prediction network's loss and this is used to update the selector network.

Figure 2: Left: The feature importance for each of 20 randomly selected patients in the MAGGIC dataset. Right: The average feature importance for different binary splits in the MAGGIC dataset.

Figure 4: Selected 16 x 16 patches by INVASE on Kaggle Dogs vs. Cats dataset -Cats

Figure 5: Selected 16 x 16 patches by INVASE on Kaggle Dogs vs. Cats dataset -Dogs

Figure 6: Selected 16 x 16 patches by INVASE on Oxford Pet dataset -Cats

Figure 7: Selected 16 x 16 patches by INVASE on Oxford Pet dataset -Dogs

Algorithm 1 Pseudo-code of INVASE

Table content: <table/>

Relevant feature discovery results for Synthetic datasets with 11 featuresAs demonstrated by Table1, our method is capable of detecting relevant features on a global level (Syn1, Syn2 and Syn3) as well as on an instance-wise level (Syn4, Syn5 and Syn6) outperforming all other methods in both cases (both global and instance-wise methods). The particularly poor performance of some global feature selection methods in Syn1, Syn2 and Syn3 (where there is no instance-wise relevance) is due to the non-linearity of the relationship between features and labels, further details can be found in the Appendix.

Table content: <table><tr><td colspan="2">4.1.2 DISCOVERY</td><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Dataset</td><td colspan="2">Syn1</td><td colspan="2">Syn2</td><td colspan="2">Syn3</td><td>Syn4</td><td colspan="2">Syn5</td><td colspan="2">Syn6</td></tr><tr><td colspan="12">Metrics (%) TPR FDR TPR FDR TPR FDR TPR FDR TPR FDR TPR FDR</td></tr><tr><td>INVASE</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td>92.0</td><td>0.0</td><td colspan="2">99.8 10.3 84.8</td><td>1.1</td><td>90.1</td><td>7.4</td></tr><tr><td>L2X</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td>69.4</td><td>30.6</td><td colspan="5">79.5 21.8 74.8 26.3 83.3 16.7</td></tr><tr><td>LIME</td><td>13.8</td><td colspan="2">86.2 100.0</td><td>0.0</td><td>98.1</td><td>1.9</td><td colspan="5">40.7 49.4 41.1 50.6 50.5 49.5</td></tr><tr><td>Shapley</td><td>60.4</td><td>39.6</td><td>93.3</td><td>6.7</td><td>90.9</td><td>9.1</td><td colspan="5">65.2 31.9 62.9 33.7 71.2 28.8</td></tr><tr><td>Knockoff</td><td>10.0</td><td>70.0</td><td>8.7</td><td>36.2</td><td>81.2</td><td>17.5</td><td colspan="5">38.8 35.1 41.0 51.1 56.6 42.1</td></tr><tr><td>Tree</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td colspan="5">54.7 39.0 56.8 37.5 60.0 40.0</td></tr><tr><td>SCFS</td><td>23.5</td><td>76.5</td><td>39.5</td><td>60.5</td><td>78.3</td><td>22.0</td><td colspan="5">48.9 52.4 42.4 51.2 56.1 43.9</td></tr><tr><td>LASSO</td><td>19.0</td><td>81.0</td><td>39.8</td><td>60.2</td><td>78.3</td><td>21.7</td><td colspan="5">49.9 50.9 45.5 48.2 56.4 43.6</td></tr></table>



Table content: <table><tr><td/><td/><td>Syn4</td><td/><td/><td colspan="2">Syn5</td></tr><tr><td>Group</td><td>1</td><td>2</td><td/><td>1</td><td/><td>2</td></tr><tr><td colspan="8">Metrics (%) TPR FDR TPR FDR TPR FDR TPR FDR</td></tr><tr><td>INVASE</td><td colspan="2">99.5 24.6 100.0</td><td>0.4</td><td>69.2</td><td>1.6</td><td>99.8</td><td>0.6</td></tr><tr><td colspan="2">L2X (k = 3) 71.1 28.9</td><td>57.2</td><td>4.6</td><td colspan="3">65.5 34.5 55.4</td><td>7.7</td></tr><tr><td colspan="2">L2X (k = 4) 81.0 39.2</td><td>74.9</td><td>6.3</td><td colspan="3">76.2 42.9 72.4</td><td>9.4</td></tr><tr><td colspan="2">L2X (k = 5) 89.9 46.0</td><td>84.6</td><td>15.4</td><td colspan="4">87.5 47.5 82.1 17.9</td></tr></table>

Prediction performance comparison with and without feature selection methods (L2X, LASSO, Tree, INVASE, and Global). Global is using ground-truth globally relevant features for each dataset titatively shows that instance-wise feature selection can further improves the predictive model from ground truth global feature selection.

Table content: <table><tr><td>Dataset</td><td/><td/><td colspan="2">AUROC</td></tr><tr><td/><td>w/o FS</td><td colspan="2">with Global with INVASE</td><td>with Tree</td><td>with L2X with LASSO</td></tr><tr><td>Syn1</td><td colspan="2">.578±.004 .686±.005</td><td>.690±.006</td><td colspan="2">.574±.101 .498±.005</td><td>.498±.006</td></tr><tr><td>Syn2</td><td colspan="2">.789±.003 .873±.003</td><td>.877±.003</td><td colspan="2">.872±.003 .823±.029</td><td>.555±.061</td></tr><tr><td>Syn3</td><td colspan="2">.854±.004 .900±.003</td><td>.902±.003</td><td colspan="2">.899±.001 .862±.009</td><td>.886±.003</td></tr><tr><td>Syn4</td><td colspan="2">.558±.021 .774±.006</td><td>.787±.004</td><td colspan="2">.684±.017 .678±.024</td><td>.514±.031</td></tr><tr><td>Syn5</td><td colspan="2">.662±.013 .784±.005</td><td>.784±.005</td><td colspan="2">.741±.004 .709±.008</td><td>.691±.024</td></tr><tr><td>Syn6</td><td colspan="2">.692±.015 .858±.004</td><td>.877±.003</td><td colspan="2">.771±.031 .827±.017</td><td>.727±.025</td></tr><tr><td>Dataset</td><td/><td/><td colspan="2">AUPRC</td></tr><tr><td/><td>w/o FS</td><td colspan="2">with Global with INVASE</td><td>with Tree</td><td>with L2X with LASSO</td></tr><tr><td>Syn1</td><td colspan="2">.567±.007 .690±.006</td><td>.694±.006</td><td colspan="2">.577±.102 .498±.007</td><td>.499±.008</td></tr><tr><td>Syn2</td><td colspan="2">.799±.005 .878±.005</td><td>.886±.004</td><td colspan="2">.878±.004 .817±.031</td><td>.591±.037</td></tr><tr><td>Syn3</td><td colspan="2">.861±.003 .905±.002</td><td>.907±.003</td><td colspan="2">.904±.002 .860±.012</td><td>.890±.002</td></tr><tr><td>Syn4</td><td colspan="2">.572±.019 .794±.006</td><td>.804±.004</td><td colspan="2">.681±.031 .672±.025</td><td>.536±.025</td></tr><tr><td>Syn5</td><td colspan="2">.665±.019 .796±.005</td><td>.797±.006</td><td colspan="2">.765±.003 .719±.011</td><td>.680±.040</td></tr><tr><td>Syn6</td><td colspan="2">.709±.018 .870±.005</td><td>.886±.004</td><td colspan="2">.779±.027 .835±.017</td><td>.757±.036</td></tr></table>

Selection probability of overall and patient subgroups by INVASE in MAGGIC dataset. (Mean ± Std)

Table content: <table><tr><td>Overall</td><td>Male</td><td>Diabetes</td><td/><td colspan="2">Hypertension</td><td>Smoker</td><td>Heart Failure</td></tr><tr><td colspan="2">43.5±10.7</td><td colspan="2">53.2±10.8</td><td colspan="2">46.6±9.3</td><td>41.0±12.1</td><td>51.8±11.1</td></tr><tr><td>42.5±18.4</td><td>Female</td><td colspan="5">Non-diabetes Non-hypertension Non-smoker No Heart Failure</td></tr><tr><td colspan="2">40.8±15.6</td><td colspan="2">39.3±8.0</td><td colspan="2">40.0±9.3</td><td>43.2±7.0</td><td>39.6±6.9</td></tr><tr><td>Datasets</td><td colspan="2">Metrics</td><td/><td>AUROC</td><td>AUPRC</td><td>AUROC</td><td>AUPRC</td></tr><tr><td/><td colspan="2">Labels</td><td/><td colspan="2">3 year</td><td>5 year</td></tr><tr><td>MAGGIC</td><td colspan="2">INVASE</td><td colspan="4">.722±.005 .655±.010 .740±.005 .867±.006</td></tr><tr><td/><td colspan="6">Without INVASE .720±.006 .639±.009 .730±.006 .855±.004</td></tr><tr><td/><td colspan="2">Labels</td><td/><td colspan="2">5 year</td><td>10 year</td></tr><tr><td>PLCO</td><td colspan="2">INVASE</td><td colspan="4">.637±.007 .329±.013 .673±.007 .506±.006</td></tr><tr><td/><td colspan="6">Without INVASE .629±.008 .324±.011 .657±.006 .485±.008</td></tr></table>

Prediction performance for MAGGIC and PLCO dataset.

Table content: <table/>

in the main manucript (i.e. we are adding 89 additional noisy signals that have no effect on the label). Relevant feature discovery for synthetic datasets with 100 features

Table content: <table><tr><td>Dataset</td><td colspan="2">Syn1</td><td colspan="2">Syn2</td><td colspan="2">Syn3</td><td>Syn4</td><td>Syn5</td><td>Syn6</td></tr><tr><td colspan="2">Metrics (%) TPR</td><td>FDR</td><td colspan="6">TPR FDR TPR FDR TPR FDR TPR FDR TPR FDR</td></tr><tr><td>INVASE</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td colspan="2">66.3 40.5 73.2 23.7 90.5 15.4</td></tr><tr><td>L2X</td><td>6.1</td><td>93.9</td><td>81.4</td><td>18.6</td><td>57.7</td><td>42.3</td><td colspan="2">48.5 46.4 35.4 60.8 66.3 33.7</td></tr><tr><td>LIME</td><td>0.0</td><td colspan="2">100.0 100.0</td><td>0.0</td><td>92.7</td><td>7.3</td><td colspan="2">43.8 47.4 42.3 50.1 50.1 49.9</td></tr><tr><td>Shapley</td><td>4.4</td><td>95.6</td><td>95.1</td><td>4.9</td><td>88.8</td><td>11.2</td><td colspan="2">50.2 43.4 49.9 44.2 62.5 37.5</td></tr><tr><td>Knock off</td><td>0.0</td><td>64.9</td><td>3.7</td><td>71.2</td><td>74.9</td><td>24.9</td><td colspan="2">28.2 59.8 33.1 59.4 46.9 53.0</td></tr><tr><td>Tree</td><td>49.9</td><td colspan="2">50.1 100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td colspan="2">40.7 49.5 56.7 37.5 58.4 41.6</td></tr><tr><td>SCFS</td><td>2.5</td><td>97.5</td><td>5.3</td><td>94.7</td><td>74.9</td><td>25.1</td><td colspan="2">27.0 74.6 30.6 62.1 38.3 61.7</td></tr><tr><td>LASSO</td><td>2.5</td><td>97.5</td><td>4.0</td><td>96.0</td><td>75.3</td><td>24.7</td><td colspan="2">28.3 73.2 36.0 56.9 45.9 54.1</td></tr></table>

INVASE is much faster (10 times) than LIME and Shapley methods and comparable with L2X; we see that INVASE takes approximately 50% longer to run than L2X, which can be accounted for by the addition of a 3rd network (the baseline network) in INVASE that is not present in L2X. Note, however, that this baseline network can be trained in parallel with the predictor network and we believe that doing so would lead to both INVASE and L2X having the same run-time.

Table content: <table><tr><td colspan="2">Methods INVASE</td><td>L2X</td><td>Shapley</td><td>LIME</td></tr><tr><td>Train</td><td colspan="3">1327.69s 939.82s 12801.21s</td><td>-</td></tr><tr><td>Test</td><td>0.38s</td><td>0.78s</td><td>0.06s</td><td>18931.98s</td></tr></table>

Comparison of CPU clock time across different instance-wise feature selection methods on average across Syn1 to Syn6 with 100 features and 10,000 samples on training/testing, respectivelyHYPER-PARAMETER ANALYSISIn the following experiment, we provide results for various values of the hyper-parameter, λ, in the Syn4, Syn5, and Syn6 100-dimensional setting. Table9gives the results in terms of TPR and FDR. Note that in the other experiments, we select the hyper-parameter λ which maximizes the predictor accuracy in terms of AUROC.

Table content: <table><tr><td>Datasets</td><td colspan="2">Syn4</td><td colspan="2">Syn5</td><td colspan="2">Syn6</td></tr><tr><td colspan="7">λ / Metris (%) TPR FDR TPR FDR TPR FDR</td></tr><tr><td>0.1</td><td colspan="6">98.0 94.3 90.0 93.4 99.2 92.3</td></tr><tr><td>0.3</td><td colspan="6">93.7 87.9 84.2 88.9 96.9 86.7</td></tr><tr><td>0.5</td><td colspan="6">99.0 43.1 88.3 50.6 99.6 31.7</td></tr><tr><td>1</td><td colspan="6">66.3 40.5 73.2 23.7 90.5 15.4</td></tr><tr><td>2</td><td>0.0</td><td>0.0</td><td>25.4</td><td>4.1</td><td>67.1</td><td>3.6</td></tr><tr><td>5</td><td>0.0</td><td>0.0</td><td>7.5</td><td>2.7</td><td>7.6</td><td>2.5</td></tr><tr><td>10</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr></table>

Relevant feature discovery results for various values of the hyper-parameter λ in the Syn4, Syn5, and Syn6 100-dimensional setting.

Table content: <table/>

+2|X 8 |+ X 9 + exp(-X 10 ). Relevant feature discovery results for complex synthetic datasets (Syn4A, 5A, 6A) with 100 features

Table content: <table><tr><td colspan="5">6 i=3 X 2 i -4), otherwise, logit =-10×sin 2X 7 Dataset Syn4A Syn5A Syn6A</td></tr><tr><td colspan="5">Metrics (%) TPR FDR TPR FDR TPR FDR</td></tr><tr><td>INVASE+</td><td>77.5 14.5 85.9</td><td>8.8</td><td>89.9</td><td>7.3</td></tr><tr><td>L2X</td><td colspan="4">65.0 39.3 48.0 57.4 74.4 35.5</td></tr><tr><td>LIME</td><td colspan="4">56.3 49.2 58.2 48.8 58.9 47.8</td></tr><tr><td>Shapley</td><td colspan="4">71.8 39.8 71.0 41.3 68.9 38.2</td></tr><tr><td>Knock off</td><td colspan="4">59.8 62.6 55.0 49.9 65.0 40.0</td></tr><tr><td>Tree</td><td colspan="4">61.3 46.9 75.6 39.4 66.9 40.0</td></tr><tr><td>SCFS</td><td colspan="4">52.8 66.9 55.3 50.6 50.4 51.8</td></tr><tr><td>LASSO</td><td colspan="4">61.0 61.2 55.0 50.0 53.9 48.8</td></tr></table>

Correlation between features and labels in Synthetic datasets with 100 features. Ground truth (in the global sense) relevant features are given in bold. Features with correlation > 0.05 are highlighted in red.As can be seen in Table15, among 33 relevant features, only 9 features have more than 0.05 (linear) correlation with the label. In particular, using a linear model, it is very hard to discover the relevant features. However, Knock-off (based on LASSO and linear correlations), LASSO, and SCFS are linear models, resulting in a poor performance in our experiments. The above table results are directly reflected in the results given in the main manuscript.

Table content: <table><tr><td colspan="2">CORRELATIONS BETWEEN FEATURES AND LABELS IN THE SYNTHETIC AND</td></tr><tr><td colspan="2">SEMI-SYNTHETIC EXPERIMENTS</td></tr><tr><td colspan="2">Variables Syn1 Syn2 Syn3 Syn4 Syn5 Syn6</td></tr><tr><td>X 1</td><td>0.003 0.008 0.006 0.009 0.007 0.006</td></tr><tr><td>X 2</td><td>0.001 0.005 0.006 0.005 0.015 0.005</td></tr><tr><td>X 3</td><td>0.006 0.011 0.001 0.017 0.016 0.010</td></tr><tr><td>X 4</td><td>0.006 0.003 0.003 0.002 0.000 0.002</td></tr><tr><td>X 5</td><td>0.003 0.015 0.022 0.004 0.017 0.028</td></tr><tr><td>X 6</td><td>0.003 0.004 0.005 0.002 0.004 0.005</td></tr><tr><td>X 7</td><td>0.013 0.009 0.481 0.002 0.242 0.235</td></tr><tr><td>X 8</td><td>0.010 0.008 0.012 0.003 0.010 0.022</td></tr><tr><td>X 9</td><td>0.001 0.003 0.239 0.002 0.115 0.121</td></tr><tr><td>X 10</td><td>0.002 0.003 0.308 0.003 0.149 0.144</td></tr><tr><td>X 11</td><td>0.014 0.012 0.004 0.028 0.018 0.002</td></tr><tr><td colspan="2">Variables Syn1 Syn2 Syn3 Syn4 Syn5 Syn6</td></tr><tr><td>X 1</td><td>0.028 0.030 0.070 0.011 0.044 0.026</td></tr><tr><td>X 2</td><td>0.002 0.011 0.009 0.001 0.001 0.012</td></tr><tr><td>X 3</td><td>0.018 0.079 0.008 0.038 0.006 0.046</td></tr><tr><td>X 4</td><td>0.005 0.113 0.006 0.056 0.001 0.055</td></tr><tr><td>X 5</td><td>0.006 0.034 0.032 0.013 0.016 0.036</td></tr><tr><td>X 6</td><td>0.019 0.114 0.027 0.099 0.018 0.004</td></tr><tr><td>X 7</td><td>0.005 0.010 0.367 0.000 0.262 0.272</td></tr><tr><td>X 8</td><td>0.020 0.030 0.112 0.023 0.075 0.082</td></tr><tr><td>X 9</td><td>0.025 0.022 0.299 0.006 0.216 0.200</td></tr><tr><td>X 10</td><td>0.07 0.043 0.328 0.027 0.222 0.206</td></tr><tr><td>X 11</td><td>0.009 0.006 0.034 0.046 0.018 0.058</td></tr></table>

Correlation between features and labels in MAGGIC datasets. Ground truth relevant features are described in bold. Features with correlation > 0.05 are described in red

Table content: <table/>

Peer review (iclr2019):
 This paper proposes a new instance-wise feature selection method, INVASE. It is closely related to the prior work L2X (Learning to Explain). There are three differences compared to L2X. The most important difference is about how to backpropagate through subset sampling to select features. L2X use the Gumbel-softmax trick and this paper uses actor-critic models. The paper is written well. It is easy to follow the paper. The contribution of this paper is that it provides a new way, compared to L2X, to backpropagate through subset sampling in order to select features. The authors compare INVASE with L2X and several other approaches on synthetic data and show outperforming results. In the real-world experiments, the authors do not compare INVASE with other approaches. Regarding experiments, instance-wise feature selection is often applied on computer vision or natural language process applications, where global feature selection is not enough. This paper lacks experiments on CV or NLP applications. For the MAGGIC dataset, I expect to see subgroup patterns. The patterns that authors show in Figure 2 are very different for all randomly selected 20 patients. The authors do not explain why it is preferred to see very different feature patterns for all patients instead of subgroup patterns. I have questions about other two differences from L2X, pointed by the authors. First, the selector function outputs a probability for selecting each feature .hat{S}^.theta(x). In the paper of L2X, it also produces a weight vector w_.theta(x) as described in section 3.4. I think the .hat{S}^.theta(x) has similar meaning as w_.theta(x) in L2X. In the synthetic data experiment, the authors fix the number of selected features for L2X so that it forces to overselect or underselect features in the example of Syn4. Did the author try to relax this constraint for L2X and use w_.theta(x) in L2X to select features as using .hat{S}^.theta(x) in INVASE? Second, I agree with the authors that L2X is inspired by maximizing mutual information between Y and X_S and INVASE is inspired by minimizing KL divergence between Y|X and Y|X_S. Both intuitions lead to similar objective functions that INVASE has an extra term .log p(y|x) and .lambda ||S(x)||. INVASE is able to add a l_0 penalty on S(x) since it uses the actor-critic models. For the .log p(y|x) term, as the author mentioned, it helps to reduce the variance in actor-critic models. This .log p(y|x) term is a constant in the optimization of S(x). In Algorithm 1, 12, the updates of .gamma does not depend on other parameters related to the predictor network and selector network. Could the authors first train a baseline network and use it as a fixed function in Algorithm 1? I don*t understand the meaning of updates for .gamma iteratively with other parameters since it does not depend on the learning of other parameters. Does this constant term .log p(y|x) have other benefits besides reducing variance in actor-critic models? I have another minor question about scaling. How does the scaling of X affect the feature importance learned by INVASE? Note: I have another concern about the experiments. Previous instance-wise variable selection methods are often tested on CV or NLP applications, could the authors present those experiments as previous works?