# A novel evaluation methodology for supervised Feature Ranking algorithms

Jeroen G. S. Overschie
Faculty of Science and Engineering
University of Groningen
Groningen, 9747 AG
<EMAIL_ADDRESS>
Acknowledgements go to George Azzopardi and Ahmad Alsahaf.

(July 1, 2021)

###### Abstract

Both in the domains of Feature Selection and Interpretable AI, there exists a
desire to ‘rank’ features based on their importance. Such feature importance
rankings can then be used to either: (1) reduce the dataset size or (2)
interpret the Machine Learning model. In the literature, however, such Feature
Rankers are not evaluated in a systematic, consistent way. Many papers have a
different way of arguing which feature importance ranker works best. This
paper fills this gap, by proposing a new evaluation methodology. By making use
of synthetic datasets, feature importance scores can be known beforehand,
allowing more systematic evaluation. To facilitate large-scale experimentation
using the new methodology, a benchmarking framework was built in Python,
called fseval. The framework allows running experiments in parallel and
distributed over machines on HPC systems. By integrating with an online
platform called Weights and Biases, charts can be interactively explored on a
live dashboard. The software was released as open-source software, and is
published as a package on the PyPi platform. The research concludes by
exploring one such large-scale experiment, to find the strengths and
weaknesses of the participating algorithms, on many fronts.

_Keywords_ Feature Ranking, Feature Selection, Evaluation Metrics, Benchmark,
Machine Learning

###### Contents

  1. 1 Introduction
  2. 2 Motivations
    1. 2.1 Motivations for Feature Ranking
      1. 2.1.1 Feature Selection
      2. 2.1.2 Interpretable Artificial Intelligence
    2. 2.2 Motivations for Evaluating Feature Ranking methods
  3. 3 Related work
    1. 3.1 Evaluation metrics
    2. 3.2 Synthetic datasets and apriori information
    3. 3.3 The gap in the current literature
  4. 4 Methods for Feature Ranking
    1. 4.1 Terminology
    2. 4.2 Types of Feature Rankings
      1. 4.2.1 Feature importance
      2. 4.2.2 Feature support
      3. 4.2.3 Feature rankings
    3. 4.3 Feature Selection Taxonomy
      1. 4.3.1 Filter methods
      2. 4.3.2 Wrapper methods
      3. 4.3.3 Embedded methods
      4. 4.3.4 Hybrid methods
    4. 4.4 Types of features
    5. 4.5 Constructing Feature Rankings
      1. 4.5.1 Regularized Linear Regression
      2. 4.5.2 Relief-Based Feature Selection algorithms
  5. 5 Evaluating Feature Rankings
    1. 5.1 Cross-Validation and Bootstrapping
      1. 5.1.1 Cross-Validation
      2. 5.1.2 Bootstrapping
    2. 5.2 Validation estimators
    3. 5.3 Apriori knowledge on relevant features
    4. 5.4 Stability
      1. 5.4.1 Stability of feature importance vectors
      2. 5.4.2 Stability of feature support vectors
      3. 5.4.3 Stability of feature ranking vectors
    5. 5.5 Time complexity
    6. 5.6 Statistical integrity
  6. 6 Building a benchmarking pipeline
    1. 6.1 Architecture
      1. 6.1.1 The main script
      2. 6.1.2 The ‘Rank and Validate’ pipeline
    2. 6.2 Components
      1. 6.2.1 Datasets
      2. 6.2.2 Cross-Validation
      3. 6.2.3 Resampling
      4. 6.2.4 Estimators
      5. 6.2.5 Storage providers
      6. 6.2.6 Callbacks
    3. 6.3 Execution
      1. 6.3.1 Multiprocessing
      2. 6.3.2 Distributed computing
    4. 6.4 Data- collection and visualization
  7. 7 Experiments
    1. 7.1 Experiment setup
    2. 7.2 Experimental results for the ‘Synclf hard’ dataset
      1. 7.2.1 ReliefF performance
      2. 7.2.2 Performance of multiple rankers
    3. 7.3 Experimental result for all datasets
    4. 7.4 Learning curves and time complexity
    5. 7.5 Online dashboard
  8. 8 Discussion
  9. 9 Conclusion
    1. 9.1 Concluding note
    2. 9.2 Limitations and Future work
  10. A Experiment line-up
  11. B Additional results

##  1 Introduction

In this day and age, more data is available than ever before in many domains
(Sagiroglu and Sinanc, 2013). In the biomedical domain sensory devices such as
MRI or PET scanners are getting ever more accurate - requiring more storage
space to store higher resolution data. In the financial domain, markets are
operating at increasingly low time intervals - requiring storage and analysis
of data at a higher time resolution than before. The internet too, sees
increasing amounts of traffic world-wide and is producing immense amounts of
data every day. Applications of Machine Learning are common in all these
domains: training predictive models by learning from example data is able to
give us interesting insights that can improve both economy and the quality of
human lives. By the nature of models that learn from examples, performance is
often better when a larger amount of examples is available. But it shall get
clear that larger quantities of data presents itself not as solely beneficial;
but rather- a mixed blessing.

The field of Machine Learning has seen vast increases in dataset sizes: both
in terms of sample size and amount of dimensions. Although the availability of
more data presents practitioners with opportunities to create better
performing models, more data will have to be processed - causing an increased
computational burden in the learning process. This increase is nonlinear: due
to the curse of dimensionality, the computational burden can get large,
quickly. Even though the field has for long relied on computer processing
speed steadily increasing, obeying Moore’s law, the rate of advancement will
inevitably start declining - and in fact already has (Theis and Wong, 2017).
Self-evidently, many technological advancements can still be realised, either
in the silicon world or in a post-silicon world, in which perhaps forms of
quantum computing might become predominant (Britt and Humble, 2017). But what
is certain, is that besides the technological opportunities the chip-making
industry still has, there also exist strict physical limitations as to how
fast computer processing can get. So, besides leveraging faster hardware, we
are also going to have to make our software smarter. Instead of increasing
just our computational power, methods for reducing the computational burden in
the first place are desired. Thus, the need for preprocessing techniques and
data reduction algorithms is instigated.

Feature selection is such a domain that focuses on reducing the overall
learning processing burden (Guyon and Elisseeff, 2003). By figuring out what
dimensions are relevant to the learning task, which ones are redundant, and
which ones are completely noisy with respect to the learning task, a smaller
dataset can be obtained by means of a preprocessing step. This is contrary to
the domain of dimensionality reduction, in which data is projected onto a
space of smaller dimensionality - but losing the exact representation of the
data distribution. In Feature Selection, we aim to obtain a binary decision
about which dimensions to keep, in the original data space. Aside from
reducing dataset size by cutting off dimensions, in some cases the
generalization ability of a learning model can even be improved: the learning
model can better learn the data distribution by using less but more meaningful
dimensions with less noise. This makes the benefit of learning on a subset of
the available features two-fold: model fitting and prediction can be both
faster and more accurate. To select relevant dataset features, a wide variety
of strategies exist. One is to assign a scoring to each dimension and keep
only the most relevant ones - such a ranking is called a feature ranking.

Feature ranking is a broader domain in contrast to feature selection, in which
the sole purpose is not to only reduce the dataset size, but to construct a
hierarchical order on the importance of features given a specific learning
task (Duch et al., 2004). Many techniques can be employed to create such
feature rankings, of which some are substantially faster than others, but
might yield sub-optimal results: the choice of a suitable feature ranking
algorithm is not trivial in most scenarios. Once such a feature ranking has
been constructed, it can be used for various applications. First, a feature
selection can be made by removing features that rank below a certain feature
importance score threshold; one such naive method would be to cut off any
feature that ranks below the average feature importance score. Secondly, the
feature ranking can also be used for better Machine Learning model
interpretation; in which the feature importance scores help humans better
understand the predictions and reasoning of the models - by knowing which
features the model found important, it can better be understood how the model
made the decision that it has. This second application is part of the bigger
domain of Interpretable Artificial Intelligence, or more commonly
Interpretable Artificial Intelligence (AI), which has in recent times become
ever more relevant (Ghosh and Kandasamy, 2020). Interpretable AI aims to
explain models that were before considered ‘black box’ models, making them
more transparent to the user.

Evaluation on the performance of feature ranking algorithms has been conducted
in many different ways. Most authors used a ‘validation’ estimator, which was
trained on a chosen subset of the dataset to then see how the estimator
performed given this feature subset. A ranker is desired, then, that ranks the
most predictive features highest, preferably in a reasonable amount of time.
In this way, we get a feature subset that is as small as possible, that gives
us the highest possible predictive power. This evaluation technique, however,
might not be systematic enough. Across papers, many different validation
estimators are used, which make the results across papers subsequently
incomparable to one another. Researchers might also benefit from more
extensive and systematic evaluation by use of synthetic datasets. By manually
controlling many aspects of the dataset, such as noise levels, the complexity
of the data distribution to be learned and the amount of informative features,
a comprehensive evaluation on the feature ranking algorithm behavior and
characteristics can be made. In this way, by employing synthetic datasets, the
exact informative features to be ranked as relevant can be known apriori, i.e.
before conducting the feature ranking operation. Such new and possibly useful
metrics can be employed to evaluate feature rankings independent of any
validation estimator.

In this paper, a comprehensive comparative experiment on feature rankers is
conducted using both real-world and synthetic datasets, employing new
evaluation metrics on the synthetic datasets by knowing the relevant feature
apriori. Both classical and more recent feature ranking algorithms are
included, using methods that originally reside in both the statistical and
feature selection domains. By systematically generating synthetic datasets
that are specifically designed to vary in various relevant data properties,
various characteristics of the feature ranking algorithms can be assessed and
estimated. To employ such a large-scale benchmark, a software framework was
built to facilitate such testing. The framework was released as open-source,
freely available software.

The research question which is to be answered is as follows. “How do we
evaluate Feature Ranking algorithms in a way that is systematic,
comprehensive, and emphasizes the differences between the algorithms in a
meaningful way?” Currently, Feature Ranking and Feature Selection algorithms
are evaluated in many ways. This research aims to find a methodology that is
both meaningful and applicable to many algorithms. The goal is to arrive at a
scientifically sound benchmarking method for feature ranking methods, taking
advantage of apriori information wherever possible. Up to our knowledge, no
current literature exists that is aimed at the evaluation of Feature- Ranking
and Selection algorithms, that takes such a comprehensive approach and
proposes novel evaluation metrics. An important sub-goal in constructing such
a benchmarking method is the development of a pipeline implementing such a new
benchmarking method and analysing its results.

The scope of this research can be defined as follows. Feature Ranking methods
are evaluated that work on tabular datasets and require example data including
prediction targets, i.e. only supervised algorithms are considered.
Furthermore, all considered datasets are tabular, that is, no underlying data
structures such as linked- or streaming data is assumed. The considered
dataset prediction tasks are regression and classification \- which limits the
considered ranking methods to these tasks as well. Eight out of fourteen
classification datasets perfectly balanced classes, the other six have varying
levels of imbalance (see Appendix). All considered Feature Ranking methods are
global rankers. Meaning that all methods construct rankings for the full
dataset, i.e., no instance-based methods are included. Lastly, all three
ranking types are supported and are considered in the research, i.e. the
research considers feature importance-, feature support- and feature ranking
vectors.

The contribution of this paper is multiple-fold. Firstly, the inclusion of
many feature ranking algorithms and many datasets makes it possible to make
more meaningful comparisons between feature ranking algorithms, which would
not have been possible across existing papers due to the lack of a single
evaluation standard. Secondly, the proposal of a new evaluation standard also
makes it possible for other authors to conduct experiments in reproducible
manner. This subsequently enables readers to compare results across papers -
saving time but also allowing more thorough comparative analysis on which
feature ranker is best suited for a given dataset. Third, the new evaluation
standard was implemented and packaged in an easy-to-use open-source software
package, distributed as a Python pip package manager distribution on the PyPi
platform111https://pypi.org/project/fseval/.

Chapters in this paper are structured as follows. First, motivations for both
Feature Ranking and the evaluation thereof are given in Chapter 2. Second, an
analysis on previous work in the literature is conducted, in Chapter 3. Third,
various methods for creating feature rankings are presented, such that a grasp
can be obtained on the overall mechanics of the methods. This is done is
Chapter 4. Fourth, an insight into how feature rankings are evaluated is
gained; and a new method is presented afterwards, in Chapter 5. Next, comments
are made on the construction of a benchmarking pipeline in Chapter 6. Then,
this new standard is applied in an experiment, which setup and results are
explained in Chapter 7. Lastly, the paper is concluded by a discussion in
Chapter 8 and a conclusion in Chapter 9.

##  2 Motivations

###  2.1 Motivations for Feature Ranking

The desire to create a feature ranking from some features in a dataset stems
from two main applications: the domain of feature selection and the domain of
interpretable AI. To better understand why one would want to create a feature
ranking, a brief exploration is made on both domains, motivating the concept
throughout.

####  2.1.1 Feature Selection

The field of Machine Learning is plagued by the curse of dimensionality
(Köppen, 2009): as the dimensionality of datasets grows larger, the
computational burden for learning algorithms gets exponentially larger. This
is due to the fact that when data gets of increasingly higher dimensionality,
the volume of the data shape grows faster than the amount of samples grows
along, causing the data volume to become sparse. There are many real-world
domains that naturally deal with datasets of such shape. In the biomedical
world, for example, collecting data samples might require arduous amounts of
human effort (Hu et al., 2018). There might be the case where one data sample
represents one human and any such sample is very laborious to collect: but the
samples that are collected are of very high dimensionality. In this case, we
deal with a scenario where the amount of dimensions far exceeds the amount of
samples ($n\gg p$ where $n$ is the amount of samples and $p$ the amount of
dimensions).

