How good Neural Networks interpretation methods really are? A
  quantitative benchmark
Saliency Maps (SMs) have been extensively used to interpret deep learning
models decision by highlighting the features deemed relevant by the model. They
are used on highly nonlinear problems, where linear feature selection (FS)
methods fail at highlighting relevant explanatory variables. However, the
reliability of gradient-based feature attribution methods such as SM has mostly
been only qualitatively (visually) assessed, and quantitative benchmarks are
currently missing, partially due to the lack of a definite ground truth on
image data. Concerned about the apophenic biases introduced by visual
assessment of these methods, in this paper we propose a synthetic quantitative
benchmark for Neural Networks (NNs) interpretation methods. For this purpose,
we built synthetic datasets with nonlinearly separable classes and increasing
number of decoy (random) features, illustrating the challenge of FS in
high-dimensional settings. We also compare these methods to conventional
approaches such as mRMR or Random Forests. Our results show that our simple
synthetic datasets are sufficient to challenge most of the benchmarked methods.
TreeShap, mRMR and LassoNet are the best performing FS methods. We also show
that, when quantifying the relevance of a few non linearly-entangled predictive
features diluted in a large number of irrelevant noisy variables, neural
network-based FS and interpretation methods are still far from being reliable.
Antoine Passemiers, Pietro Folco, Daniele Raimondi, Giovanni Birolo, Yves Moreau, Piero Fariselli
# How good Neural Networks interpretation methods really are? A quantitative benchmark.

###### Abstract

Saliency Maps (SMs) have been extensively used to interpret deep learning models decision by highlighting the features deemed relevant by the model. They are used on highly nonlinear problems, where linear feature selection (FS) methods fail at highlighting relevant explanatory variables. However, the reliability of gradient-based feature attribution methods such as SM has mostly been only qualitatively (visually) assessed, and quantitative benchmarks are currently missing, partially due to the lack of a definite ground truth on image data. Concerned about the appohenic biases introduced by visual assessment of these methods, in this paper we propose a synthetic quantitative benchmark for Neural Networks (NNs) interpretation methods. For this purpose, we built synthetic datasets with nonlinearly separable classes and increasing number of _decoy_ (random) features, illustrating the challenge of FS in high-dimensional settings. We also compare these methods to conventional approaches such as mRMR or Random Forests. Our results show that our simple synthetic datasets are sufficient to challenge most of the benchmarked methods. TreeShap, mRMR and LassoNet are the best performing FS methods. We also show that, when quantifying the relevance of a few non linearly-entangled predictive features diluted in a large number of irrelevant noisy variables, neural network-based FS and interpretation methods are still far from being reliable.

## 1 Introduction

Decision processes are different in machines and humans. The training of Neural Networks (NNs) is solely guided by the minimizing of a loss function, which incidentally and only partially aligns with human understanding of the problem, since the optimization of the loss function relies on the labeling of training samples rather than articulate expert insights, eventually making these models diverge from human expectations. An extreme example of such divergence are _Clever Hans_ predictors, which exclusively rely on dataset artifacts [1]. More often, their decision processes are influenced by contextual features which are not necessary causal, and which are thus acting as confounders (e.g. presence of wooden hurdles in horse pictures). Therefore, _Clever Hans_ predictors can appear proficient on training data but remain prone to potentially high generalization error in unseen settings or on independent test data (e.g. absence of wooden hurdles in horse pictures).

For these reasons, interpretability has become an emerging crucial aspect of Machine Learning (ML), and solutions for Explainable Artificial Intelligence (XAI) [1] are now also ethical requirements posed by institutions like the European Union [2]. In many _real-life_ applications, a clear understanding of the model decision-making process is indeed necessary for its adoption in a production environment. Interpretability is even more relevant in the context of nonlinear techniques with strong prediction power and modeling capabilities such as Deep Neural Networks (DNN) [1], since linear models are notoriously easier to explain and can therefore be complemented with simple but sound and effective FS methods (LASSO [3], shrunken centroid method [4]) that exploits the additivity of input variables [5].

Due to the intrinsic human-readability properties of simpler models (i.e. linear models, decision trees), they are often preferred in critical settings. On the other hand, complex nonlinear models such as DNNs exploit non-trivial
interactions between input features, making them a model-of-choice for solving difficult tasks that require higher abstraction (e.g. Genome interpretation, protein folding and human-level real-time control [6; 7; 8; 9]). In theory, DNNs can ignore irrelevant features during the training phase, but this depends on the optimizer and the loss function landscape, which can be highly non-convex.

In recent years, gradient-based _a posteriori_ interpretation methods for DNNs such as Saliency Maps (SMs) rapidly gained popularity [10]. The most common approaches include Integrated Gradients [11], DeepLift [12], Input \(\times\) Gradient [13], SmoothGrad [14] and Guided Backpropagation [15]. They can be easily applied to any NN model with libraries such as Captum [16].

These methods have been developed mainly in the context of computer vision, and their ability to identify the _salient_ pixels for the classification of images was primarily qualitatively assessed by visually inspecting the obtained SMs in relation to the input images, because there was no _ground truth_ for what should be considered salient, since this concept is intrinsically model-dependent. To the best of our knowledge, quantitative benchmarks of the reliability of SM interpretations is currently missing, especially on non-image data. At the same time, several studies showing puzzling SM interpretation results and perplexing behaviors have been already published, raising concerns about these methods [17; 18; 19].

