# How good Neural Networks interpretation methods really are? A quantitative
benchmark.

Antoine Passemiers1,* Pietro Folco2,* Daniele Raimondi1,† Giovanni Birolo2
Yves Moreau1 Piero Fariselli2

###### 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 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.

††footnotetext: 1ESAT-STADIUS, KU Leuven, Leuven, Belgium††footnotetext:
2Department of Medical Sciences, University of Torino, Torino,
Italy††footnotetext: *These authors contributed equally to this
work.††footnotetext: †Corresponding author<EMAIL_ADDRESS>
##  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

###  2.1 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:

$\lvert\sqrt{(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}}-0.35\rvert\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:

$\displaystyle\lvert\sqrt{(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}}-0.35\rvert\leq
0.0704$ $\displaystyle(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:

$\displaystyle\lvert\sqrt{(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}}-0.35\rvert\leq
0.0479$ $\displaystyle(x_{2}-0.5)(0.5-x_{3})\geq 0.0598$ $\displaystyle
x_{4}+x_{5}+\mathcal{N}(\mu=0,\sigma=0.2)\geq 1.4074$

The four outcomes are shown in Fig.1.

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.

####  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:

  * •

Features $X_{i}$ that are highly relevant as they are causal for the target
variable $Y$: $X_{i}\rightarrow\dotsc\rightarrow 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\dotsc\leftarrow X_{j}\rightarrow\dotsc\rightarrow Y$

  * •

Irrelevant features

###  2.2 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.

###  2.3 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.

###  2.4 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 Candès 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).

###  3.1 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$).

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.

###  3.2 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$.

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.

###  3.3 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 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.

###  3.4 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.

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.

###  3.5 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.

Dataset | % causal features | % correlated features | Performance (%)
---|---|---|---
Method | Best p | Best 2p | Best p | Best 2p | AUROC | AUPRC
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 | - | -
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 | - | -

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.

###  3.6 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 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.

Dataset | RING | XOR | RING+XOR | RING+XOR+SUM | DAG
---|---|---|---|---|---
Method | Best k | Best 2k | Best k | Best 2k | Best k | Best 2k | Best k | Best 2k | Best k | Best 2k
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
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

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.

###  3.7 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.

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.

##  4 Discussion

###  4.1 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 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.

###  4.2 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.

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-paralelipipeds. 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.

###  4.3 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 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.

## References

  * [1] Sebastian Lapuschkin, Stephan Wäldchen, Alexander Binder, Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller.  Unmasking clever hans predictors and assessing what machines really learn.  Nature communications, 10(1):1–8, 2019.
  * [2] Andreas Holzinger, Peter Kieseberg, Edgar Weippl, and A Min Tjoa.  Current advances, trends and challenges of machine learning and knowledge extraction: from machine learning to explainable ai.  In International Cross-Domain Conference for Machine Learning and Knowledge Extraction, pages 1–8. Springer, 2018.
  * [3] Seung-Jean Kim, Kwangmoo Koh, Michael Lustig, Stephen Boyd, and Dimitry Gorinevsky.  An interior-point method for large-scale $\ell_{1}$-regularized least squares.  IEEE journal of selected topics in signal processing, 1(4):606–617, 2007.
  * [4] Robert Tibshirani, Trevor Hastie, Balasubramanian Narasimhan, and Gilbert Chu.  Diagnosis of multiple cancer types by shrunken centroids of gene expression.  Proceedings of the National Academy of Sciences, 99(10):6567–6572, 2002.
  * [5] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learning: data mining, inference and prediction.  Springer, 2 edition, 2009.
  * [6] Daniele Raimondi, Massimiliano Corso, Piero Fariselli, and Yves Moreau.  From genotype to phenotype in arabidopsis thaliana: in-silico genome interpretation predicts 288 phenotypes from sequencing data.  Nucleic acids research, 50(3):e16–e16, 2022.
  * [7] John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Žídek, Anna Potapenko, et al.  Highly accurate protein structure prediction with alphafold.  Nature, 596(7873):583–589, 2021.
  * [8] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.  Human-level control through deep reinforcement learning.  nature, 518(7540):529–533, 2015.
  * [9] Alex Graves, Greg Wayne, Malcolm Reynolds, Tim Harley, Ivo Danihelka, Agnieszka Grabska-Barwińska, Sergio Gómez Colmenarejo, Edward Grefenstette, Tiago Ramalho, John Agapiou, et al.  Hybrid computing using a neural network with dynamic external memory.  Nature, 538(7626):471–476, 2016.
  * [10] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  arXiv preprint arXiv:1312.6034, 2013.
  * [11] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  CoRR, abs/1703.01365, 2017.
  * [12] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences.  CoRR, abs/1704.02685, 2017.
  * [13] Avanti Shrikumar, Peyton Greenside, Anna Shcherbina, and Anshul Kundaje.  Not just a black box: Learning important features through propagating activation differences.  CoRR, abs/1605.01713, 2016.
  * [14] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda B. Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  CoRR, abs/1706.03825, 2017.
  * [15] Jost Tobias Springenberg, Alexey Dosovitskiy, Thomas Brox, and Martin Riedmiller.  Striving for simplicity: The all convolutional net.  arXiv preprint arXiv:1412.6806, 2014.
  * [16] Narine Kokhlikyan, Vivek Miglani, Miguel Martin, Edward Wang, Jonathan Reynolds, Alexander Melnikov, Natalia Lunova, and Orion Reblitz-Richardson.  Pytorch captum.  https://github.com/pytorch/captum, 2019.
  * [17] Julius Adebayo, Justin Gilmer, Michael Muelly, Ian Goodfellow, Moritz Hardt, and Been Kim.  Sanity checks for saliency maps.  Advances in neural information processing systems, 31, 2018.
  * [18] Ahmed Alqaraawi, Martin Schuessler, Philipp Weiß, Enrico Costanza, and Nadia Berthouze.  Evaluating saliency map explanations for convolutional neural networks: a user study.  In Proceedings of the 25th International Conference on Intelligent User Interfaces, pages 275–285, 2020.
  * [19] Weili Nie, Yang Zhang, and Ankit Patel.  A theoretical explanation for perplexing behaviors of backpropagation-based visualizations.  In International Conference on Machine Learning, pages 3809–3818. PMLR, 2018.
  * [20] Vadim Borisov, Johannes Haug, and Gjergji Kasneci.  Cancelout: A layer for feature selection in deep neural networks.  In International Conference on Artificial Neural Networks, pages 72–83. Springer, 2019.
  * [21] Yang Young Lu, Yingying Fan, Jinchi Lv, and William Stafford Noble.  Deeppink: reproducible feature selection in deep neural networks.  arXiv preprint arXiv:1809.01185, 2018.
  * [22] Ismael Lemhadri, Feng Ruan, Louis Abraham, and Robert Tibshirani.  Lassonet: A neural network with feature sparsity.  Journal of Machine Learning Research, 22(127):1–29, 2021.
  * [23] Dinesh Singh, Héctor Climente-González, Mathis Petrovich, Eiryo Kawakami, and Makoto Yamada.  Fsnet: Feature selection network on high-dimensional biological data.  arXiv preprint arXiv:2001.08322, 2020.
  * [24] Abubakar Abid, Muhammad Fatih Balin, and James Zou.  Concrete autoencoders for differentiable feature selection and reconstruction.  arXiv preprint arXiv:1901.09346, 2019.
  * [25] Adriana Romero, Pierre Luc Carrier, Akram Erraqabi, Tristan Sylvain, Alex Auvolat, Etienne Dejoie, Marc-André Legault, Marie-Pierre Dube, Julie G. Hussin, and Yoshua Bengio.  Diet networks: Thin parameters for fat genomics.  In International Conference on Learning Representations, 2017.
  * [26] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [27] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [28] Aravindh Mahendran and Andrea Vedaldi.  Salient deconvolutional networks.  In European Conference on Computer Vision, pages 120–135. Springer, 2016.
  * [29] Christoph Molnar.  Interpretable machine learning.  Lulu. com, 2020.
  * [30] Javier Castro, Daniel Gómez, and Juan Tejada.  Polynomial calculation of the shapley value based on sampling.  Computers & Operations Research, 36(5):1726–1730, 2009.
  * [31] Leo Breiman.  Random forests.  Machine learning, 45(1):5–32, 2001.
  * [32] Francois Chollet et al.  Keras, 2015.
  * [33] Emmanuel Candès, Yingying Fan, Lucas Janson, and Jinchi Lv.  Panning for gold: ‘model-x’ knockoffs for high dimensional controlled variable selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3).
  * [34] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [35] S Francisca Rosario and K Thangadurai.  Relief: feature selection approach.  International journal of innovative research and development, 4(11), 2015.
  * [36] Hanchuan Peng, Fuhui Long, and Chris Ding.  Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy.  IEEE Transactions on pattern analysis and machine intelligence, 27(8):1226–1238, 2005.
  * [37] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger.  Estimating mutual information.  Physical review E, 69(6):066138, 2004.
  * [38] Giovanna Nicora, Francesca Vitali, Arianna Dagliati, Nophar Geifman, and Riccardo Bellazzi.  Integrated multi-omics analyses in oncology: a review of machine learning methods and tools.  Frontiers in oncology, 10:1030, 2020.
  * [39] Parminder S Reel, Smarti Reel, Ewan Pearson, Emanuele Trucco, and Emily Jefferson.  Using machine learning approaches for multi-omics data analysis: A review.  Biotechnology Advances, 49:107739, 2021.
  * [40] Daniele Raimondi, Gabriele Orlando, Wim F Vranken, and Yves Moreau.  Exploring the limitations of biophysical propensity scales coupled with machine learning for protein sequence analysis.  Scientific reports, 9(1):1–11, 2019.
  * [41] George Cybenko.  Approximation by superpositions of a sigmoidal function.  Mathematics of control, signals and systems, 2(4):303–314, 1989\.
  * [42] Leo Breiman, Jerome H Friedman, Richard A Olshen, and Charles J Stone.  Classification and regression trees.  Routledge, 2017.