To avoid such a curse of dimensionality, feature selection can be used. The
goal is to reduce as much features as possible while retaining the most
relevant information in the dataset. The benefits are a decreased
computational burden; though in some cases the learning algorithm
generalization performance might actually be improved over a scenario when no
feature selection is used. This is most often due to the removal of noise that
would have interfered with the learning process.

One such example process of feature selection by using a feature ranking can
be seen in Figure 1. As can be seen in the figure, a ranking can be used to
reduce the overall dataset size by using a threshold operation. In such a
threshold operation, any feature that ranks below the given threshold
$\epsilon$ is removed. Many, if not most, feature selection methods rely on
feature ranking methods under the hood to create a feature subset of reduced
size.

Figure 1: An example process of feature selection by thresholding a feature
ranking at some point $\epsilon$. The resulting feature subset $\mathbf{X}$
presumably contains the strongest predictors given the target variable.
Feature names come from a real dataset on antarctic penguins (Horst et al.,
2020).

For this reason, the more general problem of constructing feature rankings is
omni-relevant in the feature selection domain. To go from feature ranking to
feature selection, the only extra step necessary is a sufficient threshold
point, which many algorithms arrive at using different methods. A
straightforward way to cut off low-ranking features is to exclude any feature
that ranks below the mean score - or, alternatively, exclude only the highest
ranking quartile, etcetera. More sophisticated schemes exist, however. One
might also choose to iteratively run a feature ranking process and remove the
lowest ranking score in each iteration. Such a process is called Recursive
Feature Elimination (RFE). That said, it is clear that feature ranking and
feature selection go in unison.

####  2.1.2 Interpretable Artificial Intelligence

The domain of Interpretable AI has arisen due to the need to better understand
and reason about the decisions that any Machine Learning model makes.
Traditionally, the sole goal in building a learning algorithm would be to best
predict some target variable, given some set of samples to learn from. With
the advent of sophisticated learning algorithms such as Neural Networks and
especially Deep Neural Networks, however, the many computational layers that
separate inputs from answers tend to obfuscate the decision process (Rai,
2020). Whereas in classical statistical models practitioners sought to unveil
a clear relationship between the independent- and dependent variables, some
models in the field of Machine Learning have grown so complicated that no
human can reason on its output. Because such algorithms tend to compare
metaphorically to a black box one cannot possibly see through, one also refers
to such algorithms as black boxes.

Whilst at the same time models have gotten progressively more sophisticated
and thus less transparent, the market for deploying Machine Learning has been
steadily getting bigger. Many institutions seek to benefit from the
possibilities of recent developments in AI and Machine Learning (ML)\-
including many companies, schools or the government. In all such applications,
for every decision made by a Machine Learning model, an argumentation as how
the model came to such a conclusion is desired. Even, in some scenarios the
decision to be made at hand weighs so heavily that a ML model that cannot
explain itself might not be usable at all. For example, one might deploy a ML
model that scores employees based on their performance. The goal is then to
fire the lowest ranking employees and keep highest performing ones. But fully
trusting such a model might be a dangerous practice. If a model provides only
bare explanations on how it came to such a decision, employers might have a
hard time interpreting the model and employees might have a hard time
understanding its decisions. As a matter of fact, faulty models of such kind
have already been deployed out in the open - with considerable complications
with respect to model transparency as a result (O’Neil, 2016).

For this reason, decisions made by computers are required to be explainable -
allowing practitioners to assess not only ML model decisions, but also assess
the overall usefulness of the model in general. After all, many learning
models are designed such to always generate a decision, no matter how noisy
the input. That said, the black box of AI can generally be opened up in two
ways. First, one could opt for a less sophisticated in the first place - one
that reveals its inner workings and supports its decisions by an elaborate
scheme, much in line with classical statistics. Although a simpler model could
sometimes suffice in places where practitioners currently opt for more
complicated ones, a solution for models of any kind of flexibility is desired.
Therefore, a second option is to delve into any such model to reveal a
reasoning about its decision-making process - which is the option that the
field of Interpretable AI bothers itself with.

Feature rankings are one of the facets which can be used to facilitate a
better understanding of a Machine Learning model (Hansen and Rieger, 2019): by
unveiling which variables the ML model finds important, a better understanding
can be gained on its decision-making process. In some scenarios, it might, for
example, be the case that the Machine Learning model in question unexpectedly
weighs a feature as very important, even though a human expert can know
apriori that the feature is not of value to the predictive task at hand. In
this way, faulty models can be detected and prevented from being used in
critical settings. In the Interpretable AI jargon, the process of scoring and
explaining feature importance goes under different names. In the literature,
feature rankings are referred to as feature impact scores, feature importance
scores, or feature effects, which are synonymous in this case. All terms
indicate the process of quantifying feature relevance with respect to the
learning task at hand. The general process of an interpretable AI algorithm
can be seen as in Figure 2.

Figure 2: An example process of Interpretable AI. In this case, the input to
the interpretable AI algorithm is a trained black box model, which has an
otherwise hard to understand decision making process. The interpretable AI
algorithm then extracts feature importance scores for better interpretability.

As could be seen in Figure 2, one way of making a black box model more
interpretable is to extract a feature ranking from the trained model - an
interpretable AI algorithm reads in the black box model and tries to explain
it (Lundberg and Lee, 2017). Other interpretable AI methods take a more direct
approach (Arik and Pfister, 2020), in which the interpretable AI algorithm
functions as a prediction method itself. The latter approach is employed only
in more recent years; in which learning algorithms are designed from the start
to be interpretable. However, independent of whether a direct- or indirect
approach to explaining the decision process is taken, the desire to rank
features is applicable to both. That said, we can therefore conclude that the
construction of truthful feature rankings are an important aspect in
interpretable AI.

###  2.2 Motivations for Evaluating Feature Ranking methods

Given the need for Feature Ranking methods, the desire to evaluate such
methods systematically is clear. Practitioners throughout the field of Machine
Learning desire to know the behavior of different methods under a large range
of varying conditions, such as the sample size-, dimensionality-, or domain of
the dataset. One also needs to take into account the available computational
resources and the desired degree of accuracy of the solution. Finding out what
evaluation methods and metrics best represent the optimality of the solution
is essential, for else the presentation of scientific findings might be
misleading compared to real-world performance.

In this way, individuals and institutions are motivated to actively research
Feature Ranking methods and evaluate them such to determine which method works
best, under what conditions. As is more often the case in the comparison of
Machine Learning algorithms, however, is the fact that methods are not always
easily analytically compared. That is, solely a theoretical analysis does not
suffice to create a meaningful comparison. Although one can reason about the
performance of many methods before conducting concrete data analyses, a real-
world experiment is essential to any reputable publication.

##  3 Related work

The domain of feature ranking and selection has a large availability of
literature, spread out over many subtopics. What is more rare, however, is to
find papers that explicitly research and reason about the usage of certain
evaluation metrics. In general, papers tend to stick to a certain evaluation
method when the majority employs the given technique - but a chance for
conducting a more thorough analysis might be missed nonetheless.

In the feature selection domain, the evaluation and comparison of feature
selection algorithms is a non-trivial problem. Among a wide range of metrics,
no consensus exists among researchers, leaving many papers to present outcomes
in different ways (Guyon and Elisseeff, 2003). In absence of a single
consistent evaluation pipeline across the field, many scholars adhere to
methods that are ‘widely used’ (Solorio-Fernández et al., 2020) (Li et al.,
2017).

Recommendations for metrics have been given in previous papers, most often
when discussing future work. Arguments are made for relevant aspects to
evaluate, such as in Chandrashekar (2014) (Chandrashekar and Sahin, 2014):

> "a feature selection algorithm can be selected based on the following
> considerations: simplicity, stability, number of reduced features,
> classification accuracy, storage and computational requirements"

Of these aspects most proposals focus mainly on number of reduced features,
classification accuracy and computational requirements. In the regression
case, the classification accuracy would be replaced by a commonly used
regression counterpart, the R2-score. Let us explore what validation
estimators and corresponding metrics are used in papers across the field.
Afterwards, evaluation aspects are covered that are not present in the
aforementioned set.

###  3.1 Evaluation metrics

A validation estimator is often used to evaluate supervised feature selection
methods; assessing the quality of a feature subset by running some predictor
over the feature subset selected by the feature selection algorithm, obtaining
the easily interpretable prediction accuracy metric in the classification
case. Predictors often used in the literature include k- Nearest Neighbors
(k-NN) (Al-Tashi et al., 2020) (Mafarja et al., 2020), Support Vector Machines
(SVMs) (Chandrashekar and Sahin, 2014), Decision Trees (DTs) (Li et al., 2017)
and Naïve Bayes (NB) (Koller and Sahami, 1996). Metrics used are often
classification accuracy or in some cases average error rate (Khurma et al.,
2020), validated using some $n$-fold cross validation, commonly 5- or 10-fold.
See Table 1.

Name method |  Validation estimators |  Acc- uracy |  Stab- ility |  Time $t(s)$ |  Time $\Omega(n)$ |  Apriori info
---|---|---|---|---|---|---
FOCUS (Almuallim and Dietterich, 1991) | DT | ✓ |  |  | ✓ |
Relief (Kira and Rendell, 1992) | DT | ✓ |  | ✓ | ✓ |
Relief-F (Kononenko, 1994) |  PCC |  |  |  |  | ✓
INTERACT (Zhao and Liu, 2007) |  DT, SVM | ✓ |  | ✓ | ✓ |
Fisher (Gu et al., 2012) | k-NN | ✓ |  |  | ✓ |
MutInf (Zaffalon and Hutter, 2014) | NB | ✓ | ✓ | ✓ | ✓ |
Joint MutInf Maximization (Bennasar et al., 2015) |  NB, k-NN | ✓ | ✓ |  |  |
Interaction Weight-based FS (Zeng et al., 2015) |  DT, IB1, PART | ✓ |  | ✓ | ✓ | ✓
Infinite FS (Roffo et al., 2015) | SVM | ✓ | ✓ | ✓ | ✓ |
MultiSURF (Urbanowicz et al., 2018a) | - |  |  |  | ✓ | ✓
Table 1: A comparison table of evaluation metrics used in Feature- Ranking or
Selection paper proposals. Many different evaluation metrics are used,
illustrating there exist no consensus on definite meaningful evaluation
metrics in the field.

Stability, on the other hand, is not widely used in theoretical or
quantitative argumentation. The stability is defined as the ability of an
algorithm to produce consistent results given small changes in the sample
data. In this context, this can be phrased as the sensitivity of the feature
selector to data perturbations. Even though stability was recommended as a
relevant evaluation metric by Chandrashekar (2014), not many papers explicitly
argue for the stability of their method; the metric is called an “overlooked
problem” in Chandrashekar (2014). In many papers this metric is still regarded
as a future work for solidifying any experimental results - the development of
algorithms that achieve both high classification accuracy and high stability
is still seen as a ‘challenging’ problem by Tang et al (2015) (Tang et al.,
2014).

The trend seems to be turning though, with more authors becoming aware of the
importance of stability. Our reliance on machine learning is ever-increasing,
and so does the demand for interpretability and reliability of the algorithms.
Take for example a biomedical application, in which feature selection is used
to select genes in a gene sequencing analysis. Any expert in this domain field
will feel more confident if an algorithm produces stable results given a
varying sample population. In this way, algorithm stability is of much
relevance to real-world applications of ML. Stability has been long taken into
account into prediction tasks, but not so much in feature selection - see for
example Table 1.

Scalability is another point of interest that only recently caught more
attention. The extra demand of algorithms to allow for parallel execution has
been imminent as data grew tremendously large. Even, multi-core processing can
lack in terms of performance, hence introducing the need for algorithms that
can run in distributed fashion. Distributing a dataset over multiple machines
poses challenges to some existing methods, though. Some current methods
require the full dimensionality of the dataset to be available in-memory (Tang
et al., 2014). Yet, other methods require each sample to be visited multiple
times, e.g., to apply a sample re-weighting strategy in order to converge. For
these reasons, distributing any dataset workload on to multiple workers is a
non-trivial problem; no generalized solution exists for cutting the dataset
into chunks.

It is up to individual algorithms to find suitable ways of supporting parallel
solutions and more importantly, support cases where data is too large to fit
in-memory, i.e., apply distributed computing. Recent strategies overcome the
issue of working with fewer samples by retaining only those samples that are
most representative of the data - eliminating the need for working with a full
sample population. Although few in number, there exist proposals for
distributed dimensionality reduction methods (Li et al., 2020), using divide-
and-conquer techniques. Aggregating disjoint results would make for a
performance similar to that of a centralized solution.

###  3.2 Synthetic datasets and apriori information

Synthetic datasets are employed in many papers in the literature. Whereas some
datasets are injected by synthetically generated probe variables, others use
completed generated benchmark datasets, such that benchmarks can be conducted
in an even more controlled environment. About such partially synthetic
datasets has been spoken of in literature since long, using datasets that are
real but altered by injecting more data. In (Guyon and Elisseeff, 2003) the
authors speak of ‘probe variables’, which are used to discard any variable
that scores lower than any of the probes. If the probe variables are set to be
random variables, a simple way is obtained to apply a threshold to cut off
features from the selected feature subset, i.e., by introducing known noise
into the dataset we can construct more thoughtful cut-off thresholds.