The interpretation of non-linear model is an extremely active field of research and, besides _a posteriori_ interpretation methods such as SMs, in recent years several DNNs architectures incorporating FS capabilities have been developed, such as CancelOut [20], DeepPINK [21], LassoNet [22], FSNet [23], Concrete Autoencoder [24] and Diet-Net [25].

In this paper, we benchmarked the most recent nonlinear FS and interpretation methods for DNNs on synthetic, simplistic, yet challenging artificial datasets, comparing them with older and more conventional FS approaches. Each of our datasets consists of few (2-7) variables that jointly and nonlinearly correlate with the output class, as well as a variable number of irrelevant random features (decoys). To allow fair assessment, both relevant and irrelevant variables have been sampled with the exact same variance: each feature follows a uniform marginal distribution in 4 out of 5 of our datasets. In our fifth dataset, features have been simply standardized. In this way, methods that exploit variance signatures, such as Principal Component Analysis are not applicable. These datasets have been constructed in such a way that the classes cannot be segregated by linear decision boundaries, making linear FS methods totally unsuitable for the problem at hand. The synthetic nature of the datasets grants us complete knowledge of the predictive signal and allows us to quantitatively benchmark the ability of nonlinear FS methods to detect nonlinearly and jointly relevant features in controlled sample-to-feature ratio experimental settings. As baseline methods, we used additional FS methods that are not based on NNs, like Random forests (RFs) or minimum redundancy maximum relevance (mRMR).

Our results show that the DNN-based FS methods tested are not able to extract relevant features if they are diluted in a random set of noisy variables. Conversely, feature relevances extracted from RFs are more reliable on average. We obtained similar results while benchmarking several SM a posteriori interpretation methods. These results indicate that the field of FS and DNNs interpretation needs to further refine the available methods, and suggest that a standardized _quantitative validation_ of the newly proposed methods should be used to assess their performance and limits, to give users a more realistic idea of the situations in which these approaches are actually reliable.

## 2 Methods

### Benchmark datasets

#### 2.1.1 Uniformly-distributed variables

We first built four datasets representing archetypal nonlinear binary classification tasks. Each dataset contains \(n=1000\) observations and \(m=p+k\) features, uniformly distributed in the \([0,1]\) interval. \(p\) and \(k\) denote the number of predictive and irrelevant features, respectively. Each dataset was then built by attributing a label to data points according to a nonlinear combination of the predictive features. All the remaining features are effectively random with respect to the labels, and thus act as decoys when it comes to feature selection. The number of positives is equal to the number of negatives, to prevent any artifact due to class imbalance. Here we describe the different characteristics of each outcome (label). They are also visually shown in Fig. 1.

* RING: positive labels are associated to the points that form a bi-dimensional ring, defined by the features in positions \(j\in\{0,1\}\). The total number of predictive features is 2. Points where assigned to the positive class when: \[|\sqrt{(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}}-0.35|\leq 0.1151\]
* XOR: the bi-dimensional space formed by the features in position \(j\in\{0,1\}\) is divided into 4 identical regions. Points in the upper left and lower right quadrants are labelled as positive samples. The total number of predictive features is 2. Points were considered positive when: \[(x_{0}-0.5)(0.5-x_{1})\geq 0\]
* RING+XOR: the samples that are either positive samples in the RING dataset (considering the features in position \(j\in\{0,1\}\)) or positive samples in the XOR dataset (considering the features in positions \(j\in\{2,3\}\)) are positive. This dataset thus contains 4 predictive features, in positions \(j\in\{0,1,2,3\}\). Points were considered positives when they satisfied any of the following: \[|\sqrt{(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}}-0.35|\leq 0.0704\] \[(x_{2}-0.5)(0.5-x_{3})\geq 0.0337\]
* RING+XOR+SUM: the samples that are positive in the RING, XOR datasets or such that \(x_{4}+x_{5}+\epsilon>0.5\) are positive samples. \(x_{j}\) denotes the \(j\)-th feature and \(\epsilon\) is Gaussian noise sampled from \(\mathcal{N}(\mu=0,\sigma=0.2)\). Like in the RING and XOR datasets, the corresponding predictive features \(\{0,1,2,3\}\) do not contain noise. In RING+XOR+SUM the predictive features are in positions \(\{0,1,2,3,4,5\}\). See Suppl. Fig. 4 for a visual explanation. Points were considered positives when they satisfied at least one of the following inequalities: \[|\sqrt{(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}}-0.35|\leq 0.0479\] \[(x_{2}-0.5)(0.5-x_{3})\geq 0.0598\] \[x_{4}+x_{5}+\mathcal{N}(\mu=0,\sigma=0.2)\geq 1.4074\]

The four outcomes are shown in Fig.1.

#### 2.1.2 Gaussian graphical model

While for sufficiently large values of \(m\) the datasets defined above are already challenging, we built an additional dataset (called DAG) characterized also by confounding effects. Indeed, these effects occur in most real-life settings, and are very likely to _misguide_ the different models towards learning predictive patterns from irrelevant associations. This dataset has been generated by a directed Gaussian graphical model. The exact procedure is detailed in Suppl. Mat. 1.1. By construction, features can be categorized based on their degree of relevance:

Figure 1: The four datasets, with their predictive features shown by pairs. Orange and blue correspond to the positive and negative classes, respectively. Each column is associated with one dataset, and each row corresponds to a distinct features pair.
* Features \(X_{i}\) that are highly relevant as they are causal for the target variable \(Y\): \(X_{i}\rightarrow\ldots\to Y\)
* Features \(X_{i}\) that are weakly relevant as they correlate with the target variable only due to indirect effects, like in forks: \(X_{i}\leftarrow\ldots\gets X_{j}\rightarrow\ldots\to Y\)
* Irrelevant features

### Benchmark procedure

We assessed the reliability of FS methods on nonlinear ML tasks with a growing degree of difficulty by incrementally diluting the relevant features in the uniformly-distributed datasets. An exponential increase of the number \(k\) of decoy features is added in each run: \(m\in\{2,4\}\cup K\) for XOR and RING datasets, \(m\in\{4\}\cup K\) for RING+XOR and \(m\in\{6\}\cup K\) for RING+XOR+SUM, where \(K=\{8,16,32,64,128,256,512,1024,2048\}\). For each run, we set the number of samples to \(n=1000\).

For each run and dataset, we assessed the performance of both predictors and FS methods. First, for each embedded FS method that relies on a predictive model, we evaluated the latter using a 6-fold cross-validation procedure. Then we computed the AUROC and AUPRC of each step accordingly. We also reported these metrics for a baseline NN without prior or posterior FS. In the paper, we simply refer to it as "Neural Network".

Second, we evaluated the ability of each FS algorithm to rank the predictive features higher than the decoys. More specifically, we quantified the latter as the percentage of predictive features among the highest \(p\) and \(2p\) top-ranked features, where \(p\) corresponds to the number of truly predictive features in each dataset (\(p=2\) for XOR and RING, \(p=4\) for RING+XOR, \(p=6\) for RING+XOR+SUM, \(p=7\) or \(81\) in DAG depending on the definition used). We refer to them as best \(p\) and best \(2p\) scores for short. To ensure fair comparison, and avoid assigning good performance to badly-designed FS methods (where feature importances are influenced by their position/indices in the data matrix), we randomly permuted the columns of the input data matrix in each fold of the \(k\)-fold cross-validation procedure. For feature attribution methods (e.g. saliency maps), only the points from the held-out sets of the 6-fold cross-validation have been used to select features.

### Machine Learning models

To ensure a fair comparison between NN-based FS methods, we tried to reuse the same neural architecture whenever possible. By default, NNs have been implemented with Pytorch [26]. The data has been centered by replacing each input vector \(x_{i}\) by \(2x_{i}-1\), so each feature ranges between \(-1\) and \(1\). The default model is a three-layer perceptron with LeakyReLU activation functions, with a 0.2 slope for negative values. The 2 hidden layers have \(16\) neurons each. Additionally, L2 regularisation on the parameters is used, with a regularisation parameter of \(10^{-2}\). The model was trained with the Adam optimizer [27] for up to \(1000\) epochs, with a learning rate of \(0.005\) and a batch size of 64. A scheduler guided the optimisation of each model by adjusting the learning rate over time, decreasing it by a \(0.9\) factor every time the total loss function stagnates for 10 epochs (with a cooldown of 5 epochs). Gaussian noise, with a standard deviation of 0.05, was added to the inputs in order to regularise the models further. In order to minimize overfitting risks and get the best out of each NN-based approach, we implemented an early stopping criterion. We saved the NN parameters at each epoch and kept the ones that minimize the loss function on a held-out set composed of 20% of the training data. The training was interrupted prematurely if the validation loss had not decreased during the last 5 epochs.

All feature attribution methods are based on the architecture described above. However, due to the peculiarities of some embedded FS methods or the constraints of their implementation, these latter methods might build on subtle variants of this architecture. When relevant, these differences are explained in the next section.

### Feature selection methods

The benchmarked algorithms belong to three categories: feature attribution, embedded, and filter methods. We now describe them in details.

#### 2.4.1 Feature attribution methods

Feature attribution methods propose an _a posteriori_ interpretation of the method \(M\) by approximately reconstructing the decision process followed by \(M\) in order to produce the prediction \(y_{i}\). Gradient-based interpretation methods for NNs like Saliency Maps (SMs) belong to this category. Given a trained NN model \(M\), the forward pass \(M(x_{i})\) of sample \(x_{i}\) is computed, alongside with the gradient \(\partial M(x_{i})/\partial x_{i}\) of the target output \(y_{i}=M(x_{i})\) with respect to the input \(x_{i}\). The gradient values identify which positions in \(x_{i}\) are the most relevant for the prediction. Indeed,
because the gradient points towards the direction of steepest descent, the highest components of the gradient indicate which input variables require the least change to produce the largest variation in the output \(y_{i}\). In this study we used the Captum [16] library to benchmark the Integrated Gradients [11], Saliency [10], DeepLift [12], Input \(\times\) Gradient [13], SmoothGrad [14], Guided Backpropagation [15]. From the same library, we also benchmarked non-gradient-based interpretation methods like Deconvolution [28], Feature Ablation, Feature Permutation [29] and Shapley Value [30] interpretation approaches. For SmoothGrad, data was injected with random noise sampled from a zero-centered Gaussian distribution with 0.1 standard deviation (Captum's implementation of NoiseTunnel). The operation has been repeated \(50\) times per input vector. For the Integrated Gradients, DeepLift, Feature Ablation and Shapley Value Sampling methods, the \(0\) vector was supplied as baseline point. The baseline point is used for different purposes depending on the method. For example, the basepoint provided for the Integrated Gradients method defines the point from which to compute the integral and smooth the feature attribution vector. Because all these gradient-based methods only provide instance-level feature importances, we computed the overall feature importances as the average of absolute values of the instance-level importances. Because it was _a priori_ unclear to us whether these instance-level scores should be computed on the training set or the validation set (with their difference explained by the generalization error), we compared the two settings in section 3.7.

#### 2.4.2 Embedded FS methods

The second category of FS approaches are the embedded methods, which operate _at training time_, such as FSNet [23], Concrete Autoencoder [24], CancelOut [20], Random Forest [31], DeepPINK [21] and LassoNet [22]. In particular, FSNet and Concrete Autoencoder jointly _rank_ features and train the model.

FSNet's architecture is composed of a selector and an encoder network, which branches into a classifier and a final sub-network comprising a decoder and a reconstruction layer. The selector consists in a matrix multiplication operation, involving a low-rank matrix of shape \(k\times 2p\) obtained from a weights predictor, followed by a softmax activation (after addition of Gumbel noise to the logits): \(2p\) features are selected, and we computed both best \(p\) and best \(2p\) scores using the feature importances as proposed in [23]. The weights predictor of the selector is composed of a fully-connected layer with no activation. We chose the identity function as activation for both the encoder and decoder, as no further refinement of the input features was deemed necessary. The reconstruction layer is the reverse operation of the selector network, and performs a matrix multiplication analogously. The corresponding low-rank matrix is predicted from a predictor module composed of a fully-connected layer with no activation. The classifier has the same architecture as described in the previous section, except that its input size is restricted to \(2p\). The whole model was trained for \(2000\) epochs. \(30\) bins (latent size of \(10\)) were used to compute the inputs frequencies necessary to predict the low-rank matrices.

CancelOut is composed of the common architecture described in the previous section (trained in the same manner), preceded by a CancelOut layer. We experimented with 2 variants of the method: 1) with a Sigmoid activation and CancelOut weights regularisation (we denote the corresponding model by CancelOut Sigmoid for short), and 2) with a Softmax activation layer and without regularization (CancelOut Softmax). In the former case, we used a \(\lambda_{1}=0.2\) coefficient for the variance and a \(\lambda_{2}=0.1\) for the regularisation term. Because L1 regularisation encourages importances weights to converge to 0.5 (sigmoid(0) = 0.5), we replaced it by the sum of CancelOut weights, unlike the original implementation. Indeed, this choice of regularisation better promotes sparsity among feature importances. In both cases, CancelOut weights were initialised with the same value \(\beta=1\). The model has been trained for \(300\) epochs, as the convergence of CancelOut weights requires longer time than the optimisation of the classifier alone.

All models built for evaluating the LassoNet [22] feature selection method contained \(32\) hidden neurons and ReLU activation functions, as only the latter were available among activation functions in the lassonet Python package.

To evaluate the Concrete Autoencoder (CAE), we used the

concrete-autoencoder Python package [24], and implemented the underlying classifier with Keras [32] in accordance with the architecture described in the previous section. The CAE has been trained for 300 epochs, with initial temperature \(10\) and final temperature \(0.01\). The CAE has been trained twice, once for selecting \(p\) and once for selecting \(2p\) features. The model has been trained for \(10\) epochs. Each hidden layer is composed of \(32\) neurons, and followed by a 20% dropout.

DeepPINK requires knockoff features, that we generated in two different ways depending on the nature of the dataset. For the DAG dataset, we relied on the framework of Model-X Knockoff features designed by Candes et al. [33]. Details of Model-X Knockoff features are explained in Suppl. Mat. 1.2. For the remaining 4 datasets, knockoff features were generated by simply sampling a uniform distribution (the marginal distribution of the original data matrix \(X\) is multivariate but uniform).
Random forests were grown with the Scikit-learn Python package [34] and are composed of \(500\) trees each. We used the default parameters. We selected features from trained random forests using the widely-used impurity-based feature importance scores, as implemented in Scikit-Learn. More specifically, the importance of a feature is defined as the total reduction in Gini impurity caused by all node splits involving that feature, averaged across all trees in the forest.

#### 2.4.3 Filter methods

The last category of methods that we considered is the filter-method algorithms, such as Relief [35] or minimum redundancy maximum relevance (mRMR) [36]. These methods do not assign predictions; therefore no AUROC or AUPRC values are reported. mRMR relies on statistical measures such as mutual information [37], which requires a discretisation of the input variables. Therefore, we divided each feature into 20 equally-sized bins, which we are sufficiently thin for capturing the nonlinear dependencies, and sufficiently large for robust estimation (\(\sim\) 40 observations per bin). AUROC and AUPRC were computed for FS methods based on a predictive model, and reported as a function of the input feature size.