Completely synthetic datasets, on the other hand, can allow for more
sophisticated metrics to be used. Possibilities for new evaluation metrics
are, for example described in (Solorio-Fernández et al., 2020): “Evaluation in
terms of the redundancy of the selected features” and “Evaluation in terms of
the correctness of the selected features”, the latter of which requires us to
know what features are informative a priori \- which is accomplished with
synthetic generation. Controlling all facets relevant to the quantitative
analysis manually makes for a Simulation study, which is argued for in
(Urbanowicz et al., 2018b) as follows:

> “Simulation studies such as these facilitate proper evaluation and
> comparison of methodologies because a simulation study can be designed by
> systematically varying key experimental conditions, and the ground truth of
> the dataset is known i.e. we know which features are relevant vs.
> irrelevant, we know the pattern of association between relevant features and
> endpoint, and we know how much signal is in the dataset.”

Indeed, there seems to be a trend toward including synthetically generated
datasets in experiments. In a review paper (Bolón-Canedo et al., 2013) the
authors argue that synthetically generated datasets can yield statistically
sound results because of the fact no inherent noise or redundancy will
obstruct the experiment process. In other papers simulation studies are
conducted as well (Cai et al., 2020) (Tang et al., 2020) (Li et al., 2020),
concluded by a small section depicting a ‘real data’ analysis to conclude the
point. For these reasons, a recommendation is made to include simulation
studies in any comprehensive benchmark on feature ranking methods.

###  3.3 The gap in the current literature

With respect to the above summarized works, there are some aspects missing in
the evaluation process. Because this research is set out to fill in the
highlighted missing parts, it is important to get a clear idea of the entire
set of missing aspects. Therefore, the concerned literature ‘gap’ is
summarized as follows.

Many papers in the literature evaluate Feature Ranking and Feature Selection
algorithms in different ways. Many validation estimators are used across
papers, causing the results to become incomparable. Moreover, not every paper
evaluates the stability or scalability of the algorithm, like shown in Table
1. Whereas the stability is a quantification of the algorithm’s robustness
against data permutations, the scalability means both storage- and time
complexity. Lastly, there exist opportunities for systematic evaluation using
apriori information on the feature importances. Few authors currently utilise
this opportunity.

Therefore, the above problems are to be addressed in this paper. A concrete
quantification of the algorithm’s stability, scalability and performance is
desired. Thereby, also synthetic data is considered, in which the relevant
features are known apriori. This paper also fills in the gaps left by some
papers which do only theoretically describe their experimental setup: the new
methodology was implemented in a pipeline, available as open-source software.
However, first, a look is taken at how to make Feature Rankings at all.

##  4 Methods for Feature Ranking

In the following, general theory related to the construction of feature
rankings is discussed. The theory is required to be discussed because in order
to best understand the evaluation process, an understanding of the
construction process is a must.

###  4.1 Terminology

Among the subject of reducing dataset dimensions, there exists a common
terminology that is used among the literature. Whilst some terminology is
synonymous, other seemingly related terms mean different concepts entirely.

Feature Selection and Feature Ranking are two terms often used
interchangeably. Since often feature selection is done by first constructing a
feature ranking and then cutting off features ranked lower than some threshold
$\epsilon$, the construction of a feature ranking is often times an integral
part of feature selection. Although feature selection methods exist that do
not construct feature rankings, the two are synchronous in many ways, and are
in this paper thus related to one another.

The term Feature Selection is used to indicate the general process of
obtaining a feature subset with reduced size without transforming the data.
Note that any feature selection method might transform the data in the
algorithm however it likes internally - the stated terminology is only
concerned with the eventual output of the feature selection algorithm. In this
paper, a broad perspective is taken and not only feature selection methods but
feature ranking methods generally are considered.

Feature Ranking is a broader term compared to Feature Selection, mapping onto
more domains than just Feature Selection. Because in the construction of a
Feature Ranking no assumptions are made on the desired data subset, a score is
assigned to each of the dataset features, giving each dimension an
‘importance’ score. Such scores can also be used to interpret and clarify a
Machine Learning model: thus making such ranking processes useful to the
interpretable AI domain.

Feature projection and Feature Selection are both processes relating to the
concept of dimensionality reduction (Cunningham, 2007), however, there exists
an important distinction between them. Whilst in the process of feature
selection, relevant dimensions are sought and selected without altering their
input values, in the process of feature projection (also called feature
extraction) data transformations are applied, mapping the original data onto a
lower-dimensional space. Common methods of feature projection are Principle
Component Analysis (PCA) for the supervised case and Linear Discriminant
Analysis (LDA) for the unsupervised case. Both PCA and LDA take a statistical
approach to detecting feature interactions, which not always results in an
optimal feature set for prediction. Rather, machine learning techniques can be
used to select a more optimal subset. Although the two methods are different,
the two aim at achieving the same goal and are thus encountered in similar
contexts.

###  4.2 Types of Feature Rankings

To better understand in what form a Feature Ranking might be defined and how
it relates to the process of Feature Selection, exact mathematical definitions
of various types of feature rankings are given first. At all times, the
prediction- or estimation of any quantity is denoted with a hat notation, i.e.
‘^’, whilst the ‘true’ value of the quantity has no such hat.

####  4.2.1 Feature importance

Feature importance scores are defined as a vector of $p$ dimensions,
containing real-valued numbers. Let us define the vector like so:

$\hat{\boldsymbol{w}}=(\hat{w}_{1},\hat{w}_{2},\ldots,\hat{w}_{p-1},\hat{w}_{p}),$
(1)

where $\hat{\boldsymbol{w}}\in\mathbb{R}^{p}$. Such a ranking is obtained, for
example, by running a feature ranker on the dataset and having it assigned as
score to each dimension. The vector is assumed to be normalized, i.e. each
vector element is divided by the vector sum:

$\hat{\boldsymbol{w}}=\frac{\hat{\boldsymbol{r}}}{\sum^{p}_{i=1}\hat{r}_{i}},$
(2)

given a scoring vector $\hat{\boldsymbol{r}}$, which indicates feature
importance on an arbitrary scale. It is self-evident that
$\hat{\boldsymbol{w}}$ has the property that $\sum^{p}_{i=1}\hat{w}_{i}=1$,
i.e. is a probability vector. An example such vector can be:

$\hat{\boldsymbol{w}}=(0.20,0.8,0.0),$ (3)

in which it is clear that the ranking algorithm found the second feature to be
the most important. In the case where multiple feature importance vectors are
considered, e.g., in the case where $B$ bootstraps (Section 5.1.2) are
considered, the vectors are stacked in a matrix, i.e.:

$\mathbf{\hat{W}}\in\mathbb{R}^{B\times p},$ (4)

which denotes $B$ feature importance $p$-dimensional vectors arranged in a
matrix of reals. The ground-truth $w$ will remain a vector, since there will
still be only one such vector available per dataset.

####  4.2.2 Feature support

Feature support indicates whether certain dimensions are chosen to be included
a feature subset; i.e. the vector marks elements as being chosen by the
feature selection process. In this definition, the feature support vector is
synonymous with the definition of a feature subset. Although some feature-
ranking and selection processes approximate a suitable feature support vector
directly, an algorithm can also make use of a threshold $\epsilon$ to generate
a feature support vector from a feature importance vector. The feature support
vector is a boolean-valued vector of $p$ dimensions.

$\hat{\boldsymbol{s}}=(\hat{s}_{1},\hat{s}_{2},\ldots,\hat{s}_{p-1},\hat{s}_{p}),$
(5)

where $\hat{\boldsymbol{s}}\in\mathbb{B}^{p}$. Note $\mathbb{B}$ is the
boolean-valued vector space, i.e. its elements lie in the set $\\{0,1\\}$. An
example such vector can be:

$\hat{\boldsymbol{s}}=(1,1,0),$ (6)

which is the feature support mask obtained from thresholding the feature
importance vector $\hat{\boldsymbol{w}}$ in the above example (Equation 3)
using the threshold $\epsilon>0.0$, causing one feature to be dropped from the
feature subset. Just like for the feature importance vector, the feature
support vectors can be arranged in a matrix in the boolean space, like so:

$\mathbf{\hat{S}}\in\mathbb{B}^{B\times p},$ (7)

given $B$ feature support vectors.

A sparse representation can also be constructed, given a feature support
vector. A set is created containing only the indices of the selected features,
causing a more sparse representation of the feature subset in case of high
dimensionality and relatively few selected features. The support vector
$\boldsymbol{s}$ and its prediction $\hat{\boldsymbol{s}}$ can be readily
converted back- and forth into such a sparse representation. We define the
sparse representation as the set $\hat{\mathbb{S}}$:

$\hat{\mathbb{S}}=\\{i\mid i\in\mathbb{Z}\wedge\hat{s}_{i}=1\\},$ (8)

where $|\hat{\mathbb{S}}|=d$, i.e. $d$ dimensions were selected in the feature
subset. Note, that the sparse feature subset is represented as a set instead
of a vector, meaning that the sequence is no longer considered ordered. To
show a concrete example, the vector from Equation 6 is converted to the
following set:

$\hat{\mathbb{S}}=\\{2,1\\}\text{ where }|\hat{\mathbb{S}}|=2,$ (9)

containing the indices of the selected features as defined in Eq 6, in no
particular ordering whatsoever. A sparse set $\hat{\mathbb{S}}$ can be
converted back to a feature support vector $\hat{\boldsymbol{s}}$ like so:

$\hat{s}_{i}=\begin{cases}1&\text{if }i\in\hat{\mathbb{S}}\\\
0&\text{otherwise}\end{cases}$ (10)

Denote $B$ such sparse feature support sets as $\hat{\mathbb{S}}^{boot}$:

$\hat{\mathbb{S}}^{boot}=\\{\hat{\mathbb{S}}_{1},\hat{\mathbb{S}}_{2},\ldots,\hat{\mathbb{S}}_{B-1},\hat{\mathbb{S}}_{B}\\},$
(11)

meaning $B$ sparse support sets were arranged in the superset
$\hat{\mathbb{S}}^{boot}$.

####  4.2.3 Feature rankings

Feature rankings are similar to feature importance scores, but with less
precision. Whereas in a feature importance vector each element is approximated
with a real-valued score, in an ordinary feature ranking the only considered
facet is the order of the features in terms of importance. Although in most
cases a feature importance vector is constructed first, after which a support
or ranking vector can be constructed, in some cases only a ranking is
available - e.g. in the case of RFE. A feature ranking is constructed by
assigning each dimension a rank, anywhere in the integer set
$\\{1,2,\ldots,p-1,p\\}$. Such, the vector can be expressed as:

$\hat{\boldsymbol{r}}=(r_{1},r_{2},\ldots,r_{p-1},r_{p}),$ (12)

where $\hat{\boldsymbol{r}}\in\mathbb{Z}^{p}$, the integer-space. An example
such vector can be:

$\hat{\boldsymbol{r}}=(2,3,1),$ (13)

which is the feature ranking obtained from converting the feature importance
vector $\hat{\boldsymbol{w}}$ in the above example (Eq 3) to a ranking. Such
rankings are easily converted to importance vectors using Equation 2, allowing
one to use the same statistical machinery as for feature importance vectors. A
higher rank number means a greater importance.

Such an integer-valued ranking can be utilised to select some $k$ best
features to use in a subsequent learning task, i.e. to perform a feature
subset selection. This is similar to selecting features using a feature
importance vector. Now, however, a subset is constructed not based on a
threshold value $\epsilon$, but by including a certain number of best features
in the subset. Even though the feature- importance and ranking vectors both
have the ability to select feature subsets, the feature importance vector
carries more meaning, because it more precisely quantifies the relative
importance of each feature. This extra meaning can be made to good use during
the evaluation process of the feature rankings, especially when using
synthetically generated datasets, in which the ground-truth feature relevance
is available.

###  4.3 Feature Selection Taxonomy

To better understand what type of feature selectors exist, and how they relate
to the types of feature rankings that can be constructed, a taxonomy is
considered. Feature Selection can be constructed by running a separate
statistical operation on the dataset, before running any learning algorithms,
or as part of a learning algorithm itself. In some cases, a learning algorithm
that is itself very sophisticated and time-consuming might still be worthwhile
to use as a feature selection pre-processing step. This is because if time is
saved by having the prediction estimator process less features, one might
still have to spend less time in his learning process. In this way, one might
enjoy gains in processing efficiency by using a feature selector.

Due to this distinction in the approach used in a feature ranking algorithm, a
subdivision can be made to separate methods into more specific categories. As
such, a common taxonomy in the field is created: subdividing feature ranking
methods into the categories of filter-, wrapper- and embedded methods
(Chandrashekar and Sahin, 2014).

####  4.3.1 Filter methods

Filter methods use some scoring mechanism to compute ‘usefulness’ for each
feature, without applying a learning algorithm. Having applied some ranking
criterion, often a feature ranking is produced, after which some thresholding
operation can be applied to select features. Although filter methods are often
computationally light and do not overfit due to the absence of a learning
algorithm, filter methods might miss out on more complex feature interactions,
causing a non-optimal subset as a result. Also, choosing a suitable threshold
to use can be difficult.

Examples of filter methods include the Fast Correlation Based Filter (Yu and
Liu, 2003), which uses a quick and easy-to-compute statistical measure to
select features according to some predefined threshold $\epsilon$ (in Yu et
al. denoted as $\delta$). To illustrate the type of statistical quantities
generally used in filter methods, the statistical quantity in Yu et al (2003)
can be denoted like so:

$SU(X,Y)=2\left[\frac{IG(X\mid Y)}{H(X)+H(Y)}\right],$ (14)

where $IG(X|Y)$ is the information gain between two random variables $X$ and
$Y$ and $H(X)$ is the entropy of a variable: which are both metrics coming
from the information-theoretical domain. The measure $SU(X,Y)$, then, is the
symmetrical uncertainty of two features $X$ and $Y$, where a feature $Y$ is
more correlated to $X$ than to $Z$ if $IG(X,Y)>IG(Z,Y)$. In this way, a
ranking can be constructed considering the measure $SU(X,Y)$, where feature
are sought with high $SU$ scores. In a second phase of the algorithm, the
scoring table is traversed yet again, to eliminate possible redundant features
included in the selected feature subset.

####  4.3.2 Wrapper methods

Wrapper methods, on the other hand, use some learning algorithm to determine a
suitable subset of features. A search is conducted over the space of possible
feature subsets, eventually selecting the subset that has the highest
validation score in the test set using a chosen learner as a predictor.
Characteristics that define wrapper methods lend themselves similar
characteristics to traditional optimisation problems; although an exhaustive
search might yield an optimal solution, such a solution might not always be
feasible due to its great time complexity. For this reason, in some
applications a filter is applied first, before running a wrapper method.

Examples of such methods are numerous. Straight-forward methods include the
range of sequential feature selection methods, which aim to start (1) either
with the full subset of dataset features or (2) with an empty set of features.
The two approaches are called Forward Feature Selection and Backward Feature
Elimination, respectively. To then facilitate a forward- or backward iteration
step it is customary to use an estimator of some kind to retrieve a scoring on
the features, selecting either the best- or eliminating the worst scored
feature. RFE is one such backward-elimination method, which might, for
example, use a SVM to construct estimation scores (Maldonado and Weber, 2009).
Another option is to perform an exhaustive search, in which every feature
combination is tried such to obtain the optimal feature subset, given the
learning task and the estimator used. To obtain such a ranking, the chosen
estimator might use an arbitrary method to compute it, be it a relatively
simple learning step or a sophisticated model evaluation. Although such
methods can perform reasonably well, such methods tend to be more time-
consuming than filter methods.

####  4.3.3 Embedded methods

Embedded methods seek to combine the training task and feature selection.
Given some suitable learner, features are weighted during the training
process, producing either a feature ranking or a feature subset afterwards.
e.g. some learners compute feature importance scores as part of their training
process, which can then be used in combination with some threshold to select
relevant features. Having already trained the model, subsequent prediction
tasks can benefit from increased prediction speed by using less data.

Examples of embedded methods are Least Absolute Shrinkage and Selection
Operator (LASSO) and Ridge Regression, which perform feature selection
alongside the process of finding optimal regression coefficients. For further
explanations, see Section 4.5.1. Another notable embedded method is a DT,
which constructs, inherent during its fitting stage, a measure of importance
on each variable. It does so by computing the probability of reaching a
certain node - and determining the decrease in node impurity caused by
weighting this probability value. To facilitate this computation, the
probability of reaching a certain node is computed by considering the number
of samples that reach the node during its decision-phase, divided by the total
number of samples. Such, the leaves and depth obtained during the construction
phase of the DT can be used to determine a measure of feature importance,
‘embedded’ into its learning phase.

####  4.3.4 Hybrid methods

A hybrid method is any method that is not classifiable by a single category,
but rather lends from multiple categories. Hybrid methods can, for example,
combine filter and wrapper methods (Hsu et al., 2011), by first applying a
computationally efficient filter and refining the result by using a wrapper
method. Another paper (Das, 2001) describes its approach as hybrid due to both
adding- and removing features in the feature selection process - exhibiting
both forward- and backward selection at the same time. Lately, research was
also put into examining Ensemble feature selection methods (Bolón-Canedo and
Alonso-Betanzos, 2019), which combines the outputs of multiple selectors and
decides useful features accordingly using some voting committee. Ensemble
methods can be seen as hybrids or are seen as a category on its own.

###  4.4 Types of features

An important consideration in choosing a suitable feature selection method for
any task is what kind of structure the concerning data has, if any at all.
Data might exhibit tree, graph, or grouped structures, which is essential
information when detecting feature interactions and determining useful
features. Support for structured data is relatively new in the field and has
not been an extensive point of concern for many feature selection algorithms
in the past. Many traditional feature selection algorithms focused primarily
‘conventional flat’ features (Li et al., 2017), in which the assumption is
made that the data is independent and identically distributed (i.i.d.).
However, this assumption is widespread among many machine learning algorithms,
since in many applications datasets are normalized to fit the i.i.d. condition
before they are used.

Conventional data are opposed to more complex data structures, i.e., datasets
with ‘structured features’, as coined by the proposed taxonomy from Li et al
(2017) Li et al. (2017), but also to linked data and streaming data. In
streaming data, the quality of an initially selected feature subset can be
improved as more data comes in, and a feature selection algorithm is able to
benefit from a larger distribution of samples. Adapting existing feature
selection algorithms to fit the demands of streaming data proved to be a
nontrivial problem. Nowadays, many companies and institutions have to deal
with data volumes that easily exceed the boundaries of in-memory storage
capacity - limiting the data scientist to train on only subsets of the entire
datasets. Smart sampling is therefore needed to retain a representative sample
distribution.

The scope of this research is limited to only the most common type of
features, conventional-, flat features (i.i.d.), or also tabular data.

###  4.5 Constructing Feature Rankings

Although the number of existing feature ranking methods are numerous, a small
set of example methods are explored to get a better understanding of the range
of methods that do exist. First of all, a way of constructing a feature
ranking that originates from classical statistics is examined. Next, a more
algorithmic-type of approach is examined, which was specifically designed for
the feature selection domain.

####  4.5.1 Regularized Linear Regression

One of the most fundamental methods in statistics is linear regression. It can
be solved both analytically and numerically: in which the optimal approach is
dependent on the amount of dataset dimensions at hand - where the amount of
dataset dimensions $p$ gets very large, the analytic solution gets slower
compared to an approximate method like Stochastic Gradient Descent. Recall
that we can analytically solve Linear Regression by minimizing the Residual
Sum-of-Squares cost function (Hastie et al., 2009):

$\text{R}(\boldsymbol{\beta})=(\mathbf{Y}-\mathbf{Z}\boldsymbol{\beta})^{\intercal}(\mathbf{Y}-\mathbf{Z}\boldsymbol{\beta}),$
(15)

in which $\mathbf{Z}$ is our design matrix. Regression using this loss
function is also referred to as ‘Ordinary Least Squares’. The mean of the cost
function R over all samples is called Mean Squared Error, or MSE. Our design
matrix is built by appending each data row with a bias constant of 1 - an
alternative would be to first center our data to get rid of the intercept
entirely. To now minimize our cost function we differentiate R with respect to
$\boldsymbol{\beta}$, giving us the following unique minimum:

$\hat{\boldsymbol{\beta}}=(\mathbf{Z}^{\intercal}\mathbf{Z})^{-1}\mathbf{Z}^{\intercal}\mathbf{Y},$
(16)

which results in the estimated least-squares coefficients given the training
data, also called the normal equation. We can classify by simply multiplying
our input data with the found coefficient matrix:
$\hat{\mathbf{Y}}=\mathbf{Z}\hat{\boldsymbol{\beta}}$. Now, in the case where
our model is fit using multiple explanatory variables, we are at risk of
suffering from multicolinearity \- the situation where multiple explanatory
variables are highly linearly related to each other causing non-optimal
fitting of the model coefficients.

In Ridge regression, we aim to tamper the least squares tendency to get as
‘flexible’ as possible to fit the data best it can. This might, however, cause
parameters to get very large. We therefore like to add a penalty on the
regression parameters $\boldsymbol{\beta}$; we penalise the loss function with
a square of the parameter vector $\boldsymbol{\beta}$ scaled by new
hyperparameter $\lambda$. This is called a shrinkage method, or also:
regularization. This causes the squared loss function to become:

$\text{R}(\boldsymbol{\beta})=(\mathbf{Y}-\mathbf{Z}\boldsymbol{\beta})^{\intercal}(\mathbf{Y}-\mathbf{Z}\boldsymbol{\beta})+\lambda\boldsymbol{\beta}^{\intercal}\boldsymbol{\beta},$
(17)

where we can see that $\boldsymbol{\beta}^{\intercal}\boldsymbol{\beta}$
denotes the square of the parameter vector, thus supplementing the loss
function with an extra penalty. This is called regularization with an $L^{2}$
norm; which generalization is called Tikhonov regularization, which allows for
the case where not every parameter scalar is regularized equally. If we were
to now derive the solutions of $\boldsymbol{\beta}$ given this new cost
function by differentiation w.r.t. $\boldsymbol{\beta}$:

$\hat{\boldsymbol{\beta}}^{\text{ridge
}}=\left(\mathbf{Z}^{T}\mathbf{Z}+\lambda\mathbf{I}\right)^{-1}\mathbf{Z}^{T}\mathbf{Y},$
(18)

in which $\lambda$ will be a scaling constant that controls the amount of
regularization that is applied. Note $\mathbf{I}$ is the $p\times p$ identity
matrix - in which $p$ are the amount of data dimensions used. An important
intuition to be known about Ridge Regression, is that directions in the column
space of $\mathbf{Z}$ with small variance will be shrunk the most; this
behavior can be easily shown be deconstructing the least-squares fitted vector
using a Singular Value Decomposition (SVD).

LASSO regression is a slightly modified variant of Ridge Regression, where
instead of an $L^{2}$ norm an $L^{1}$ norm is used instead. This problem can
be denoted as the following minimization problem:

$\hat{\boldsymbol{\beta}}^{\text{lasso
}}=\underset{\boldsymbol{\beta}}{\operatorname{argmin}}\left\\{\frac{1}{2}\sum_{i=1}^{N}\left(y_{i}-\beta_{0}-\sum_{j=1}^{p}x_{ij}\beta_{j}\right)^{2}+\lambda\sum_{j=1}^{p}\left|\beta_{j}\right|\right\\},$
(19)

in which the similarities between Ridge Regression are easily seen. The
squared $\sum_{1}^{p}\beta^{2}_{j}$ ridge penalty is replaced by by the
$L^{1}$ norm penalty of $\sum_{1}^{p}|\beta_{j}|$.

Feature selection can be employed using both LASSO\- and Ridge regression.
Because the dimensions whose coefficients are shrunk the most are presumably
the least relevant to the prediction task, non-contributing features can be
cut off from the design matrix using a threshold point. In fact, because LASSO
does not square the weights vector but takes the absolute value, some
coefficients might even be shrunk to near-zero values: removing the need for
defining a threshold at all, since the coefficients have zero contribution
already. In line with the feature subset definition given in Section 4.2.2,
such a feature support set can be constructed by thresholding the LASSO
coefficients like so:

$\hat{\mathbb{S}}^{lasso}=\\{j\mid\lvert\beta_{j}\rvert\geq\epsilon\\},$ (20)

which will result in a feature support set containing only the indices of
features where the coefficients were at least larger than $\epsilon$.

Although the former methods are suited for regression tasks only,
regularization schemes are employed widely and have similar mechanics in these
applications. Examples are numerous and include regularized Logistic
Regression, regularized SVMs and regularized Neural Networks: employing
regularization in any Machine Learning model is standard practice. In this
way, we gained insight into a fundamental tool to estimate feature importance
and shrinking the model coefficients using a regularization term - and more
importantly, the fact that it can be employed for feature selection.

####  4.5.2 Relief-Based Feature Selection algorithms

The Relief-family of feature selection algorithms originates from a seminal
paper by Kira et al. (Kira and Rendell, 1992), introducing the original
version of the Relief algorithm. Over time, many variations on the Relief
algorithm have been made, most notably ReliefF (Kononenko, 1994). In fact, so
many variations on the algorithm were made that one can speak of Relief-Based
Feature Selection Algorithmss (RBAs) in the literature (Urbanowicz et al.,
2018a). Due to the overall architectural design of the algorithm, RBAs manage
to rank features and optionally select a feature subset within a reasonable
time-complexity domain and is often competitive when it comes to validation
estimator performance.

The algorithm works by usage of an instance-based learning method, considering
one data sample at a time and updating statistics on feature importance. Given
a training dataset $\mathbf{Z}$, $n$ samples, $p$ dimensions and a relevancy
threshold $\epsilon$, Relief can compute a measure of feature importance in a
finite amount of time. The essential concept is to compute the $p$-dimensional
Euclidean distance for a randomly picked instance, iterating over all dataset
samples and computing the nearest ‘miss’ and nearest ‘hit’ for every instance.
Such a miss- or hit is defined to be either one of the positive- or negative
dataset instances. To now compute how much ‘relevance’ should be added for
each feature given the closest hit- and miss, a subroutine is used, in which
for every feature the diff between the randomly picked instance and nearest-
hit and miss is computed and squared.

Once the iterations are complete, the summed weights for each feature are
stored in a vector $\boldsymbol{w}^{relief}$. Once all $n$ samples have been
traversed, the feature weights vector $\boldsymbol{w}^{relief}$ is normalized
by dividing by the amount of samples $n$. In this way, a ranking is
constructed - which can be converted into a feature subset by applying the
threshold $\epsilon$ to the ranking. i.e., using the notation introduced in
4.2, a feature subset can be constructed like so:

$\hat{\mathbb{S}}^{relief}=\\{i\mid\boldsymbol{w}^{relief}_{i}\geq\epsilon\\},$
(21)

given some threshold $\epsilon$.