## 3 Results

Feature Selection (FS) is widely used across many fields of science [38; 39; 40]. The goal of FS algorithms is to identify or rank features in function of the _predictive signal_ they carry with respect to a prediction label. In this study we benchmarked 20 FS methods on 5 synthetic datasets (RING, XOR, RING+XOR, RING+XOR+SUM and DAG) providing nonlinear binary classification tasks (see Methods). These datasets represent classical ML problem, such as the XOR problem, the discrimination of points lying on a ring-shaped subspace and combinations thereof (see Methods for more details).

### Random forests outperform other methods on RING by a large margin

In the RING dataset, positive labels are associated to the points lying on a bi-dimensional ring defined by features in positions \(0,1\) (see Suppl. Fig. S1).

Top panels in Fig. 2 show the AUROC and AUPRC of all trained models, as a function of the total number of features \(m=p+k\). It is noteworthy that all models but the Random Forest have AUROCs and AUPRCs approaching a random predictor (50%) for values of \(m\) greater than 32. This effect is consistent with the percentages of relevant features reported in the bottom graphs of panel Fig. 2. Both Random forests and TreeSHAP perfectly re-identified the relevant features, even when \(m=2048\). However, their decaying performance (as \(m\) increases) suggest that tree-based models may lose their FS capabilities in extremely high-dimensional (\(m\gg 2048\)).

### Neural networks are better suited for solving the XOR problem

As shown in Fig. 3, NN models obtained overall better results, and mRMR underperformed, even for a small number of irrelevant features \(k\). In particular, LassoNet reached maximal best \(p\) and best \(2p\) scores for each \(m\), and \(>90\%\) AUROC and AUPRC for each \(m\leq 512\). This relatively high predictive performance of LassoNet can be attributed to the internal cross-validation procedure used to find the optimal degree of sparsity of the NN parameters. The FS methods, besides LassoNet, that found the highest proportions of relevant features are Relief, CancelOut Sigmoid, and feature attribution methods like DeepLift or Deconvolution. Most methods undergo drastic performance losses when \(m\geq 128\).

### RF, mRMR and LassoNet are the best performing methods on RING+XOR

This dataset comprises four relevant features (\(p=4\)). Similarly to what we observed for the RING dataset, we see from Fig. 4 and Tab. 2 that Random Forests and TreeSHAP (which interprets learned Random Forests) perform the best, along with mRMR. Random Forests achieve high AUROC and AUPRC values and correctly rank the features in most settings. In particular, it constitutes the only FS method capable of ranking over the \(60\%\) of relevant features among the top 4 features for \(m\leq 2048\). For any other method to achieve a similar percentage of relevant features selected \(>80\%\), the total number of features needs to be decreased to \(m=16\), thus requiring to make the problem orders of magnitude simpler. Among these remaining methods, mRMR dominates regardless of the number of features. Overall, only RF, mRMR and LassoNet result in best \(p\) and best \(2p\) scores \(>50\%\) when \(m>256\). All remaining methods except LassoNet appear to perform close to random for \(m\geq 256\). LassoNet starts performing random for \(m\geq 1024\)
Figure 3: Performance of the different models and feature selection methods on the XOR dataset. (Top) AUROC and AUPRC of each trained model as a function of the number of features. (Bottom) Percentage of relevant features selected by each FS method in the top \(p\) and \(2p\), respectively. Shaded areas correspond to the best \(p\) and best \(2p\) scores of a dummy FS method that performs worse than random. Methods have been sorted by decreasing order of average performance in the legend.

Figure 2: Performance of the different models and feature selection methods on the RING dataset. (Top) AUROC and AUPRC of each trained model as a function of the number of features. (Bottom) Percentage of relevant features selected by each FS method in the top \(p\) and \(2p\), respectively. Shaded areas correspond to the best \(p\) and best \(2p\) scores of a dummy FS method that performs worse than random. Methods have been sorted by decreasing order of average performance in the legend.
### Ring+XOR+SUM is challenging for all methods

There are six relevant features in this dataset, two of them being linear with some additive noise. The prediction problem appears to be quite difficult, with AUROC and AUPRC values \(<\) 0.85 even for \(p=6\) when no decoy feature is added (see Fig. 5). Contrary to the previous experiments, overall performance drops gradually with \(m\). This can be explained by the slightly larger number of relevant features (\(p=6\)) as well as their heterogeneity, allowing models to detect the features that are the easiest to pick, given the characteristics of each model (XOR features for NNs, RING features for mRMR and tree-based models, SUM features for most models). RF/TreeSHAP and mRMR are the best FS methods, with the best \(p\) scores \(\geq 40\%\) and the best \(2p\) scores \(\geq 60\%\) for each value of \(m\). Concrete autoencoder ranks last as FS method, and DeepPINK ranks last as a predictive model. Finally, we observe that CancelOut Sigmoid clearly outperforms its softmax variant by a large margin.

### Disentangling causal from spurious effects on DAG is too challenging for FS methods

In the fifth and last dataset in this benchmark, which has been generated by a graphical model (see Methods), relevant features can be defined either as features that are causal for the observed variable \(Y\), or features that are expected to correlate with \(Y\) due to indirect (confounding) effects. We refer to this dataset as DAG. In Table 1, we reported the best \(p\) and \(2p\) scores based on these two definitions. We can see that, in both settings, TreeSHAP and its underlying Random Forest model are outperforming other methods. Concrete autoencoder, FSNet, CancelOut (softmax) and DeepPINK fail at detecting relevant features. Among feature attribution methods, Input \(\times\)Gradient, Feature permutation and Shapley value sampling slightly improve over the other methods. Overall, none of the benchmarked approaches seems to be relatively better at disentangling indirect correlations from causal effects on the DAG dataset.

### Benchmark summary

Complementary to the presented figures, we summarised in Tab. 2 the best \(p\) and best \(2p\) scores for all FS methods on each dataset. For readability purposes, we only report the mean scores, computed across all values of \(m\). The table is organised in two parts: instance-level feature attribution methods (_a posteriori_ gradient-based methods such as Saliency Maps), and the remaining approaches. Best performing methods from each category are highlighted in bold. Among model-based and _a priori_ FS techniques, TreeSHAP and its underlying model Random Forests both

Figure 4: Performance of the different models and feature selection methods on the RING+XOR dataset. (Top) AUROC and AUPRC of each trained model as a function of the number of features. (Bottom) Percentage of relevant features selected by each FS method in the top \(p\) and \(2p\), respectively. Shaded areas correspond to the best \(p\) and best \(2p\) scores of a dummy FS method that performs worse than random. Methods have been sorted by decreasing order of average performance in the legend.
outperform all approaches by a large margin on all datasets but XOR. Overall, mRMR appears to be the second best-performing technique, despite its underperformance on XOR. Complementary to mRMR, LassoNet achieved maximal performance on XOR, while getting average results on the remaining datasets.

Instance-level _a posteriori_ methods do not show significant differences, except on the DAG dataset. In the latter setting, Input \(\times\) Gradient, Feature Ablation and Shapley value sampling seem to perform relatively better.

In order to compare the methods in situations where the curse of dimensionality is exacerbated, we reported the same results on sub-sampled versions of the same datasets, with \(n\in\{250,500\}\). These results are shown in Supp. Tables 1 and 2. We observed that mRMR, tree-based methods and LassoNet consistently outperform other methods in the same settings.

### Bootstrapping influences the quality of feature attribution scores

Since practitioners may use _a posteriori_ SM methods in different ways, and since it was unclear whether using points from the training set can indeed improve the quality of feature ranking, we investigated whether the latter can benefit from bootstrapping and the use of the training set. In the first setting, we computed one feature importance vector per instance from the held-out set after training, and averaged them across the whole held-out set. In a second stage, we added a bootstrapping component by repeating this step 10 times and re-training the model each time on a random sample composed of 80% of the training set (sampling with replacement). The final importance vectors have been obtained by averaging across the 10 runs. In the third setting, we removed bootstrapping but computed the feature importance vectors on the points from the training set only, and left the held-out set unused.

## 4 Discussion

### Feature selection and modeling quality are interdependent

Care should be taken in interpreting the results presented in our study. Indeed, each FS method necessarily relies on some modeling assumptions, and the relevance of selected features is heavily impacted by the model's adequacy to the data. In particular, the inference process does not guarantee to yield the optimal solution (the set of parameters that produces the least generalisation error). Despite the flexibility of NNs, as described by the universal approximation

Figure 5: Performance of the different models and feature selection methods on the RING+XOR+SUM dataset. (Top) AUROC and AUPRC of each trained model as a function of the number of features. (Bottom) Percentage of relevant features selected by each FS method in the top \(p\) and \(2p\), respectively. Shaded areas correspond to the best \(p\) and best \(2p\) scores of a dummy FS method that performs worse than random. Methods have been sorted by decreasing order of average performance in the legend.
Figure 6: Average best \(2p\) score on each of the 5 datasets, using three different approaches to infer feature importances from instance-level feature attribution methods. In the first setting, only the points from the held-out set have been used to rank features. In the second case, bootstrapping was performed by randomly sampling (with replacement) 80% of the training set before training the NN. In the last setting, only the points from the training set have been used.
theorem [41] for example, the optimal architecture choice (with smallest generalisation error) is unknown and can only be found by cross-validation. The eventual lack of regularisation, coupled with the presence of a large number of input features, therefore increasing the effective number of parameters, is likely to drive the model towards learning from irrelevant correlations. In such case, because the model's decision process relies on irrelevant features, the FS method building on this model will necessarily attach higher importance to those features. Therefore, FS techniques based on NNs might be hard to exploit in practice, as they require building the most accurate model, and guiding the inference process to the optimal solution. This goal is easier to reach when provided with sufficiently deep insights about the data. In this sense, FS resembles a chicken-and-egg problem, and because of that FS methods cannot be blindly applied on a dataset without minimal prior knowledge of the data and the limitations of the model used.

In particular, all feature attribution / SM methods considered in this study, as well as many of the other FS techniques, rely on a neural network that requires proper training. Because the number of parameters varied with the number of input features, it remains highly probable that the corresponding models either under- or overfitted the data in some situations, regardless of the presence of dropout and L2 regularisation. Overall, there is no guarantee that the Adam optimizer consistently guided the model to the globally-optimal solution, or that the generalisation error was minimal.