##  5 Evaluating Feature Rankings

Like seen in Chapter 3, the evaluation of feature- ranking and selection
methods is a nontrivial process, that has been conducted in many ways in the
literature. Because of this reason, it is desired to acquire a general
evaluation method that is applicable to all feature- ranking and selection
methods. The algorithms should be evaluated comprehensively: highlighting all
aspects relevant to the performance of the method. In this chapter, a
reasoning is given on which evaluation metrics are sensible to use, and which
evaluation methods are recommended to use for different scenarios.

###  5.1 Cross-Validation and Bootstrapping

Several steps can be undertaken to provide more reliable estimates of the
feature ranker performance. The goal is to prevent the practitioner to get
promising but misleading results, whilst keeping the complexity of the
benchmarking process at a reasonable level.

####  5.1.1 Cross-Validation

For almost any Machine Learning experiment, it is recommended to use some form
of Cross-Validation (CV). Because a model has the possibility to get ‘learn’
the training data, evaluating the performance of an estimator on just the set
of training data is dangerous practice. Evaluation metrics can be misleading,
suggesting better model performance than is actually the case. For example, an
improperly cross-validated model that is employed outside the realm of the
training data might turn out to have very poor generalization performance.
Thus, also in the context of benchmarking feature ranking and selection
algorithms, validation must be performed by holding out a set of the available
data samples at hand.

The most simple form of CV is a training/testing split, in which the
practitioner holds out one set of the example data for use in the testing
phase - it is to be remain unseen by the model during the training phase. More
robust methods include training the model multiple times on various training
datasets and evaluating on the held out testing data. This can be done using
5-fold or 10-fold CV, where, for example, in the case of 5-fold CV
$\nicefrac{{4}}{{5}}$th of the data is reserved for training and
$\nicefrac{{1}}{{5}}$th for testing, repeated for each fold - so 5 times.

What is especially important in the process of conducting CV, is that no
operations are performed before having split the data, that can influence the
experimental results. For example, in the scenario where one want to select
variables to be included in some feature subset and run some prediction
estimators afterward, a practitioner must be careful to split the data before
the variable selection process. This is because otherwise, the predictive
estimators might gain an unfair advantage due to having already seen the left
out samples, i.e. they have already gained an advantage from the variable
selection that was performed on the full dataset. This might cause skewed and
misleading estimations of the error rate of the estimators - causing faulty
models. Therefore, it is at all times important to keep in mind the ‘right’
way of performing CV: split first before running any operations that regard
the dataset samples (Ambroise and McLachlan, 2002).

####  5.1.2 Bootstrapping

Bootstrapping, on the other hand, is a similar but different process.
Performing a bootstrap is a classical method for estimating statistical
quantities regarding the learning process, such as variance, prediction error,
or bias. The process works by resampling the dataset with replacement $B$
times, such to create $B$ different permutations of the dataset. Then, the
learning process is repeated for each of the $B$ bootstrap datasets, i.e.,
refitting the estimators for each of the permuted datasets. If, for example,
the designated dataset is denoted like
$\mathbf{Z}=\left(z_{1},z_{2},\ldots,z_{b}\right)$ with each sample
$z_{i}=\left(x_{i},y_{i}\right)$, then the $b$-th bootstrap permutation of the
dataset can be denoted as $\mathbf{Z}^{*b}$. Such, estimates can be made of,
for example, the variance of some statistical quantity $S$ computed over the
dataset $\mathbf{Z}$:

$\widehat{\operatorname{Var}}[S(\mathbf{Z})]=\frac{1}{B-1}\sum_{b=1}^{B}\left(S\left(\mathbf{Z}^{*b}\right)-\bar{S}^{*}\right)^{2},$
(22)

which can be seen to be the unbiased average of the statistical quantity $S$
over the $B$ bootstrap permutations of $\mathbf{Z}$. Note that the average
value of the statistic $S$ is computed like so:

$\bar{S}^{*}=\sum_{b}S\left(\mathbf{Z}^{*b}\right)/B,$ (23)

i.e. by averaging $S$ over the $B$ bootstraps. Like such, a reasonable
estimate can be made of any statistical quantity computed over the dataset, as
long as the quantity can be computed for any permutation on the dataset and
enough computational resources are available to run $B$ repeated experiments
on each of the permutations.

The practice of bootstrapping comes useful to the evaluation of feature-
ranking and selection algorithms, allowing the practitioner to better estimate
quantities that would previously be less reliable estimates. Examples of such
quantities are the feature ranking algorithm stability, variance, or fitting
time. Especially in assessing the algorithm stability, there exists an
interest to know how the designated algorithm functions under conditions of
varying samples. For this, bootstrapping is especially useful, since the exact
data generating distribution at hand is often not available- meaning no more
samples can be drawn from the distribution to generate a larger data
population. In the lack of a data generating distribution, resampling with
replacement offers a solution.

###  5.2 Validation estimators

A straight-forward way to evaluate the performance of a feature selection
algorithm is to run the ranking algorithm, and subsequently run a ‘validation’
estimator on the selected feature subset. The quality of the selected feature
subset is then quantified through the performance of the validation estimator:
when more informative and relevant features are selected, validation estimator
performance presumably goes up. The validation estimator can be configured to
be any classifier or regressor, dependent on the prediction task at hand,
though often choices are made from a common set of estimators as can be seen
in Table 1.

This idea can be extended to feature rankings. Many feature selection
algorithms are in fact feature ranking algorithms and do not provide a built-
in mechanism to perform feature subset selection - the algorithms allow a user
to define the desired amount of features to be selected using a hyper-
parameter, which then uses its feature ranking internally to construct a
feature subset, like explained in Section 4.2. It is therefore desired to also
benchmark such algorithms in a systematic way: allowing more flexibility with
respect to the evaluation process.

Figure 3: The process of running a validation estimator on a feature subset.
In the case of a feature ranking, a validation estimator is fit for some
various amount of feature subsets. Often, feature subsets with some number of
the highest ranked features are evaluated.

As can be seen in Figure 3, the process of running a validation estimator is
straight-forward. Once a feature subset has been determined by thresholding a
feature ranking or using the feature subset computed by a feature selection
algorithm, a validation estimator is trained and evaluated on the held-out
test set. Afterwards, a suitable evaluation metric is used depending on the
learning task at hand: classification or regression. Note - the research scope
was confined to just classification and regression, not clustering. Exactly
how many validation estimators are fit to evaluate a feature ranking depends
on the ranking type (Section 4.2). In the case of a feature support vector,
only the feature subset might have to be evaluated, whilst for feature
importance or ranking vectors one might want to evaluate the first $k$ best
feature subsets. Generally, each of the vectors are evaluated like so:

  * •

Feature importance (Section 4.2.1). Fit $k$ validation estimators, for the
first $k$ best feature subsets, starting with a feature subset including only
the highest ranked feature and subsequently including lower-ranked features.
$k$ might be chosen depending on the dataset size at hand.

  * •

Feature support (Section 4.2.2). A feature subset was already selected by the
feature selection algorithm and such the feature subset can be directly
evaluated by the validation estimator.

  * •

Feature ranking (Section 4.2.3). A feature ranking can fit $k$ validation
estimators; similarly to the feature importance vector.

Notably, in the case where a feature ranking algorithm supports computing both
feature importance or ranking vectors and the feature support vector, simply
the $k$ best feature subsets are evaluated by the validation estimator. It is
thereby assumed that the selected feature subset is included in one of the $k$
best feature subsets. If for some reason the selected feature subset
$\hat{\mathbb{S}}$ is not in any of the $k$ best feature subsets, the feature
subset will be evaluated by the validation estimator separately.

A bootstrap can be used to generate more accurate estimations of the
validation process. Especially for validation estimators that are sensible to
varying permutations of the training data this process is useful - since
otherwise validation metrics could be particularly misleading. In this way,
given either an accuracy- or R2 score for classification and regression,
respectively, an averaged bootstrapped score can be computed. This is done
simply by plugging in either the R2 score or accuracy (or any other metric
suitable for regression- or classification) into Equation 23.

###  5.3 Apriori knowledge on relevant features

In the scenario where we know apriori which features are relevant, more direct
evaluation can be applied on the constructed feature ranking. Such apriori
ground-truth information can be gained in a couple of ways: (1) by utilising
human domain expert knowledge, (2) by augmenting real-world datasets with
noisy dimensions and assuming all other dimension to be uniformly relevant or
lastly (3) by generating synthetic datasets with known feature importance
levels (see Section 6.2.1. In the context of this section, the assumption is
made that the feature importance scores are simply ‘known’ apriori: no
restriction is made on exactly how feature relevance ground truth was
retrieved.

Knowing the relevant dataset features apriori can be useful information to the
feature ranking evaluation process. This is because even though the goal of
feature- ranking and selection algorithms is to separate the useful from the
irrelevant features, the performance of such a ranking is nowadays most
commonly evaluated by training yet another estimator; the ‘validation’
estimator (Section 3 and Section 5.2). This makes the evaluation scores
dependent on the validation estimator, requiring publications to have run the
same validation estimator to make the results comparable between papers. Also,
for some datasets, the chosen validation estimator might be sophisticated
enough to make up for any noisy dimensions that were added - reducing the
differences between the feature rankings, which is instead desired to be
amplified to determine which feature ranker is the best for which dataset.
Lastly, a sophisticated validation estimator might also perform a form of
feature selection on its own, which might also skew the performance of faulty
feature rankings that included lots of noisy dimensions in its feature subset.
For these reasons, it is interesting to investigate the possibilities to find
robust and reliable evaluation methods that do not require a validation
estimator.

Using the ground-truth relevant features, an accompanying evaluation metric
can be constructed. Taking into consideration how different feature ranking
algorithms construct various types of feature rankings as described in Section
4.2, a suitable metric can be created for each.

Feature importance (Section 4.2.1) scores are defined to be a real-valued $p$
dimensional vector, $\hat{\boldsymbol{w}}$. In order to evaluate the closeness
of the vector $\hat{\boldsymbol{w}}$ to the dataset ground truth
$\boldsymbol{w}$, a range of metrics can be used. A first approach would be to
regard the estimated vector as an estimation of a continuous target vector,
i.e., to regard the problem as a regression type of task. In such a
perspective, each predicted feature importance would pose as a data sample. In
this way, usual metrics related to the regression task can be used.

The R2-score is one such metric that can be used in this context. This can be
formalized like so:

$\displaystyle R^{2}$ $\displaystyle=1-\frac{\mathrm{RSS}}{\mathrm{TSS}}$ (24)
where $\displaystyle\mathrm{RSS}$
$\displaystyle=\sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)^{2}=\sum_{i=1}^{p}\left(w_{i}-\hat{w_{i}}\right)^{2},\text{
and}$ $\displaystyle\mathrm{TSS}$
$\displaystyle=\sum_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}=\sum_{i=1}^{p}\left(w_{i}-\bar{\boldsymbol{w}}\right)^{2},$

which can be used to evaluate the closeness of the predicted feature
importance vector $\hat{\boldsymbol{w}}$ to the ground-truth feature
importance vector $\boldsymbol{w}$.

The logistic-loss, or cross-entropy score is another metric which can be
employed to quantify the quality of a feature importance vector
$\hat{\boldsymbol{w}}$. In this scenario, however, the target variables are
regarded to be binary labels, which therefore means the predicted feature
importance vector $\hat{\boldsymbol{w}}$ is to be compared against
$\boldsymbol{s}$ instead of $\boldsymbol{w}$. Remember, from Section 4.2.2,
that $\boldsymbol{s}$ is the $p$-dimensional vector indicating with a boolean
whether or not the feature is relevant, i.e.
$\boldsymbol{s}\in\mathbb{B}^{p}$. Such, a zero indicates a feature is not
relevant; and a one indicates the feature is informative or relevant. In this
way, a feature ranker approximates the true binary relevance labels
$\boldsymbol{s}$ with probability values arranged in the probability vector
$\hat{\boldsymbol{w}}$. The logistic loss can be formalized like so:

$\displaystyle L_{\log}(y,p)$ $\displaystyle=-\log\operatorname{Pr}(y\mid
p)=-(y\log(p)+(1-y)\log(1-p))$ (25) $\displaystyle\text{substituting for
}s_{i}\text{ and }\hat{w_{i}}$ $\displaystyle L_{\log}(s_{i},\hat{w_{i}})$
$\displaystyle=-\log\operatorname{Pr}(s_{i}\mid\hat{w_{i}})=-(\hat{w_{i}}\log(\hat{w_{i}})+(1-s_{i})\log(1-\hat{w_{i}}))$
averaged over all vector components $\displaystyle
L_{\log}(\boldsymbol{s},\hat{\boldsymbol{w}})$
$\displaystyle=-\frac{1}{p}\sum_{i=1}^{p}(\hat{w_{i}}\log(\hat{w_{i}})+(1-s_{i})\log(1-\hat{w_{i}}))$

Feature support (Section 4.2.2) is evaluated differently. Just like in the
logistic-loss metric, the ground truth relevant labels are encoded as binary
labels, i.e. the vector $\boldsymbol{s}$ is used as the ground-truth. Now,
however, the predicted targets are also encoded as binary labels, i.e. the
vector $\hat{\boldsymbol{s}}$ is used. This means that the problem is to be
viewed in the supervised classification perspective and all metrics
accompanying such task are available.

Classification accuracy is such a metric, measuring simply the average number
of correct predictions were made between the predicted- and true target
labels. It can be defined as such:

$\displaystyle\operatorname{accuracy}(y,\hat{y})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{1}\left(\hat{y}_{i}=y_{i}\right)$
(26) $\displaystyle\text{substituting for }s_{i}\text{ and }\hat{s_{i}}$
$\displaystyle\operatorname{accuracy}(s,\hat{s})$
$\displaystyle=\frac{1}{p}\sum_{i=1}^{p}\boldsymbol{1}\left(\hat{s}_{i}=s_{i}\right),$

where $\boldsymbol{1}(x)$ is the indicator function (Davis, 2004). Using this
metric, the amount of useful features included in a feature subset is rewarded
with higher accuracy scores accordingly.

Feature rankings (Section 4.2.3) can be evaluated similarly to feature
importance scores - thereby also requiring normalization. Presume that besides
the relevance of the features is known as a binary value, also the order of
relevance is known, i.e. which features are more relevant than others. This is
very similar to the feature importance scores: though the difference is that
the feature importance scores must first be converted to a ranking such to
allow for meaningful comparison. In this reasoning, both the predicted- and
the ground truth feature ranking vectors, which are $\hat{\boldsymbol{r}}$ and
$\boldsymbol{r}$ respectively, can be normalized using Equation 2 to obtain
$\hat{\boldsymbol{w}}$ and $\boldsymbol{w}$, respectively.

In this way, the normalized vectors $\hat{\boldsymbol{w}}$ and
$\boldsymbol{w}$ can again be considered probability vectors, such that
metrics like the R2-score can be used: just like for the feature importance
scores.

One has to keep in mind, however, not to intermix the feature importance and
feature ranking scores with each other: even though the same metric is used to
convert to summarize the predicted feature- importance and ranking vectors
into a single scalar, the scorings are built from different vectors to begin
with. The R2-score coming from the feature importance vectors might have an
unfair advantage due to the fact that they have floating-point precision on
their approximations, whilst the feature ranking vectors $\boldsymbol{r}$ are
converted from the integer domain to be normalized into floating-point
numbers.

###  5.4 Stability

The stability of any algorithm is an important facet of the total method
performance - which must not be overlooked or forgotten. This is for in many
applications, only robust algorithms can be systematically relied on. Although
it is only natural for algorithms to vary in fitting behavior under different
permutations of the sample set, it is at all times desired to get an algorithm
that is as stable as possible. To evaluate the stability of feature ranking
algorithms, an easy method is to take several bootstrap permutations of the
dataset (Section 5.1.2), to simulate drawing new samples from the data
generating distribution at hand. Assuming such a bootstrapping procedure to be
in place, several metrics can be used to quantify stability.

####  5.4.1 Stability of feature importance vectors

Feature importance scores are defined to be $p$ dimensional vectors in the
domain of real numbers $\mathbb{R}$ (Section 4.2.1). Accordingly, the matrix
containing $B$ such vectors is the $B\times p$ dimensional matrix in
$\mathbb{R}$, denoted as $\hat{\mathbf{W}}$. A straight-forward way to assess
the stability of the feature importance matrix is to compute the variance for
each of the $p$ dimensions over $B$ bootstraps, i.e. Eq 22 is used by plugging
in a column of the feature importance matrix $\hat{\mathbb{W}}$ as the
statistical quantity $S$:

$\widehat{\operatorname{Var}}[\hat{\mathbf{W}}_{:,i}]=\frac{1}{B-1}\sum_{b=1}^{B}\left(\hat{W}_{b,i}-\bar{\hat{\mathbf{W}}}_{:,i}^{*}\right)^{2},$
(27)

where $\bar{\hat{\mathbf{W}}}_{:,i}^{*}$ is the average feature importance
score for the $i$-th feature over $B$ bootstraps. In this way, the variance
over each dimension can be computed for a single feature ranker.

To summarize the variances over all dimensions, their summation might be
taken, i.e., the variances are summarized as the scalar
$\sum_{i=1}^{p}\widehat{\operatorname{Var}}[\hat{\mathbf{W}}_{:,i}]$.
Alternatively, one might want to weight such summation, in order to reflect
the desire to penalize instabilities in the ranking of relevant features more
so than instabilities in the ranking of irrelevant features. In the case where
the ground-truth feature importance scores are given, i.e. $\boldsymbol{w}$ is
known, this vector might be used to apply a weighting scheme to the summation
of the variances. One can do so by taking the inverse normalized ground truth
feature importance score, i.e. $\frac{1}{w_{i}}$ for the $i$-th feature. Such,
the summarization of $\hat{\mathbf{W}}$ now becomes
$\sum_{i=1}^{p}\frac{1}{w_{i}}\widehat{\operatorname{Var}}[\hat{\mathbf{W}}_{:,i}]$.

####  5.4.2 Stability of feature support vectors

This time around, a way to quantify the stability of a feature subset is
desired, i.e., the feature support vector $\hat{\boldsymbol{s}}$. A measure
for the quantifying stability of set permutations that comes easily to mind
might be the Hamming distance between multiple pairs of algorithm runs, i.e.
ran on the same dataset with varied sample populations. Indeed, in (Dunne et
al., 2002) a measure based on Hamming distance is proposed. Reports also exist
on numerous other approaches, like an entropy based measure (Křížek et al.,
2007), a measure based on the cardinality of intersection measure (Kuncheva,
2007) and a measure based on correlation coefficients (Kalousis et al., 2007).

Desired properties of any measure are important to define in concrete manner.
To compare the strength of any of these measures, the general objective for
exactly what information we want to convey in a stability metric should be
verbalized first. In (Mohana Chelvan and Perumal, 2016), three desired
stability metric properties are expressed: (1) Monotonicity, (2) Limits and
(3) Correction for chance. Given the proposed measures in the literature, only
few measures satisfied all properties. In another paper (Nogueira et al.,
2018), however, the authors extend the set of desired properties to include
another two: (4) the stability estimator must be Fully defined and (5) Maximum
Stability $\leftrightarrow$ Deterministic Selection \- meaning that a maximum
stability value should be achieved if-and-only if all feature sets are exactly
identical.

A measure for quantifying feature subset stability was proposed in (Nogueira
et al., 2018), using the set of newly proposed desired properties. Having
explored a statistically sound method to define a metric satisfying all five
properties, the authors proposed a novel stability estimator, adjusted to this
paper’s terminology:

$\hat{\Phi}(\hat{\mathbb{S}}^{boot})=1-\frac{\frac{1}{p}\sum_{i=1}^{p}\sigma_{i}^{2}}{\mathbb{E}\left[\frac{1}{p}\sum_{i=1}^{p}\sigma_{i}^{2}|H_{0}\right]}=1-\frac{\frac{1}{p}\sum_{i=1}^{p}\sigma_{i}^{2}}{\frac{\bar{k}}{p}\left(1-\frac{\bar{k}}{p}\right)},$
(28)

where $\hat{\Phi}$ resembles the stability estimate of the feature subsets
arranged in $\hat{\mathbb{S}}^{boot}$ (Eq 11). Remember, that each set in
$\hat{\mathbb{S}}^{boot}$ is a feature subset containing the indices of the
selected features. Each of the sets in $\hat{\mathbb{S}}^{boot}$ is created by
running the feature ranker on a permutation of the datasets, i.e. a resampling
of a dataset equal probability distribution and dataset properties.
Furthermore, the authors define $\sigma_{i}^{2}$ as the unbiased sample
variance of the selection of the $i^{th}$ feature and $\bar{k}$ as the average
number of features selected over the $B$ feature sets. Given this paper’s
goal, there is no necessity for going into further mathematical details - for
this we refer to the paper itself.

What is important is that this new measure satisfies all desirable properties
for quantifying stability, as was proven in the paper. Accompanying the novel
definition are instructions for computing confidence intervals and for
performing a hypothesis testing for comparing various method stabilities.
Exploring feature selection stability values given a set of parameter choices
not only allows for choosing better hyperparameters, but also allows for
comparing stabilities over various feature selection methods.

####  5.4.3 Stability of feature ranking vectors

Lastly, the stability of feature ranking vectors is also to be quantified.
Since the stability measure for feature support was all about working with
sets instead of vectors, the best option is to again normalize to a feature
importance vector and compute the variance thereof. The summarization of the
variances might again be weighted using the ground-truth feature importance
vector $\boldsymbol{w}$, if it is available. See Section 5.4.1.

###  5.5 Time complexity

Another metric to be taken into account is algorithm complexity, which
manifests itself in three interlinked aspects: time-, storage- and algorithm
simplicity. Like the classical principle Occam’s razor implies - there at all
times exists a preference for simpler models over more complex ones,
especially in the case both accomplish the same feat. So, an understandable
model is preferred that performs limited computational steps in order to
restrain time- and storage complexity from rising too high.

Although measuring algorithm time- and storage requirements does provide some
understanding, more insightful would be a theoretical description of the
complexity in Big O notation. Theoretical complexities are harder to obtain
though - leaving many authors to resort solely to measurements. Nonetheless,
complexity analysis is recommended to be part of any feature selection method
comparison, for it is a critical aspect to consider.

###  5.6 Statistical integrity

Given that the above discussed evaluation methods are computed, a statistical
test is ought to be applied to provide convincing evidence for any algorithm’s
superiority. Like explained in (Demšar, 2006), many papers make implicit
hypotheses acclaiming improved performance over existing methods. Although the
chosen metrics might have been appropriate to statistically show performance
gains as significant, the results are less reliable when left not validated by
a statistical test. Therefore, a statistical verification step is required in
the feature selection evaluation process.

Comparing multiple feature selectors is a non-trivial problem, which can be
seen as the problem of comparing multiple classifiers. To compare multiple
classifiers, (Demšar, 2006) recommends the Wilcoxon signed ranks test for
comparison of two classifiers and the Friedman test for comparison of more
classifiers given multiple datasets. Accompanying the Friedman test, it is
recommended to perform corresponding post-hoc tests, such as the Nemenyi test
(Nemenyi, 1963). The usually popular ANOVA test was advised against because,
given the context of machine learning, ANOVA assumptions are violated, e.g.
ANOVA assumes samples are drawn from normal distributions and the requirement
of random variables having equal variance.

A Wilcoxon signed ranks test and the Friedman test, are for these reasons
recommended for two- or more methods respectively, to statistically verify
significant differences in performance of feature selection algorithms over
the others using some summarizing scalar per feature selection method per
dataset. The Friedman test is recommended to be accompanied with a post-hoc
test, like Nemeyi’s post-hoc test.

##  6 Building a benchmarking pipeline

To facilitate conducting a comprehensive experiment, an accompanying software
pipeline is to be built. Such a pipeline must be able to run feature ranking
methods, select feature subsets, and validate the quality of the feature
subsets using some validation estimator. After both the feature ranking- and
validation step, evaluation metrics have to be computed - assessing the
performance of the feature ranking both directly- and indirectly. After the
metrics have been computed, results should be stored either on-disk or be
uploaded to the cloud, after which they should be made interpretable by
visualizing them. Furthermore, because the ‘direct’ evaluation of feature
rankings is only possible when a ground-truth reference of the relevant
features exists (Section 5.3), it is also desired to be able to generate
synthetic datasets built in to the framework. In this way, a practitioner
using the pipeline can easily design synthetic datasets suited for the purpose
at hand.

Challenges lie mainly in the scalability of the system. To conduct an
experiment of moderate size, a single machine can be used for all stages of
the pipeline: run the algorithms, collect and visualize the data. In such a
small-scale context, the data- collection and visualization stage could as
well be manual- no huge amounts of data have to be processed, so such a
collection- and visualization stage could be run from a single script. When a
larger experiment is desired to be run, however, a more sophisticated system
is desired. Practitioners must not be restricted in running the pipeline on a
single CPU or single machine, i.e., the system should support both
horizontal\- and vertical scaling. To design a system that can run experiments
in such scalable and distributed ways, the system must be modular to
facilitate distributing jobs over multiple machines, although it is also
necessary to aggregate results of individual jobs as part of a pipeline step
itself - requiring a piece of code to be aware of every job having finished.

An implementation of such a pipeline was made in Python, for the purposes of
this paper. Many of the challenges described above were addressed and overcome
by using an architecture that was built from the ground up to be scalable.
Both the data- collection and aggregation processes were automated - after
having run the pipeline for a collection of feature ranking methods and
datasets, results are automatically stored, aggregated, and visualized. In
this chapter, the reader is explained the general architecture of the system
and the components that make the system work - after which the usage of the
pipeline is shown along with its capabilities for scaling. The chapter is
meant to be explanatory at first and provide concrete instructions and
examples afterward.

###  6.1 Architecture

To first get an idea of the structure of the system, comments are made on its
architecture. The pipeline has two main components: an encapsulating pipeline
‘run’ script and the pipeline itself. The pipeline run script is to be called
the main script: allowing us to differentiate easily between the two. Whilst
the main script is designed to be generic and work for any ML related
benchmarking task, the pipeline itself is to be tailor made for the specific
task at hand. Ideally, one could swap out the pipeline for another (custom
built) one, but still using the main script mechanics: allowing a practitioner
to create various but related tasks using a single framework.

####  6.1.1 The main script

The sole purpose of the main script is to do three things: (1) load a dataset,
(2) split the dataset using a CV method and lastly (3) run the pipeline. In
this procedure, any dataset is loaded into the system using an adapter \-
allowing data to come in from any source by use of a plug-able architecture
(Section 6.2.1). Immediately after the dataset has been loaded a CV split is
conducted. In this way, because the cross-validation step is separated from
all the pipeline steps itself, it is less likely that a pipeline can go faulty
on the CV process - since the data splitting operation was already conducted
once a pipeline implementer gets their hands on the data. An illustration of
the main script, with as part of it the pipeline, can be seen in Figure 4.