### Tree-based modeling and decision tree induction are two separate concepts

Random forests have largely outperformed other methods on the RING, RING+XOR, RING+XOR+SUM and DAG datasets, both as a FS technique and as predictive models. However, the dataset where RFs perform comparatively (to other methods) worse is XOR, which consists of data points obeying a simple logical rule. Such rule can be perfectly captured through a piece-wise linear function. In particular, among all off-the-shelf ML models, decision trees are the optimal choice for modeling such data, as only three decision splits should be theoretically sufficient for perfect segregation of the classes. However, performance of RFs is sub-optimal as observed in Fig. 3.

\begin{table}
\begin{tabular}{l r r r r r} \hline \hline Dataset & \multicolumn{2}{c}{\% causal features} & \multicolumn{2}{c}{\% correlated features} & \multicolumn{2}{c}{Performance (\%)} \\ \cline{2-7} Method & Best p & Best 2p & Best p & Best 2p & AUROC & AUPRC \\ \hline Saliency maps & 14.3 & 21.4 & 9.5 & 14.6 & - & - \\ Integrated gradient & 14.3 & 21.4 & 9.7 & 14.6 & - & - \\ DeepLift & 14.3 & 21.4 & 9.7 & 14.6 & - & - \\ Input \(\times\) Gradient & 19.0 & 26.2 & 10.1 & 14.2 & - & - \\ SmoothGrad & 14.3 & 21.4 & 9.5 & 14.6 & - & - \\ Guided backpropagation & 14.3 & 21.4 & 9.5 & 14.6 & - & - \\ Deconvolution & 14.3 & 21.4 & 9.5 & 14.6 & - & - \\ Feature ablation & 19.0 & 26.2 & 9.7 & 14.2 & - & - \\ Feature permutation & 16.7 & 21.4 & 10.3 & 14.4 & - & - \\ Shapley value sampling & 16.7 & 26.2 & 10.3 & 14.6 & - & - \\ \hline mRMR & 16.7 & 19.0 & 6.6 & 11.5 & - & - \\ LassoNet & 14.3 & 14.3 & 6.8 & 11.1 & 71.9 & 66.0 \\ Relief & 14.3 & 14.3 & 6.8 & 9.5 & - & - \\ Concrete Autoencoder & 2.4 & 11.9 & 4.9 & 8.0 & 50.1 & 52.8 \\ FSNet & 2.4 & 4.8 & 3.5 & 7.2 & 56.7 & 55.5 \\ CancelOut (softmax) & 2.4 & 4.8 & 5.8 & 8.2 & 46.5 & 48.5 \\ CancelOut (sigmoid) & 16.7 & 21.4 & 9.7 & 13.8 & 56.2 & 56.8 \\ DeepPINK & 0.0 & 0.0 & 4.3 & 8.4 & 50.0 & 50.0 \\ Random Forest & 21.4 & 40.5 & 12.8 & **17.9** & **75.4** & **73.9** \\ TreeSHAP & **28.6** & **42.9** & **13.4** & 17.3 & - & - \\ \hline \hline \end{tabular}
\end{table}
Table 1: Best p and best 2p scores of all feature selection methods on the DAG dataset. \(p=7\) when considering only causal features as relevant (first multi-column) and \(p=81\) when also including confouders (second multi-column). Additionally, AUROC and AUPRC scores have been reported for embedded methods.
The quality of its feature selection and prediction is not linked to the sophistication of the modeling _per se_, _but rather the optimality of the decision tree induction algorithm_. Indeed, Scikit-learn's implementation is based on the heuristic CART algorithm [42], which is unlikely to infer the tree with highest information gain and minimal number of nodes. In particular in the XOR dataset, selecting one of the two relevant features as first decision split is not sufficient, as it does not produce any change in the class proportions within the newly obtained hyper-paraleipipeds. Therefore, RFs are counter-intuitively better at growing from ring-shaped data, since any split on one of the 2 corresponding features results in an increase of class purity (strictly positive information gain). On the XOR dataset, optimal inference would require a lookahead of 1 feature, or bivariate decision splits. In conclusion, the relatively lower performance of RFs on the XOR dataset can mostly be attributed to the sub-optimality of its underlying tree induction algorithm.

### Univariate feature selection remains relevant in a high dimensionality context

Although the datasets have been constructed in a way that they are highly challenging for both linear and univariate FS methods, it must be noted that (nonlinear) univariate filter methods remain highly relevant in some contexts. First, they are computationally more efficient (and embarrassingly parallel), making them competitive in a high-dimensional setting (e.g. Whole Genome Sequencing data). Second, they are still capable of capturing relevant features when they _individually_ correlate with the explained variable, in a nonlinear fashion. This is shown by the maximal performance (\(100\%\) best \(p\) score for any value of \(m\)) of mutual information (MI) on the RING dataset (as shown in Suppl. Tab. 3), suggesting that mutual information (MI) consistently detects the reduction of entropy caused by the ring-shaped function that generated the data. Indeed, this ring-shaped data is leaking information about the class at the level of individual features (the probability distribution of each feature is altered, when conditioned on the class distribution). Because real-life problems are more likely to exhibit univariate nonlinear correlations than our artificial datasets, simple information-theoretic approaches could remain highly relevant. However, data availability is a crucial prerequisite for accurate estimation of MI.

## 5 Conclusion

In this paper, we investigated the usefulness of nonlinear FS approaches in the case of high-dimensional data with a low sample-to-feature ratio. What emerges from the results is that Random Forests and Relief outperform neural network-based approaches on all datasets that exhibit nonlinear correlations between individual features and the target

\begin{table}
\begin{tabular}{l r r r r r r r r r} \hline \hline Dataset & \multicolumn{2}{c}{RING} & \multicolumn{2}{c}{XOR} & \multicolumn{2}{c}{RING+XOR} & \multicolumn{2}{c}{RING+XOR+SUM} & \multicolumn{2}{c}{DAG} \\ \cline{2-10} Method & Best k & Best 2k & Best k & Best 2k & Best k & Best 2k & Best k & Best 2k & Best k & Best 2k \\ \hline Saliency maps & 25.8 & 34.8 & 53.8 & 54.5 & 35.0 & 42.5 & 56.1 & 60.8 & 14.3 & 21.4 \\ Integrated gradient & 24.2 & 36.4 & 53.0 & 54.5 & 35.4 & 41.7 & 55.3 & 60.0 & 14.3 & 21.4 \\ DeepLift & 24.2 & 36.4 & 53.0 & 55.3 & 34.6 & 42.5 & 55.3 & 60.0 & 14.3 & 21.4 \\ Input \(\times\) Gradient & 26.5 & 34.8 & 54.5 & 54.5 & 35.8 & 43.3 & 55.8 & 60.8 & 19.0 & 26.2 \\ SmoothGrad & 25.8 & 36.4 & 53.8 & 54.5 & 36.2 & 42.9 & 56.1 & 61.1 & 14.3 & 21.4 \\ Guided backpropagation & 25.8 & 36.4 & 53.8 & 54.5 & 35.4 & 42.5 & 56.4 & 61.1 & 14.3 & 21.4 \\ Deconvolution & 25.8 & 35.6 & 53.8 & 55.3 & 35.0 & 43.3 & 55.6 & 60.8 & 14.3 & 21.4 \\ Feature ablation & 25.0 & 34.8 & 54.5 & 54.5 & 34.2 & 42.1 & 56.1 & 60.6 & 19.0 & 26.2 \\ Feature permutation & 25.8 & 35.6 & 53.0 & 54.5 & 32.9 & 41.7 & 55.6 & 60.0 & 16.7 & 21.4 \\ Shapley value sampling & 23.5 & 37.9 & 52.3 & 53.8 & 35.4 & 42.9 & 55.3 & 60.6 & 16.7 & 26.2 \\ \hline mRMR & **100.0** & **100.0** & 12.5 & 28.3 & 81.7 & 88.8 & 74.7 & 84.7 & 16.7 & 19.0 \\ LassoNet & 34.8 & 35.6 & **81.8** & **81.8** & 44.6 & 52.9 & 64.2 & 67.8 & 14.3 & 14.3 \\ Relief & 40.2 & 45.5 & 72.7 & 74.2 & 37.1 & 42.1 & 43.9 & 53.3 & 14.3 & 14.3 \\ Concrete Autoencoder & 19.7 & 24.2 & 22.7 & 27.3 & 19.2 & 30.0 & 25.8 & 36.4 & 2.4 & 11.9 \\ FSNet & 21.2 & 28.0 & 31.1 & 38.6 & 25.4 & 33.3 & 33.3 & 43.1 & 2.4 & 4.8 \\ CancelOut (softmax) & 14.4 & 25.0 & 22.0 & 28.8 & 18.8 & 29.2 & 30.0 & 38.1 & 2.4 & 4.8 \\ CancelOut (sigmoid) & 21.2 & 31.1 & 67.4 & 72.0 & 35.8 & 46.7 & 58.6 & 63.6 & 16.7 & 21.4 \\ DeepPINK & 17.4 & 25.0 & 17.4 & 24.2 & 20.4 & 32.9 & 28.1 & 37.8 & 0.0 & 0.0 \\ Random Forest & **100.0** & **100.0** & 54.5 & 59.8 & **88.8** & **95.4** & **85.3** & **90.8** & 21.4 & 40.5 \\ TreeSHAP & **100.0** & **100.0** & 40.2 & 43.9 & 81.7 & 88.3 & 78.3 & 86.1 & **28.6** & **42.9** \\ \hline \hline \end{tabular}
\end{table}
Table 2: Best k and best 2k score percentages on the 5 datasets. For the first 4 datasets, scores have been averaged over \(m\in\{2,4,6,8,16,32,64,128,256,512,1024,2048\}\).Top and bottom parts of the table correspond to instance-level feature attribution and embedded/filter FS methods, respectively. Best performing methods are highlighted in bold.
variable (all but the XOR dataset). Therefore, DNN models might not be the best choice for feature selection on datasets with these characteristics, i.e. a low density of relevant features, a prevalence of nonlinear patterns and variance homogeneity (for both predictive and decoy features). In real-life applications, we mainly encounter that relationships among the features in the data tend to be additive. However, our study indicates that when both additive and nonlinearly entangled features are present, we can only see the former unless a considerable number of samples are available compared to the input feature size.