Figure 4: A general architecture schematic of the built benchmarking
pipeline’s main script: providing data and instructions to the pipeline. Once
the designated dataset is loaded using a configured adapter, the data is fed
to the pipeline after a CV split. The pipeline is first fit using the training
data and then scored using the testing data, after which the results are
collected and stored.

An important facet facilitating the operations in the main script is the
config object (as seen at the top-left in Figure 4). The config object
contains all configuration and parameters relevant to all pipeline steps,
essentially storing all instructions required to perform a single pipeline
run. It is thereby important to store all configuration and not to leave any
parameters set at random or use temporary variables - the pipeline execution
is desired to be deterministic and therefore making its results reproducible.
In this way, experiments conducted using the pipeline are better suited for
scientific experiments - a situation in which a user desires to be able to
store exactly the state that produces a certain result and explain its
findings in minute detail.

To facilitate configuring the system in such a way, several options exist in
the Python ecosystem. One could use the built-in optparse module, or its
Python 3 counterpart argparse. The built-in packages leave much room for
improvement, however. The modules do not allow loading config from a file, do
not allow nesting, has only basic input validation and provides no tools for
more complex configuration scenarios like the pipeline at hand. A solution for
this is the Hydra library (Yadan, 2019), a framework providing powerful
utilities for configuring a Python application. The package allows one to
create hierarchical and composable configurations, to be set using either the
command-line or using yaml files. A simple yaml definition configuring the
main script to use the ‘Iris Flowers’ dataset and perform a K-fold CV can be
seen in Listing 1.

Listing 1: A simple configuration for the main script, using the ‘Iris
Flowers’ dataset and K-fold CV.

⬇

defaults:

- base_config

- base_rank_and_validate

- dataset: iris

- cv: kfold

- callbacks:

- wandb

- storage_provider: wandb

It can also be observed, in Listing 1, that there exist configuration
possibilities for ‘callbacks’ and a ‘storage provider’. These allow
customizing a back-end for sending the result data to and an adapter for
storing- or caching pipeline steps, respectively. Callbacks and the storage
provider will be elaborated upon in Section 6.2.6 and Section 6.2.5,
respectively. First, a look is taken at the pipeline implementation for
feature ranking and subset validation, which is called Rank and Validate.

####  6.1.2 The ‘Rank and Validate’ pipeline

With the main script now defined, a better look can be taken at the pipeline
implementation itself. Whereas the main script was still a general-purpose ML
benchmark runner, an implementation of the system by means of a pipeline is
task-specific. i.e., one such pipeline is meant to execute a specific ML task,
such as performing evaluation on a feature ranker. The pipeline does so in a
couple steps:

  1. 1.

Resample the dataset. Using the configured settings, the dataset is resampled.
Options that are ought to be supported are shuffling, and more importantly;
bootstrapping. It is thereby important to fixate a random ‘seed’, which allows
one to reproduce the resampling exactly.

  2. 2.

Rank features. Next, the feature ranker is set to work. The resampled data is
passed to the ranker’s fit function, after which the stored ranker might be
cached.

  3. 3.

Validate feature subset. Finally, after the ranker has been fit, its ranking
is used to validate its feature subsets, using some validation estimator.
Again, the fit validation estimator might be cached. The final results are
stored to disk and/or uploaded to the cloud.

This process is illustrated in Figure 5.

Figure 5: The ‘Rank and validate’ pipeline. The pipeline performs feature
ranking, subset selection and subset validation, in a sequence of loops, such
to facilitate bootstrapping. The pipeline is first fit and then scored.

As can be seen illustrated by rectangles in Figure 5, the above steps are
executed in a specific sequence of loops. First of all, all steps (1-3) are
passed through twice: once for the training data using the fit function and
once using the testing data using the score function. The reasoning behind the
division into these two steps is to force keeping application logic separate,
i.e., to enforce a separation of concerns. Because the fitting step is
essentially different than the scoring step, it makes sense to separate them.
This separation also allows better implementations of caching, because a
practitioner is enforced to store all estimator related state inside the
instance itself. This distinction is on par with the sci-kit learn (Pedregosa
et al., 2011) API.

For each of the two passes, the entire process is run $B$ separate times: once
for each bootstrap. The random resampling seed is changed with each bootstrap
iteration accordingly: making sure each bootstrap dataset is a different
random resampling of the dataset, but yet keep the pipeline to be
reproducible. The results from each bootstrap are then stored with each
bootstrap random state attached, allowing a practitioner to differentiate
between results and possibly aggregate- and group them.

Finally, each of the bootstrap iterations performs a number of feature subset
validations. Ideally, in the case where a dataset has $p$ dimensions and a
ranker has computed a feature importance score for each of them, it would be
desired to validate $p$ feature subsets - the first subset starting with only
the best ranked feature and subsequently including lesser-ranked features.
However, this would fast become computationally intractable: in the case where
$p$ is of considerable size, the amount of feature subsets to validate becomes
vast, fast. A compromise between the thoroughness and computational
tractability is to only validate some number of the best feature subsets up to
a certain dimension. The choice that is made in this paper is to validate at
most 50 feature subsets and possibly less if $p<50$. i.e., the $\min(p,50)$
best feature subsets are evaluated. This means for every
$k\in\\{1,\ldots,\min(p,50)\\}$ the $k$ best features are used for subset
validation.

The scoring process happens only after all estimators in the pipeline have
been fit. Once the fitting phase is complete, the testing sets will be passed
to the pipeline’s score function. Each estimator gets the opportunity to
generate some score given the testing dataset. Because both the rankers and
estimators are considered estimators in the system, a ranker will also get
such scoring opportunity. More specifically so, rankers and validators are
scored in the following ways:

  * •

Rankers are scored using the dataset ground-truth relevant features, if
available. Thereby metrics as described in Section 5.3 are used: the R2 score
and the logistic loss.

  * •

Validators are scored like described in Section 5.2. i.e., the scoring depends
on the dataset task: the R2 score is used in the case of regression, and the
accuracy score is used in case of classification.

That said, a general picture of the pipeline architecture and data flow is now
obtained. The pipeline is, however, more complex. The pipeline supports
multiprocessing, distributed computing, caching, and data visualization.
Moreover, each pipeline component has very different functionalities and can
be configured separately from the rest of the pipeline. Therefore, a look is
taken into the separate pipeline components first in Section 6.2, to then take
an in-depth view of the pipeline execution afterwards, in Section 6.3.
Finally, an elaboration is made on data visualization in Section 6.4.

###  6.2 Components

Each pipeline component can be configured using their own yaml files and has a
reasonable level of isolation with respect to its module structure. The main
components are discussed briefly.

####  6.2.1 Datasets

The system allows using datasets from numerous sources, using adapters.
Adapters are designed to be generic interfaces for fetching data from any
source, as long as it can be expressed in two matrices $\mathbf{X}$ and
$\mathbf{y}$. Although some adapters are available built-in to the framework,
new adapters can easily be configured. Any such adapter must only implement a
get_data function, which returns the two matrices $\mathbf{X}$ and
$\mathbf{y}$, given a configuration object facilitated by Hydra (see Section
6.1.1). Such adapters can implement any logic that retrieves the data
themselves, be it logic for loading a dataset from local disk or from an
internet platform.

OpenML is one such platform. OpenML (Vanschoren et al., 2014) is a platform
for sharing and organizing data, completely open to all. Due to the flexible
adapter structure of the system, and thanks to the availability of a Python
package providing interfacing support with the platform, an integration with
OpenML was made possible. Using the interface, datasets are fetched from a
remote and cached to disk once downloaded. In this way, a cached version of
the dataset is used in subsequent requests - making sure a dataset is not
unnecessarily downloaded again. Integration with the platform allows users of
the pipeline to access a large library of datasets, varying in type and
domain. The dataset library features both regression and classification
datasets, and of various target types, i.e., both univariate and multivariate.
Also, the platform has multiple existing benchmark suites that can be used,
such as the OpenML-CC18 (Bischl et al., 2019). An example definition of a
configuration file for loading an OpenML dataset is to be found in Listing 2.

Listing 2: A dataset config for loading the ‘Iris Flowers’ dataset from
OpenML.

⬇

name: Iris Flowers

task: classification

adapter:

_target_: fseval.adapters.OpenML

dataset_id: 61

target_column: class

It can be seen, in Listing 2, that a ‘target’ class can be defined. This is
passed to Hydra, which in turn instantiates the class with the given
configuration. This is one of the core features that makes possible the
modular adapter architecture that is in place.

Synthetic datasets can also be generated using the modular adapter-structure
in place. By viewing a synthetic data generator as yet another adapter,
providing data given a configuration object, both real-world and synthetic
datasets can be used and generated using a single interface. A straight-
forward way to generate synthetic datasets is to use the sci-kit learn
datasets module, which provides functions for drawing samples from a wide
variety of distributions. Aside from generic generators, there also exists
support for drawing from more specific distributions, such as ‘two
interleaving half circles’, an ‘S curve’ or a ‘swiss roll’ distribution.

Most interesting for the current purposes, however, are two functions for
generating classification- and regression datasets, make_classification and
make_regression, respectively. Besides allowing one to configure the amount of
samples and dimensions to generate, one can also define the desired
informative-, redundant-, repeated- and irrelevant features to generate. In
this way, a ground-truth for the relevant features can be constructed. In the
case of regression, one can also retrieve the coefficients, or weighting, for
each feature which are to be approximated, i.e., the desired feature
importance for each feature. This means that whilst in the case of
classification one has access to a binary ground-truth, an exact measure of
the desired feature importance can be retrieved in the case of regression. An
example definition of a synthetic classification dataset can be found in
Listing 3.

Listing 3: A dataset config generating a synthetic dataset using the sci-kit
learn make_classification function.

⬇

name: Synclf hard

task: classification

adapter:

_target_: sklearn.datasets.make_classification

class_sep: 0.8

n_classes: 3

n_clusters_per_class: 3

n_features: 50

n_informative: 4

n_redundant: 0

n_repeated: 0

n_samples: 10000

random_state: 0

shuffle: false

feature_importances:

X[:, 0:4]: 1.0

Worth noting in Listing 3 is the presence of a ‘feature_importances’
attribute. This is the exact definition of the ground-truth relevant features,
used by the pipeline. The system allows one to define the feature importances
as numpy-indexed selectors, operating on some matrix the same size as
$\mathbf{X}$. Every instance and every accompanying dimension can have a
specific ground-truth weighting, although for the current purposes only global
feature- ranking and selection are considered. In the example above, it can be
seen that the first 5 features are defined to be equally relevant: and
therefore all have a uniform weighting of 1.0.

####  6.2.2 Cross-Validation

The Cross-Validator can be configured to be any class that has a split
function, taking in the matrices $\mathbf{X}$ and $\mathbf{y}$ and returning a
generator. The generator must then be an iterator over the number of
configured folds, with each split returning the designated training- and
testing data. An example configuration for a 5-Fold CV split can be seen in
Listing 4.

Listing 4: A config for 5-Fold CV, with shuffling. The split is reproducible
due to the fixed random seed.

⬇

name: K-Fold

splitter:

_target_: sklearn.model_selection.KFold

n_splits: 5

shuffle: True

random_state: 0

fold: 0

As is seen in Listing 4, a CV technique can be configured to be a sci-kit
learn module. Another important thing to note is the fixed ‘fold’ attribute,
set in this case to 0. Because the pipeline always executes exactly one fold,
performing a K-Fold CV is done by executing the pipeline multiple times using
different values of fold.

####  6.2.3 Resampling

Built in, there are two options for resampling. There exist options for (1)
bootstrap sampling and (2) shuffle resampling. Whilst the latter only changes
the permutation of the dataset, the former performs resampling with
replacement, therefore actually changing the dataset distribution. Bootstrap
resampling is configured like as can be seen in Listing 5.

Listing 5: A config for bootstrap sampling using a sample size of 1.0, i.e,
using the full dataset.

⬇

name: Bootstrap

replace: true

sample_size: 1.00

Noteworthy, in Listing 5, is the ‘sample size’ parameter. If needed, the
dataset can be down- or up-sampled. Although up-sampling generally does not
make much sense, down-sampling can increase the stochasticity of the
bootstrapped dataset permutation. This means that, if desired, one can create
more ‘random’ distributions by setting sample_size to a lower number. Lastly,
a random seed can be fixed, like as can be seen in the pipeline schematic,
Figure 5: the resampling random_state is set to the bootstrap number. In this
way, the resampling is exactly reproducible.

####  6.2.4 Estimators

In the context of the current system, an ‘estimator’ can mean two things. An
estimator is either: (1) a feature ranker or (2) a validation estimator.
Still, the two are combined into a single estimator interface: which is
because both can share much of the same API for interacting with them. Both
implement a fitting- and scoring step, and can be of regressor- or classifier
type. Also, to support these two learning tasks in a better way, any estimator
is encapsulated in a generalized estimator interface. An example configuration
for a ranking estimator can be seen in Listing 6.

Listing 6: An estimator config for ReliefF (Kononenko, 1994).

⬇

name: ReliefF

classifier:

estimator:

_target_: skrebate.ReliefF

regressor:

estimator:

_target_: skrebate.ReliefF

estimates_feature_importances: true

estimates_feature_support: false

Both a regressor and classifier were defined in a single interface, as can be
observed in Listing 6. Although the two were in this case defined to be the
same, they might have as well pointed to different modules. Under the hood,
the system detects the configured dataset task at initialization and
instantiates either one of the two modules. Furthermore, the system has to be
told the capabilities of the estimator, e.g., whether it estimates targets,
feature importance or feature support. This is done using the boolean
attributes starting with ‘estimates_’: where in this case the ReliefF
estimator estimates just the feature importance, but no feature subset.

Validation estimators are configured in much the same way - and can also have
both a regressor and classifier defined. The interface also allows indicating
an estimator’s support for multioutput learning or the estimator requiring a
strictly positive dataset to work.

####  6.2.5 Storage providers

A storage provider is defined to be a bridge between the pipeline and a file
system, independent of whether this file system is on the local disk or in the
cloud. Because in the pipeline it is desired that the fit estimators can be
cached, a storage system must be built accordingly. Independent of exactly
where the files are stored, a storage provider can be constructed by
implementing routines for saving and restoring files. As long as the storage
provider returns a usable file object given a file path, the system is
ambivalent about how the file is loaded. A schematic of how this storage
provider works is to be seen in Figure 6.

Figure 6: An illustration of the Storage Provider functionality. The storage
provider can be used to save- and restore files, using both a local- or remote
file system.

In the example in Figure 6, it can be seen that the storage provider is used
by an estimator to save- and restore data. Such data can, for example, be a
pickle file used to serialize- and deserialize a fit estimator object.

Caching is supported in this way. By storing a reference to the file-path
right in the job configuration itself, cached files can be reused via either
the local- or remote file system. The pipeline allows configuring exactly
which parts of the pipeline allow using cached versions of their fit
estimators, such that in any subsequent run of a job, some estimators can be
reused and others overwritten. This functionality facilitates running multiple
validation estimators given a single fit feature ranker, for example.

####  6.2.6 Callbacks

The final component to discuss is the Callback. A callback can be used to send
data to a back end during the execution of the pipeline. Such data can be
configurations, objects containing metrics or entire tables with results. The
requirements for implementing such a callback are low, since any such data can
be stored as an opt-in. If a pipeline user desires so, the usage of a callback
can be completely ignored, storing only data to the disk using the storage
provider. Sending metrics to a callback can, however, be a powerful way to
send data to a database or cloud back-end right in the pipeline. In this way,
the usual data collection phase can be automatized. About an implementation of
such a callback can be read in Section 6.4. First, however, a closer look is
taken at the pipeline execution and its scalability.

###  6.3 Execution

The pipeline execution happens in a predefined number of steps. Because the
pipeline is desired to be run in parallel and in a distributed way, the
execution of such a system is nontrivial. Support for scaling in both the
horizontal- and vertical directions are required, i.e., the pipeline must be
able to scale over more machines but must also scale over more system
resources such as CPU or RAM. Systems that fulfil such desires are High-
Performance Computing (HPC) systems, which are designed to support scalable
computing (Ristov et al., 2016). An elaboration is made on an implementation
of such a scalable architecture.

####  6.3.1 Multiprocessing

To facilitate vertical scaling, a multiprocessing implementation is built. It
is thereby the goal to best utilize the system resources at hand: maximizing
the usage of the available processors. To do so, certain iterative parts of
the pipeline were wrapped in an encapsulating ‘Experiment’ module. The task of
the module is, then, to pick up wherever the pipeline would iterate over a
number of steps, and add multiprocessing support. As could be read in Section
6.1.2 and seen in Figure 5, the pipeline has two main loops: the first (1)
runs $B$ bootstraps and the second (2) validates $\min(p,50)$ feature subsets.
In this implementation, the choice was made to distribute the former loop,
i.e., to distribute the bootstraps over multiple CPU’s. An illustration making
clear the functionality of such a multiprocessing approach can be seen in
Figure 7.

Figure 7: The pipeline multiprocessing distribution process. In the pipeline
fit step the set of bootstraps is distributed over the available CPU’s.

As could be seen in Figure 7, every bootstrap resampling job is assigned to a
processor. In case more bootstraps are to be run than the number of available
processors, i.e., $B>$ # CPU’s, some CPU’s are assigned more than one
bootstrap resampling job. In HPC environments, it is common to have a dozen
CPU’s available in a single job. Given that a reasonable but still feasible
number of bootstraps to run is also a few dozen, it can be straight-forward
option to configure the number of bootstraps to run as the number of CPU’s
available in a job on the designated HPC environment.

Important to note in such a multiprocessing job, is the available scope of
variables accessible inside any subprocess. Because each subprocess runs
isolated on a CPU, it cannot access variables exclusively available to the
main thread process. In the case that such a subprocess requires access to
variables from the main thread, either one of two options can be considered.
One option is to pass down a copy of the variables in question to the
subprocess: this, however, does not apply when such variables change during
the program execution. In that case, the second option is to communicate
between the subprocess and the main thread. In this way, functionality such as
storing cache files (Section 6.2.5) can still be done in the main thread, if
it is deemed necessary.

####  6.3.2 Distributed computing

To make the pipeline horizontally scalable, an approach must be constructed
for distributing the processes over multiple machines. Although running
multiple jobs in a HPC environment can be a relatively straight-forward
practice, a choice must be made on the amount of processing to put in a single
job. For example, if a system were in place that automatically distributes
jobs over a number of processors and machines by itself, the horizontal- and
vertical scaling approaches could be tackled in a single solution. In this
scenario, an automatic job scheduler would assign very small jobs to nodes and
processors by itself; the smallest experimental unit could be running one
bootstrap on a single feature ranker. However, such an approach is rather
complex and often unsupported by existing HPC systems. Therefore, the choice
is to go with a hybrid approach - a single HPC job is assigned multiple
experiments. In this way, long queue waiting times as a consequence of running
many small jobs are prevented.

Such an approach can be accomplished by the use of a second job queue, next to
the HPC job queue. The queue of choice is RQ, a library for keeping track of
an arbitrary number of job queues in Redis which then can be executed by
workers. Workers are run on the HPC system itself and are built solely to
execute work coming from the queues. The library is built to be fault-
tolerant, scalable, and easy to use. The basic idea is that during the
enqueueing phase, a piece of Python code is serialized into a binary format
and stored in Redis into a First-In-First-Out queue. Then, once the time has
come that the job can be executed, the job execution code is deserialized from
Redis back into Python code, and executed. This functionality can be seen
illustrated in Figure 8.

Figure 8: The RQ enqueueing- and dequeueing process. Jobs are serialized on
the client and stored in Redis FIFO queue. Workers on the HPC system process
the jobs on the designated queues.

It can be seen in Figure 8, that the Redis database acts as an intermediary
between the user and the HPC system. Whilst the queues in Redis may contain
hundreds or thousands of jobs, the amount of jobs on the HPC system only need
to be several: one job for each worker. Each worker may run on its own node,
or share a node by means of a virtual machine. In overall, the RQ-based
approach has a number of benefits.

  * •

First of all, the HPC system is not overloaded by small jobs. Instead, RQ is
better designed to store and process a multitude of small jobs - preventing
overhead and long queuing times in the HPC system. This is because the only
jobs that are launched on the HPC system are the workers - which can be
configured to last an arbitrary amount of time, or to terminate once its queue
is empty.

  * •

A second benefit is that RQ has good support for handling job failures. RQ
allows automatically retrying failed jobs, with optionally set time intervals
between tries, such that in the case of a network time-out or system error the
job is be restarted without further human input. Furthermore, if the job still
fails execution after its pre-configured number of retries were completed, the
job is not just removed from all queues. Instead, the job is moved to a
separate ‘failed job’ queue. The error logs for jobs in this queue are easily
inspected using a dashboard interface for RQ.

  * •

A last benefit is the serialization- and deserialization process in RQ.
Because Redis does not just store instructions to execute the code, but the
actual execution code itself, code on the HPC system might be changed during
the time a job is in queue. This gives a practitioner more flexibility to
change program code and run the pipeline at the same time.

As a side benefit, the Hydra library used for configuring the system has
built-in support for working with RQ. By use of a plugin, the pipeline can be
configured to enqueue its jobs to RQ instead of executing right away. An
important note is that due to the serialization process the Python version
during enqueueing and dequeueing have to be exactly equal.

###  6.4 Data- collection and visualization

Finally, all results coming out of the pipeline should be collected and
visualized. Since the pipeline has support for plugging in any arbitrary
‘Callback’ (Section 6.2.6), primitives are in place for sending the data to a
data- collection and visualization back-end. Any such callback can implement
functions for processing the pipeline configuration, individual metric
objects, and result tables. Because the pipeline is built to be executed at
scale, a clear distinction is made between the pipeline execution phase and
the data collection phase: the two are separated by means of the Callback
interface.

One callback implementation is the ‘Wandb’ callback, providing integration
with the Weights and Biases platform (Biewald, 2020). Weights and Biases is a
platform for real-time experiment tracking and visualization, with built-in
primitives in Python. The platform allows you to upload tables, metrics, or
configuration objects, and subsequently visualize them in the front-end
interface. The platform also allows one to upload- and download raw files,
allowing a user to integrate remote caching functionality (Section 6.2.5). A
screenshot of the platform front-end can be seen in Figure 9.

Figure 9: The Weights and Biases dashboard front-end. The platform allows
tracking experiments, aggregating data, and visualizing results.

An important feature on the Weights and Biases platform is the ability to
collect- and aggregate data interactively in the front-end. This is done by
means of ‘Vega-Lite’ (Satyanarayan et al., 2017), a visualization grammar for
creating interactive visualizations in a front-end environment. This allows a
user to send ‘raw’ data to the dashboard, and aggregate it on the front-end
itself. Any aggregation operation can be applied, be it taking averages, sums
or computing the variances. The front-end then allows you to visualize the
aggregated data using a simple grammar. In this way, the pipeline requires
much less small adjustments and modifications for the user to be satisfied
with the data aggregation process - the data can be summarized in a later
stage after all results are already in.

##  7 Experiments

Accompanying the proposal of a new feature- ranking and selection methodology,
there is an experiment. The aim of the experiment is to show an example of
what is possible with both the pipeline, and the newly proposed evaluation
methodology. In this way, ML practitioners and authors of new methods are able
to conduct experiments themselves with the same setup and configuration:
therefore allowing comparison between the experiments.

###  7.1 Experiment setup

The experiment closely follows the outlines of Section 5 and Section 6. The
setup can be summarized like so:

  * •

The dataset at hand is split into a 80% training- and 20% testing dataset.
Whereas the training set is reserved for fitting the feature ranker and
validator, the testing set is reserved for scoring the validator.

  * •

For each Feature Ranker and dataset, $B=25$ bootstrap resampling iterations
are run.

  * •

In each bootstrap resampling iteration, $\min(p,50)$ feature subsets are
evaluated, with each including the $k$ best features.

  * •

Validation estimators are evaluated with the R2 score in case of regression
and with accuracy in the case of classification. The validation estimators
used are k-NN with $k=5$ and a DT at default sklearn settings.

  * •

In the aggregation process, in case several experimental runs are found with
the same configuration, i.e., for one Feature Ranker executed on one dataset,
only the ‘best’ run is considered. The best run is considered to be the run
with the highest average mean score over all features, considering all
bootstraps. This aggregation step will be elaborated upon in Section 7.2.

All experiments were run on a SLURM HPC environment. Specifically, the
University of Groningen has a ‘Peregrine’ compute cluster, with machines of
various types, of which the most common one is a 24 core machine powered by
two Intel Xeon E5 2680v3 CPUs running at 2.5 GHz. Per node, 128 GB memory and
1 TB internal disk space is available, but 10 GB was requested per CPU
instead. This accounts for a total of 100 GB memory for 10 CPU processes. In
total, the experimentation facilitating these results took 7682 hours of
processing time on the above-mentioned machines.

A line-up of feature rankers and datasets was constructed to conduct
benchmarking on. The full list of both is available in the Appendix Section A.
See Table 3 for the Feature Ranking line-up, and Table 4 for the datasets
line-up. The distribution of classes can be found in Appendix Figure 24. Most
datasets are balanced, though also some imbalanced ones are included.
Imbalanced datasets are not treated specially, it is up to the ranking- and
validation estimators to correct for this. It can also be seen that rankers of
various types are included: both regressors and classifiers, and rankers that
support both learning tasks. Furthermore, some rankers were included that
support multioutput datasets, i.e., datasets with multivariate targets.

###  7.2 Experimental results for the ‘Synclf hard’ dataset

To best understand the format of the experimental results and the accompanying
metrics, a look is taken at the results for one dataset. In this way, a better
understanding of the charts is gained first.

First of all, the considered dataset is the ‘Synclf hard’ dataset. It is a
synthetically generated dataset, created using the sklearn make_classification
function (Section 6.2.1). Its full configuration specification is defined in
Listing 3. The dataset is defined to have $n=10000$ samples and $p=50$
dimensions. After the CV step was conducted, 8000 samples are left for
training. The dataset has 4 relevant features and 3 target classes, which are
perfectly balanced. That said, observations are first made on running one
feature ranker on the dataset. For this, ReliefF (Kononenko, 1994) is chosen.

####  7.2.1 ReliefF performance

To start, a look is taken at a plot that explicitly plots all bootstraps. See
Figure 10.

Figure 10: Validation performance for ReliefF on the ‘Synclf hard’ dataset,
for all 25 bootstraps.

As can be seen in Figure 10, the various bootstraps have had an effect on the
validation performance per subset. Due to the random resampling with
replacement taking place in the bootstrap phase, the feature ranker has to
deal with permutations of the dataset each run. Indeed, this randomness is