# A Taxonomy of Error Sources in HPC I/O Machine Learning Models

Mihailo Isakov1, Mikaela Currier1, Eliakin del Rosario1, Sandeep Madireddy2,
Prasanna Balaprakash2, Philip Carns2, Robert B. Ross2, Glenn K. Lockwood3,
Michel A. Kinsy1  1 Secure, Trusted, and Assured Microelectronics (STAM)
Center
Ira A. Fulton Schools of Engineering, Arizona State University, Tempe, AZ
85281
<EMAIL_ADDRESS>2 Argonne National Laboratory, Lemont, IL 60439
<EMAIL_ADDRESS><EMAIL_ADDRESS>3 Lawrence Berkeley
National Laboratory, Berkeley, CA 94720
<EMAIL_ADDRESS>
###### Abstract

I/O efficiency is crucial to productivity in scientific computing, but the
increasing complexity of the system and the applications makes it difficult
for practitioners to understand and optimize I/O behavior at scale. Data-
driven machine learning-based I/O throughput models offer a solution: they can
be used to identify bottlenecks, automate I/O tuning, or optimize job
scheduling with minimal human intervention. Unfortunately, current state-of-
the-art I/O models are not robust enough for production use and underperform
after being deployed.

We analyze multiple years of application, scheduler, and storage system logs
on two leadership-class HPC platforms to understand why I/O models
underperform in practice. We propose a taxonomy consisting of five categories
of I/O modeling errors: poor application and system modeling, inadequate
dataset coverage, I/O contention, and I/O noise. We develop litmus tests to
quantify each category, allowing researchers to narrow down failure modes,
enhance I/O throughput models, and improve future generations of HPC logging
and analysis tools.

###### Index Terms:

High performance computing, I/O, storage, machine learning

##  I Introduction

As scientific applications push to leverage ever more capable computational
platforms, there is a critical need to identify and address bottlenecks of all
types. Due to the large data movements in these applications, the I/O
subsystem is often a major source of performance bottlenecks, and it is common
for applications to attain only a small fraction of the peak I/O rates [1].
These performance problems can severely limit the scalability of applications
and are difficult to detect, diagnose, and fix. Data-driven machine learning-
based models of I/O throughput can help practitioners understand application
bottlenecks (e.g., [2, 3, 4, 5, 6, 7]), and have the potential to automate I/O
tuning and other tasks. However, current machine learning-based I/O models are
not robust enough for production use [6]. A thorough investigation of _why_
these models underperform when deployed on high performance computing (HPC)
systems will provide key insights and guidance on how to address their
shortcomings. The goal of our study is to help machine learning (ML)-driven
I/O modeling techniques make the transition from theory to practice.

There are several reasons why machine learning-based I/O models underperform
when deployed: poor modeling choices [2, 7], concept drift in the data [5],
and weak generalization [6], among others. I/O models are often opaque, and
there is no established methodology for diagnosing the root cause of model
errors. In this work, we present a taxonomy of ML-based I/O modeling errors,
as shown in Figure 1. Through this taxonomy, we show that I/O throughput
prediction errors can be separated and quantified into five error classes:
inadequate (1) application and (2) system models, (3) novel application or
system behaviors, (4) I/O contention and (5) inherent noise. For each class,
we present data-driven litmus tests that estimate the portion of modeling
error caused by that class. The taxonomy enables independent study of each
source of error and prescribes appropriate ML techniques to tackle the
underlying sources of error.

Our contributions in this work are as follows:

  1. 1.

We introduce a taxonomy of ML-based I/O throughput modeling errors which
consists of five causes.

  2. 2.

We show that choice of machine learning model type, scaling model size, and
hyperparameter tuning cannot reduce all potential errors. We present two
litmus tests that quantify error due to poor application and system modeling.

  3. 3.

We present a litmus test that estimates what portion of error is caused by
rare jobs with previously unseen behavior, and apply uncertainty
quantification methods to classify those jobs as out-of-distribution jobs.

  4. 4.

We present a method for quantifying the impact of I/O contention and noise on
I/O throughput, which (1) defines a fundamental limit in how accurate ML
models can become, and (2) gives HPC system users and administrators a
practical estimate of the I/O throughput variance they should expect. We show
that underlying system noise is the dominant source of errors, and not poor
modeling or lack of application or system data.

  5. 5.

We present a framework for how the taxonomy is practically applied to new
systems and evaluate it on two leadership-class supercomputers: Argonne
Leadership Computing Facility (ALCF) Theta and National Energy Research
Scientific Computing Center (NERSC) Cori.

##  II Related work

In recent years, automating HPC I/O system analysis through ML has received
significant attention, with two prominent directions: (1) workload clustering
to better understand groups of HPC jobs and automate handling of whole groups,
and (2) I/O subsystem modeling and make predictions of HPC job I/O time, I/O
throughput, optimal scheduling, etc. Clustering HPC job logs has been explored
in [2, 8, 9] with the goal of better understanding workload distribution,
scaling I/O expert effort more efficiently, and revealing hidden trends and
I/O workload patterns. ML-based modeling has been used for predicting I/O time
[4], I/O throughput [7, 2], optimal filesystem configuration [10, 11], as well
as for building black boxes of I/O subsystems in order to apply ML model
interpretation techniques [2]. While there have been some attempts at creating
analytical models of I/O subsystems [12], most attempts are data-driven, and
rely on HPC system logs to create models of I/O [1, 7, 4, 2, 13]. Although the
challenges of developing accurate machine learning models are well known, the
nature of the domain requires special consideration: I/O subsystems have to
service multiple competing jobs, their configuration evolves over time, they
have periods of increased variability, they experience occasional hardware
faults, etc. [14, 15, 16]. Diagnosing this I/O variability, where the
performance of a job depends on external factors to the job itself has been
extensively studied [3, 14, 17, 4, 16]. Finally, the deployment of I/O models
has been shown to require special consideration as these models often
significantly underperform on new applications [5, 6]. While different sources
of model error have been studied individually, no prior work characterizes the
relative impact of different sources of error on model accuracy.

##  III Modeling HPC applications and systems

The behavior of an HPC system is governed by both complex rules and inherent
noise. By formalizing the system as a mathematical function (or, more
generally, a stochastic process) with its inputs and outputs, the process may
be decomposed into smaller components more amenable to analysis. The I/O
throughput of a system running specific sets of applications may be treated as
a data-generating process from which I/O throughput measurements are drawn.
While building a perfect model of an HPC system may not be possible, it is
useful to understand the inputs to the ‘true’ process and the process’s
functional properties. The theoretical model of the process must include all
causes that might affect a real HPC system, such as: how well a job uses the
system, hardware and software configurations over the life of the system,
resource contention between concurrent jobs, inherent application-specific and
system noise, as well as application-specific noise sensitivity. Although many
of these aspects are not observable in practice, the true data-generating
process takes these aspects into account.

We adapt the global modeling formulation from Madireddy et al. [4], and
formulate I/O throughput $\phi(j)$ of a job $j$ as:

$\phi(j)=f(j,\zeta,\omega)$ (1)

Here, $\zeta$ represents system state (e.g., filesystem health, system
configuration, node availability, etc.) and system behavior (e.g., the
behavior of other applications co-located with the modeled application during
its run, contention from resource sharing, etc.) at a given time. $\omega$
represents randomness acting on the system. The system $\zeta$ can be further
decomposed as:

$\zeta=\zeta_{g}(t)+\zeta_{l}(t,j)$ (2)

The component $\zeta_{g}(t)$ represents the global system impact on all jobs
running on the system (e.g., a service degradation that equally impacts all
jobs) and is only a function of time $t$. The component $\zeta_{l}(t,j)$
represents the local system impact on the I/O throughput of job $j$ caused by
resource contention and interactions with other jobs running on the system.
Contrary to the $\zeta_{g}(t)$ component, $\zeta_{l}(t,j)$ is job-specific and
depends on the behavior of the current set of applications running on the
system, the sensitivity of $j$ to resource contention and noise, etc. Without
loss of generality, the I/O throughput of a job from Equation 1 can be
represented as:

$\begin{split}\phi(j)&=f(j,\;\zeta_{g}(t),\;\zeta_{l}(t,j),\;\omega)\\\
&=f_{a}(j)+f_{g}(j,\zeta_{g}(t))+f_{l}(j,\zeta_{l}(t,j))+f_{n}(j,\zeta,\omega)\end{split}$
(3)

Here, $f_{a}(j)$ represents job throughput on an idealized system where the
job is alone on the system, the system does not change over time, and there is
no resource contention. $f_{g}(j,\zeta_{g}(t))$ represents how the evolving
configuration of the system (hardware provisioning, software updates, etc.)
affects a job’s I/O throughput. The $f_{l}(j,\zeta_{l}(t,j))$ component
represents the impact of resource contention and the sensitivity of $j$ to I/O
noise. Finally, $f_{n}(j,\zeta,\omega)$ represents the impact of inherent
system noise (e.g., packets dropped) on the job.

###  III-A Modeling assumptions

The task of modeling a system’s I/O throughput is to predict the behavior of
the system when executing a job from some application on some data. Modeling
I/O throughput requires modeling both the HPC system and the jobs running on
it. Machine learning models used in this work attempt to learn the true
function $\phi$ by mapping observable features of the job $j$ and the system
$\zeta$ to measured I/O throughputs $\phi(j)$. A model $m(j_{o},\zeta_{o})$ is
tasked with predicting throughput $\phi(j)$, where $j_{o}$ and $\zeta_{o}$ are
the observable job and system features.

When designing ML models, the choice of model architecture and model inputs
has implicit assumptions about the process that generates the data. When
incorrect assumptions are made about the domain, the model will suffer errors
that cannot be fixed within that modeling framework. We investigate four
common assumptions about the HPC domain.

All data is in-distribution: a common assumption that practitioners make is
that all model errors are the product of insufficiently trained models,
inadequate model architectures, or missing discriminative features. However,
some jobs in the dataset may be Out of Distribution (OoD), that is, they may
be collected at a different time or environment, or through a different
process. The model may underperform on OoD jobs due to the lack of similar
jobs in the training set and not due to lack of insight (features) into the
job. The cause of the problem is epistemic uncertainty (EU) \- the model
suffers from lack of knowledge or reducible uncertainty, since a broader
training set would make the OoD jobs in-distribution (ID). In the HPC domain,
epistemic uncertainty is present in cases of rarely run or novel jobs or
uncommon system states. Without considering the possibility that a portion of
the error is a product of epistemic uncertainty, practitioners will put effort
into tuning models instead of collecting more underrepresented jobs. Referring
to Equation 1, this assumption may be expressed as: deployment time $j_{d}$
and $\zeta_{d}$ are drawn from a different distribution from training time
$j_{t}$ and $\zeta_{t}$.

Noise is absent: all systems have some inherent noise that cannot be modeled
and will impact predictions. Aleatory uncertainty (AU) refers to inherent
uncertainty, either due to noise or missing features about jobs in the
dataset. Modeling errors due to aleatory uncertainty are different from
epistemic uncertainty because collecting more jobs may not reduce AU. HPC I/O
domain experts note that certain systems do have significantly higher or lower
I/O noise [18, 19], but quantifying noise in modern leadership computing
systems is an under-explored subject. Understanding and characterizing
system’s inherent I/O noise is key to account for the aleatory uncertainty in
the ML-based I/O models and better quantify the effect of this uncertainty on
the I/O throughput predictions. I/O modeling works rarely attempt to quantify
ML model uncertainty [20] even though an estimate of AU significantly helps in
model selection. The assumption that noise is not present in the dataset can
be expressed as follows: The practitioner assumes that the process has the
form of $\phi(j)=f(j,\zeta)$ instead of $\phi(j)=f(j,\zeta,\omega)$.

Sampling is independent: running a job on a system can be viewed as sampling
the combination of system state and application behavior and measuring I/O
throughput. Most I/O modeling works implicitly assume that multiple samples
taken at the same time are independent of each other. The system is modeled as
equally affecting all jobs running on it, that is, the placement of different
jobs on nodes, the interactions between neighboring jobs, network contention,
etc. do not affect the job. This assumption can then be expressed as: the
process has the form of $\phi(j)=f(j,\zeta_{g}(t),\omega)$, not
$\phi(j)=f(j,\zeta,\omega)$.

Process is stationary: a common assumption is that the data-generating process
is stationary, and that the same job ran at different times achieves the same
I/O throughput. As hardware fails, as new nodes are provisioned, and shared
libraries get updates, the system evolves over time. This assumption is
therefore incorrect, and ignoring it, e.g., by not exposing when a job is ran
to the ML model may cause hard-to-diagnose errors. In other words, the
assumption is that the process has the form of $\phi(j)=f(j,\omega)$ and
$f_{g}(j,\zeta_{g}(t))=0$.

##  IV Classifying I/O throughput prediction errors

No matter the problem to which machine learning is applied, a systematic
characterization of the sources of errors is crucial to improve model
accuracy. While there is no substitute for ‘looking at the data’ to understand
the root cause of the problem, this approach does not scale for large
datasets. We seek a systematic way to understand the barriers to greater
accuracy and improve ML models applied to system data.

The key questions we ask in this work are: What are the impediments to the
successful application of learning algorithms in understanding I/O? Should ML
practitioners focus on acquiring more data on HPC applications or the HPC
system? How much of the error stems from poor ML model architectures? How much
of the error can be attributed to the dynamic nature of the system and the
interactions between concurrent jobs? What fraction of jobs exhibit a truly
novel I/O behavior compared to the jobs observed thus far? At what point are
the applications too novel, so much so that users should no longer trust the
predictions of the I/O model? We now describe five classes of errors and dive
deeper into error attribution in Sections VI, VII, IX and VIII.

The lack of application and system observability, the inherent noise $\omega$,
and the OoD jobs prevent ML models from fully capturing system behavior,
causeing errors. We define the I/O throughput prediction error of a model $m$
in a job $j$ as:

$e(j)=\phi(j,\zeta,\omega)-m(j_{o},\zeta_{o})$ (4)

Following the $\phi(j)$ terms from Eq. 3 and including the out-of-distribution
error, the error can be broken down as follows:

$e(j)=e_{app}+e_{system}+e_{OoD}+e_{contention}+e_{noise}$ (5)

Here, the application modeling error $e_{app}$ is caused by a poor model fit
of application behavior ($f_{a}(j)$ component), the global system error
$e_{system}$ is caused by poor predictions of system dynamics
($f_{g}(j,\zeta_{g}(t))$ component), the out-of-distribution error $e_{OoD}$
is caused by weak model generalization on novel applications or system states,
the contention error $e_{contention}$ is caused by poor predictions of job
interactions ($f_{l}(j,\zeta_{l}(t,j))$ component), and the noise error
$e_{noise}$ is caused by the inability of any model to predict inherent noise
($f_{4}(j,\zeta,\omega)$ component). This leaves us with five classes of
errors shown at the bottom of Figure 1. While attributing such an error on a
per-job basis is difficult, we will show that estimating each component across
a whole dataset is possible.

Figure 1: Taxonomy of I/O throughput modeling errors, with examples of the
effects of each class of error - (a) Median errors of XGBoost models with
varying numbers of estimators and estimator depth.(b) I/O throughput
prediction error for sets of identical (duplicate) jobs, for 5 different
applications. (c) Distribution of I/O throughput prediction errors and time
ranges for sets of duplicate jobs. (d) Median errors before (green) and after
training (red), with error distribution shown in the back.

###  IV-A I/O Model Error Taxonomy and Litmus Tests

Errors in Equation 5 must be estimated in a specific order shown in Figure 1
due to the specifics of individual litmus tests. For example, before the
effect of aleatory and epistemic uncertainty can be separated, a good model
must be found [21]. Similarly, before global and local system modeling errors
can be separated, OoD jobs must be identified.

Application modeling errors:  ML models can have varying expressivity and may
not always have the correct structure or enough parameters to fit the
available data. Models whose structure or training prevents them from learning
the shape of the data-generating process are said to suffer from approximation
errors, which are further divided into application and system modeling errors.
Application modeling errors are caused by poor predictions of application
behavior which can be fixed through better I/O models or hyperparameter
searches. The first column of Figure 1 illustrates the impact of application
modeling errors with an example hyperparameter search over two XGBoost
parameters on the Theta dataset. The search finds that 32 trees with depth of
21 perform best, while the XGBoost defaults use 100 trees of depth 6.
Approximation errors cannot be classified as epistemic or aleatory in nature,
because no new features or jobs are necessary to remove this error. To
estimate AU and EU in the dataset, methods such as AutoDEUQ [21] require first
that an appropriate model architecture is found and trained, so approximation
errors are the first branch of the taxonomy.

System modeling errors: system behavior changes over time due to transient or
long-term changes such as filesystem metadata issues, failing components, new
provisions, etc. [22]. A model that is only aware of application behavior, but
not of system state implicitly assumes that the process is stationary. It will
be forced to learn the average system response to I/O patterns, and will
suffer greater prediction errors during periods when system behavior is
perturbed. Errors that occur due to poor modeling of the global system
component $\zeta_{g}(t)$ are called system modeling errors. To illustrate this
class of errors, two models are trained to predict I/O throughput in the
second column of Figure 1, and each model’s weekly average error is plotted
against time. The blue model is exposed only to application behavior, while
the orange model also knows the job start time. During service degradations,
the blue model has long periods of biased errors while the orange model does
not, since it knows when degradations happen.

Generalization errors:  ML models generally perform well on data drawn from
the same distribution from which their training set was collected. When
exposed to samples highly dissimilar from their training set, the same models
tend to make mispredictions. These samples are called ‘out-of-distribution’
(OoD) because they come from new, shifted, distributions, or the training set
does not have full coverage of the sample space. As an example, the third
column of Figure 1 shows model error before (green) and after (red)
deployment, with the error significantly rising when the model is evaluated on
data collected outside the training time span.

Contention errors: a diverse and variable number of applications compete for
compute, networking, and I/O bandwidth on HPC systems and interact with each
other through these shared resources [17, 23]. Although the global system
state will impact all jobs equally, the impact of resource sharing is specific
to pairs of jobs that are interacting and is harder to observe and model.
Prediction errors that occur due to lack of visibility into job interactions
are called contention errors and are shown in the fourth column of Figure 1.
Here, the I/O throughputs of a number of identical runs (the same code and
data) of different applications illustrate that some applications are more
sensitive to contention than others, even when accounting for global system
state.

Inherent noise errors: while hard to measure, resource sharing errors can
potentially be removed through greater insight into the system and workloads.
What fundamentally cannot be removed are inherent noise errors: errors due to
random behavior by the system (e.g., dropped packets, randomness introduced
through scheduling, etc.). Inherent noise is problematic both because ML
models are bound to make errors on samples affected by noise and because noisy
samples may impede model training. The fifth column of Figure 1 shows the I/O
throughput and start time differences between pairs of identical jobs. The
leftmost column contains jobs that ran at exactly the same time, which often
experience 5% or more difference in I/O throughput.

##  V Datasets and experimental setup

This work is evaluated on two datasets, one collected from ALCF Theta
supercomputer in the period from 2017 to 2020, and one collected from NERSC
Cori supercomputer in the period from 2018 to 2019. Theta collects Darshan
[24] and Cobalt logs and consists of about 100K jobs with an I/O volume larger
than 1GiB, while Cori collects Darshan and Lustre Monitoring Tools (LMT) logs,
and consists of 1.1M jobs larger than 1GiB.

Darshan is an HPC I/O characterization tool that collects HPC job I/O access
patterns on both the POSIX and MPI-IO levels, and serves as our main insight
into application behavior. It collects POSIX aggregate job-level data, e.g.,
total number of bytes transferred, accesses made, read / write ratios, unique
or shared files opened, distribution of accesses per access size, etc. MPI-IO
is a library built on top of POSIX that offers higher-level primitives for
performing I/O operations and may offer the model more insight into
application behavior. Darshan collects MPI-IO information for jobs that use
it, and all requests through MPI-IO are also visible on the POSIX level. The
Cobalt scheduler logs number of nodes and cores assigned to a job, job start
and end times, job placement, etc. Cobalt logs may be useful to the model
since the number of cores a job allotted is not visible to Darshan. Darshan
only collects the number of job processes, which is commonly equal to or
greater than the number of cores allocated to a job. LMT collects I/O
subsystem information such as storage server load and file system utilization,
and serves as our main insight into the I/O subsystem state over time. LMT
records the state of object storage servers (OSS) and targets (OST), and
metadata servers (MDS) and targets (MDT) of the Lustre scratch filesystem
every 5 seconds. Some of the features LMT collects are CPU and memory
utilization of the OSS’s and MDS’s, number of bytes transferred to and from
the OSTs, the fullness of the system, or the number of metadata operations
(e.g., open, close, mkdir, etc.) performed by the metadata targets, etc. LMT
collects per-OSS/OST/MDS/MDT logs, but since a job may be served by an
arbitrary number of these I/O nodes, only the minimum, maximum, mean and
standard deviation are exposed to the ML model. Overall, models have access to
48 Darshan POSIX, 48 Darshan MPI-IO, 37 LMT, and 5 Cobalt features.

The ML models in this work are trained using supervised learning on the task
of predicting the I/O throughput of individual HPC jobs. The error models are
optimizing is:

$e(y,\hat{y})=\frac{1}{n}\sum^{n}_{i=1}\left|log_{10}\left(\frac{y_{i}}{\hat{y}_{i}}\right)\right|$
(6)

where $y_{i}$ and $\hat{y}_{i}$ are the $i$-th job’s measured and predicted
I/O throughputs. Because $log(x)=-log(1/x)$, if a model overestimates or
underestimates the I/O throughput by the same relative amount, the error
remains the same. We use percentages to write errors, where, e.g., a -25%
error specifies that the model underestimated real I/O throughput by 25%;
however, some figures show the absolute error when model bias is not
important. While models try to lower mean error, they report median values
since some of the distributions have heavy tails that make mean estimates
unreliable.

##  VI Application modeling errors

When an ML practitioner is tasked with a regression problem (predicting
continuous values), the first model they evaluate will likely under-perform on
the task, e.g., due to inadequate data preprocessing, architecture, or
hyperparameters. Therefore, the model will suffer from approximation errors,
which can be removed by tuning the model hyperparameters. Because
inappropriate models may be unable to fit even easy tasks, this class of
errors should be solved before seeking, e.g., additional samples or sample
features.

Figure 1 shows that the approximation consists of the application and system
approximation. This section asks whether I/O models build faithful
representations of application behavior. The questions we ask are as follows:
What are the limits of I/O application modeling? In practice, do I/O models
faithfully learn application behavior? Can I/O application modeling benefit
from extra hyperparameter fine-tuning or new application features?

###  VI-A Estimating limits of application modeling

We develop an application modeling error litmus test whose objective is to
separate the application modeling error $e_{app}$ from the four other error
classes in Equation 5. To do so, we seek a ‘golden model’ that models
application behavior as accurately as possible given the inputs. When
practical ML models are compared with this golden model, an estimate of the
application modeling error can be made.

To build this ‘golden model’, we rely on a property of synthetic datasets
where the data-generating process can be freely and repeatedly sampled. When
analyzing HPC logs, it is common to see records of the same application ran
multiple times on the same data, or data of the same format. For example,
system benchmarks such as IOR [25] may be run periodically to evaluate
filesystem health and overall performance. We call these sets of repeated jobs
‘duplicate jobs’. Jobs are duplicates if they belong to the same application
and all their observable application features are identical. Note that an
application can have many different sets of duplicate jobs, likely ran with
different input parameters. Because the duplicate features fed to an ML model
are identical, the model cannot distinguish between duplicates from the same
set, and will achieve best possible accuracy on the training set if it learns
to predict the mean I/O throughput value of a set. A model that does not learn
to predict a set’s mean value suffers from application-modeling error.

Sets of duplicate jobs can be used to build a litmus test that evaluates the
median absolute error of a model which has $e_{app}=0$. Any practical model
with a median error $e^{p}$ can then learn its application modeling error
$e^{p}_{app}$ by comparing it with the median error of the ‘golden model’
$e^{g}$ as $e^{p}_{app}=e^{p}-e^{g}$. The litmus test is administered as:

1\. Sets of duplicate jobs in the dataset are found; 2. The mean I/O
throughput of each set is calculated; 3\. This mean is subtracted from each
duplicate’s I/O throughput in the set to calculate duplicate error, and
Bessel’s correction is applied [26]; 4\. The median error across all
duplicates is reported.

This average difference represents the smallest average error that a model
$m(j)$ can achieve on duplicate jobs. Assuming that duplicate jobs are drawn
from the same distribution of applications as the whole dataset, the duplicate
median absolute error represents the lower bound on median absolute error a
model can achieve on the whole dataset. Note that different applications may
have different distributions of duplicate I/O throughputs, as shown in the
fourth column of Figure 1. For this litmus test to be accurate, a large sample
of applications representative of the HPC system workload must be acquired.
When applied to Theta, 19010 duplicates (23.5% of the dataset) over 3509 sets
show a median absolute error of 10.01%. Cori has 504920 duplicates (54%) in
77390 sets show a median absolute error of 14.15%. If the litmus test is
correct, a model that has the same (but not lower!) median absolute error on
the general dataset can be found.

###  VI-B Minimizing application modeling error

The next question is whether ML models can practically reach the estimate of
the lower bound of the error. Several I/O modeling works have explored
different types of ML models: linear regression [2], decision trees [27],
gradient boosting machines [2, 27, 28], Gaussian processes [4], neural
networks [5], etc. Here, we explore two types of models: XGBoost [29], an
implementation of gradient boosting machines, and feedforward neural networks.
These model types are chosen for their accuracy, scalability, and previous
success in I/O modeling.

Neither type of model achieves ideal performance ‘out of the box’. XGBoost
model performance can be improved through hyperparameter tuning, e.g., by
exploring different (1) numbers of decision trees, (2) their depth, (3) the
features each tree is exposed to, and (4) part of the dataset each tree is
exposed to. Neural networks are more complex, since they require tuning
hyperparameters (learning rate, weight decay, dropout, etc.), while also
exploring different architectures (number, size, type of layers, and their
connectivity). In the case of XGBoost, we exhaustively explore four
hyperparameters listed above, for a total of 8046 XGBoost models. In the case
of neural networks, exhaustive exploration is not feasible due to state space
explosion, so we use AgEBO [30], a Network Architecture Search (NAS) method
that trains populations of neural networks and updates each subsequent
generation’s hyperparameters and architectures through Bayesian logic.

The leftmost column of Figure 1 shows a heatmap of an exhaustive search over
two parameters on the Theta dataset, with the other two parameters (% of
columns and rows revealed to the trees) selected from the best possible result
found. The best performing model has an error rate of 10.51% - very close to
the predicted bound. The search on Cori data arrives at a similar
configuration (omitted due to space).

In the case of neural networks, Figure 2 shows a scatter plot of test set
errors of 10 generations of neural networks on the Cori system, with 30
networks per generation. The networks are evolved using a separate validation
test to prevent leakage of the test set into the model parameters. Networks
approach the estimated error limit, and the best result achieves a median
absolute error of 14.3%. The fact that (1) after extensive tuning both types
of models asymptotically approach the estimated limit in model accuracy, and
(2) NAS does little to improve models, since models improve predictions only 6
times (gold stars). This suggests that both types of ML models are impeded by
the same barrier and that the architecture and the tuning of models are not
the fundamental issue in achieving better accuracy, i.e., that the source of
error lies elsewhere.

Figure 2: Results of the Neural Architecture Search (NAS), with the estimated
error lower bound highlighted in red.

###  VI-C Increasing visibility into applications

While hyperparameter and architecture searches approach but do not surpass the
litmus test’s estimated lower bound on error, this is not conclusive evidence
that all application modeling error has been removed and that error stems from
other sources. Possibly, there exist missing application features that might
further reduce errors. We explore two such sets of features: MPI-IO logs and
Cobalt scheduler logs.

Figure 3 shows the absolute error distribution of the tuned models on three
Theta datasets: POSIX, POSIX + MPI-IO, and POSIX + Cobalt (Cori excluded
because of the lack of Cobalt logs). None of the dataset enrichments help
reduce error, corroborating the conclusion that application modeling is not a
source of error for these models, and further insight into applications will
not help. Adding Cobalt logs does reduce the error on the training set, and
our experiments show that the job start and end time features are the cause.
Once timing features are present in the dataset, no two jobs are duplicates
due to small timing variations. Although previously the ML model was not able
to overfit the dataset due to the existence of duplicates, this is no longer
the case, and the ML model can differentiate and memorize each individual
sample. In [2] authors remove timing features for a similar reason: ML models
can learn Darshan’s implementation of I/O throughput calculation and achieve
good predictions without having to observe job behavior.

Figure 3: Error distributions for models trained on POSIX, POSIX + MPI-IO, and
POSIX + Cobalt feature sets.

##  VII Global system modeling errors

The second part of the approximation error in our taxonomy is the global
system modeling error. This error refers to I/O climate and I/O weather
effects [22] that affect all jobs running on the system, and corresponds to
the second component in Equation 3. While global and local system impacts have
complex and overlapping effects, we suggest that factorizing the impact
applied to all jobs versus the impact that is dependent on pairs of concurrent
jobs is useful for modeling purposes. The main difference between the two is
that local system impacts cannot be predicted or modeled without knowledge of
all jobs running on the system, while global system impacts can. In other
words, global system impacts can be expressed as a property of the system at a
given time, effectively compressing a part of system behavior. We now ask: How
does I/O contention impact job I/O throughput prediction? What are the limits
of global system modeling? How can I/O models approach this limit?

###  VII-A Estimating limits of global system modeling

Global system impact $\zeta_{g}(t)$ on job $j$ from Equation 3 can be
formalized as some function $\zeta_{g}(t)=g(J(t))$ where $J$ is the set of
jobs running at time $t$. Since jobs have a start and end time, given a
dataset with a dense enough sampling of $J$, $g(J(t))$ can be calculated for
every point in time. During periods of time where the file system is suffering
a service degradation, all jobs on the system will be impacted with varying
severity. A model of the system does not need to understand how and why the
degradation happened, it only needs to know the period when it lasted and how
different types of jobs were impacted. This time-based model is useless for
predicting future performance, and it’s only utility is in evaluating how much
of the degradation can be described as purely a function of time. A deployed
model does not have data on the future and will still need to observe the
system.

To evaluate the global system impact, a golden model that exhibits no global
modeling error is developed, against which we can compare other, ‘real’ ML
models. Since the global system impact of $\zeta_{g}(t)$ only depends on time
and does not need the set of all jobs $J$, only application behavior $j$ and
the job start time feature are exposed to the golden model. Here, a golden
model is an XGBoost model fine tuned on a validation set and evaluated on a
test set. Assuming that the golden model is exposed to enough jobs throughout
the lifetime of the system, it will learn the impact of $\zeta_{g}(t)$ even
without having access to the underlying system features causing that impact.
This golden model is used in the following litmus test:

1\. The timing feature is added to the Darshan-only (no Lustre or Cobalt)
dataset; 2. A hyperparameter search is performed on a validation set and
application modeling error is removed; 3. the final model is evaluated on the
test set, and it’s error is reported.

If the litmus test is correct, the golden model only suffers from the last
three classes of errors: poor generalization, local system impact, and
inherent noise. In Figure 4 we evaluate a baseline model (blue) and a model
enriched with the job start time (orange). The addition of this single feature
has a large impact on error: on Cori, the error drops $40\%$, from 16.49% down
to 10.02%, while on Theta the error drops by 30.8%. Note that to obtain high
accuracy on the POSIX+start time, a much larger model is needed, i.e., one
that can remember the I/O weather throughout the lifetime of the system.

###  VII-B Improving modeling through I/O visibility

With an estimate of minimal error achievable assuming perfect application and
global system modeling, we investigate whether I/O subsystem logs can help
models approach this limit. Since Theta does not collect I/O subsystem logs,
we analyze Cori, which collects both application and I/O logs. Figure 4 shows
the XGBoost performance of three models: a baseline where $e_{app}=0$ (blue),
the litmus test golden model where also $e_{system}=0$ (orange), and a Lustre-
enriched model (green). Cori’s median absolute error is reduced by 40%, from
16.49% down to 9.96%. The Lustre-enriched results are surprisingly close to
the litmus test’s predictions, and suggest that predictions cannot be improved
through further I/O insight since the litmus test’s prediction is reached.

Figure 4:  Error distribution of models trained on (1) POSIX, (2) POSIX + the
start time feature, and (3) Darshan and Lustre

##  VIII Generalization errors

The remaining three classes of error are caused by lack of data and not poor
modeling, as the top branch of the taxonomy shows. While I/O contention and
inherent noise errors are examples of aleatory uncertainty and are caused by
lack of insight into specific jobs, generalization errors stem from epistemic
uncertainty, i.e., the lack of other logged jobs around a specific job of
interest. To motivate this section, in the third graph of Figure 1 we show
error distribution of a model trained on data from January 2018 to July 2019.
When evaluated on held-out data from the same period, the median absolute
error is low (green line). Once the model is deployed and evaluated on the
data collected after the training period (July 2019 and after), median error
spikes up (red line).

###  VIII-A Estimating generalization error

Estimating the amount of out-of-distribution error $e_{OoD}$ is important
because any unaccounted OoD error will be classified as noise or contention.
This will make systems that run a lot of novel jobs appear to be more noisy
than they truly are. Because OoD and ID jobs will likely have a similar amount
of I/O and contention noise, it is better to have false positives (ID jobs
classified as OoD) than the other way, since false negatives contribute to
overestimating I/O noise. To estimate the impact of out-of-distribution jobs
on error $e_{OoD}$, we aim to quantify the how much of error is epistemic, and
how much is aleatory in nature, as shown in Figure 1 (upper right). The
leading paradigm for uncertainty quantification works by training an ensemble
of models and evaluating all of the models on the test set. If the models make
the same error, the sample has high aleatory uncertainty, but if the models
disagree, the sample has high epistemic uncertainty [31]. The intuition is
that predictions on out-of-distribution samples will vary significantly on the
basis of the model architecture, whereas predictions on noisy samples will
exhibit the same amount of error. Since this method relies on ensemble model
diversity, several works have explored increasing diversity through different
model hyperparameters [32], different architectures [33], or both [21]. We
choose to use AutoDEUQ [21], a method that evolves an ensemble of neural
network models and jointly optimizes both the architecture and hyperparameters
of the models. AutoDEUQ’s Neural Architecture Search (NAS) is compatible with
the NAS search from section VI, reducing the computational load of applying
the taxonomy. Figure 5 shows the distribution of epistemic (EU) and aleatory
uncertainties (AU) of Theta and Cori test sets. For both systems, aleatoric
uncertainty is significantly higher than epistemic uncertainty. Furthermore,
all jobs seem to have AU larger than some about 0.05, hinting at the inherent
noise present in the system. The inverse cumulative distributions on the top
and right show the what percentage of total error is caused by AU / EU below
that value. For example, for both systems 50% of all error is caused by jobs
with EU below 0.04, while in case of AU, 50% of error is below AU=0.25. The
low total EU is expected since the test set was drawn from the same
distribution as the training set, and increases on the 2020 set (omitted due
to space concerns).

Figure 5: Distribution of job aleatory and epistemic uncertainties for the two
systems, with marginal distributions and inverse cumulative error shown on the
margins.

Epistemic uncertainty does not directly translate into the out-of-distribution
error $e_{OoD}$ from Equation 5. When a sample is truly OoD, it may not be
possible to separate aleatory and epistemic uncertainty, since a good estimate
of AU requires dense sampling around the job of interest. Therefore, we choose
to attribute all errors of a sample marked as out-of-distribution to
$e_{OoD}$. This error attribution requires classifying every test set sample
as either in- or out-of-distribution, but since EU estimates are continuous
values, an EU threshold which will separate OoD and ID samples is required.
Although this threshold is specific to the dataset and may require tuning, the
quick drop or ‘shoulder’ in inverse cumulative error around EU=0.1 in Figure 5
makes the choice of an $e_{OoD}$ threshold robust. A litmus test that
estimates the error due to out-of-order samples has the following steps:

1\. Run NAS and collect the best performing models; 2. Estimate the aleatory
and epistemic uncertainty using AutoDEUQ; 3. Find a stable EU threshold and
classify the samples as ID and OoD; 4. Calculate $e_{OoD}$ as the sum of OoD
sample errors.

On Theta, for an EU threshold of 0.24, .7% of the samples are classified as
OoD, but constitute 2.4% of the errors, while on Cori 2.1% of error gets
removed for the same EU threshold. In other words, the selected jobs have
$3\times$ larger average error than random samples. By manually exploring the
types of jobs that do get removed, we confirm that these are typically rare or
novel applications.

##  IX I/O Contention and Inherent Noise Errors

With the ability to estimate the amount of application and system modeling
error, as well as detect outlier jobs, leftover error is caused by system
contention or inherent noise. Both of these error classes are caused by
aleatory uncertainty, since the model lacks deeper insight into jobs or the
system, as opposed the OoD case where the model lacks samples. While e.g.,
application error was predictable and explainable in terms of broad
application behavior (e.g., this application is slow because it frequently
writes to shared files), the impact of contention and noise on I/O throughput
is caused by lower level, transient effects. Though it may be possible to
observe and log such effects through microarchitectural hardware counters or
network switch logs, such logging would require vast amounts of storage per
job and would impact performance. Lack of practical logging ability makes the
last two error categories typically unobservable. Furthermore, these two
classes may only be separated in hindsight, and while I/O noise levels may be
constant, the amount of I/O contention on the system is unpredictable for a
job that is about to run.

The question we ask in this section are: how can errors due to noise and
contention be separated from errors due to poor modeling or epistemic
uncertainty? Is there a fundamental limit to how accurate I/O models can
become? Can system I/O variability be quantified?

###  IX-A Establishing the bounds of I/O modeling

To separate the contention and noise from the first three classes of error, we
develop a litmus test based on the test from Section VI. There, by observing
sets of duplicates, the performance of an ideal model was estimated for which
$e_{app}=0$. Comparing real models against this ideal model allowed us to
calculate real model’s $e_{app}$. A similar litmus test can be designed
estimate the sum of contention and noise error, where only concurrent
duplicates are observed and both application behavior $j$ and global system
behavior $\zeta_{g}$ can be held static for each duplicate set:

1\. OoD jobs are removed and sets of duplicate jobs ran at the same time
($\Delta t=0$) are collected; 2\. The mean I/O throughput of each set is
calculated; 3\. Duplicate error is calculated as before; 4\. The median error
across all duplicates is reported.

In the fifth column of Figure 1 we show the distribution of I/O throughput
differences $\Delta\phi$ and timing differences $\Delta t$ between all pairs
of Cori duplicate jobs, weighted so that large duplicate sets are not
overrepresented. The vertical strip on the left contains Cori duplicate jobs
that were ran simultaneously, largely because they were batched together.
These jobs share $j$ and $\zeta_{g}$, but may differ in $\zeta_{l}$ and
$\omega$. Due to the denser sampling around 1 minute to 1 hour range, it is
not immediately apparent how the I/O difference changes between duplicates ran
at the same time and duplicates ran with a delay. By grouping duplicates from
different $\Delta t$ ranges and independently scaling them, a better
understanding of duplicate I/O throughput distributions across timescales can
be made, as shown in Figure 6. For both systems (Cori omitted due to lack of
space), the distributions on the right contain jobs ran over large periods of
time where global system impact $\zeta_{g}$ might have changed, explaining the
asymmetric shape of some of them. The left-most distributions are similar,
since variance only stems from contention $\zeta_{l}$ and noise $\omega$.
While some distributions (e.g., the $10^{5}$ to $10^{6}$ second) show complex
multimodal behavior, all of the distributions seem to contain the initial zero
second ($0s$ to $1s$) distribution.

Figure 6: Distribution of errors for different periods between duplicate runs.

By fitting a normal distribution to the $\Delta t=0$ distribution (0s to 1s)
in Figure 6, we can both (1) learn the lower limit on total modeling error and
(2) learn the system’s I/O noise level, i.e., how much I/O throughput variance
should jobs running on the system expect. However, upon closer inspection, the
$\Delta t=0$ distribution does not follow a normal distribution. This is
surprising, since if noise is normally distributed, independent over time, and
it’s effects are cumulative, the total impact is a sum of normal
distributions, which should also be a normal distribution. The answer lies in
how the concurrent ($\Delta t=0$) duplicates are sampled. When observing
duplicates, in general, duplicate sets have between 2 and hundreds of
thousands of identical jobs in them. However, in duplicate sets with identical
start times on Theta, 70% of the sets only have two identical jobs, and 96%
have 6 jobs or less, with similar results on Cori. The issue stems from how
small (sub-30 samples) duplicate set errors are calculated: when only a small
number of jobs exist in the set, the mean I/O throughput of the set is biased
by the sampling, i.e., the estimated mean is closer to the samples than the
real mean is. This causes the set I/O throughput variance to decrease and
therefore duplicate error estimate will be reduced as well. Student’s
t-distribution describes this effect: when the true mean of a distribution is
known, error calculations follow a normal distribution. When the true mean is
not known, the biased mean estimate makes the error follow the t-distribution.
As the sample count reaches 30 and above, the t-distribution approaches the
normal distribution.

Finally, we seek to estimate the I/O noise variance of a given system. Naively
taking the variance of the t-distribution will produce a biased sample
variance $\sigma^{2}$, but it can be fixed by applying Bessel’s correction as
$\frac{n}{n-1}\sigma^{2}$. In practical terms, a job running on Theta can
expect an I/O throughput within $\pm 5.71\%$ of the predicted value 68% of the
time, or within $\pm 10.56\%$ 95% of the time. For Cori, these values are $\pm
7.21\%$ and $\pm 14.99\%$, respectively. This is a fundamental barrier not
just to I/O model improvement, but to predictable system usage in general.
Although some insight into contention can be gained through low-level logging
tools, noise cannot be overcome. I/O practitioners can use this litmus test to
evaluate the noise levels of their systems, and ML practitioners should
reconsider how they evaluate models, since some systems may be simply harder
to model.

##  X Applying the taxonomy

Figure 7: Framework for applying the taxonomy (left), and the results from
ALCF Theta and NERSC Cori systems.

We now illustrate how the proposed taxonomy can be used in practice. In Figure
7 (left), we show the steps a modeler can follow to evaluate the taxonomy on a
new system. Step 1: The modeler splits the available data into training and
test sets, and then trains and evaluates some baseline machine learning model
on the task of predicting I/O throughput. This model does not have to be
perfect, as the taxonomy should reveal what are the main sources of error and
approximately how much the quality of the model is at fault. Step 2.1: The
modeler estimates application modeling errors by finding duplicate jobs and
evaluating the mean predictor performance on every set of duplicates. Assuming
that the distribution of duplicate HPC jobs is representative of the whole
population of jobs, this step provides the modeler with a lower bound on the
application modeling error. Step 2.2: By contrasting the baseline model error
(Step 1) and the estimated application modeling error, the modeler can
estimate the percentage of error that can be attributed to poor modeling. The
modeler performs a hyperparameter or network architecture search and arrives
at a good model close to the bound. Step 3.1: The modeler estimates system
modeling errors by exposing the job start time to a golden model. This step
requires that the modeler has developed a well-performing model in Step 2.2,
that is, one that achieves close to the estimated ideal performance. The test
set error of the model serves as an estimate of the application \+ system
modeling lower bound. Step 3.2: The modeler explores adding sources of system
data to improve the performance of the baseline model up to the estimated
limit of application and system modeling. Step 4: The modeler identifies out-
of-distribution samples using AutoDEUQ, calculates OoD error that stems from
these samples, and removes them from the dataset. Step 5: The modeler
estimates the error that can be attributed to contention and noise, as well as
I/O variance of the system. This estimate is made by observing the I/O
throughput differences between sets of concurrent duplicates, i.e., duplicate
jobs ran at around the same time.

In the middle and right portion of Figure 7 we show the average baseline model
error (inner pink circle segment) of both ANL Theta and NERSC Cori systems,
and how that error is broken down into different classes of error. We do not
focus on the cumulative error value of the two systems; instead, we focus on
attributing the baseline model error into the five classes of errors in the
taxonomy (middle circle segments of the pie chart), and how much improved
application and system modeling can help reduce the cumulative error (outer
segments of the pie chart). The inner blue section of the two pie charts
represents the estimated application modeling error, as arrived at in Step
2.1. The outer blue section represents how much of the error can be fixed
through hyperparameter exploration, as explored in Step 2.2. The inner green
section represents the estimated system modeling error, derived in Step 3.1.
Note that the total percentage of system modeling error is relatively small on
both systems; i.e., I/O contention, filesystem health, hardware faults, etc.,
do not have a dominant impact on I/O throughput. The outer green circle
segment represents the percentage of error that can be fixed by including
system logs (LMT logs in our case), as described in Step 3.2. Only the Cori
pie chart has this segment, as Theta does not collect LMT logs. On Cori, the
inclusion of LMT logs helps remove most of the system modeling errors,
pointing to the conclusion that including other logs (i.e., topology,
networking) may not help to significantly reduce errors. The inner red segment
represents the percentage of error that can be attributed to out-of-
distribution samples of the two systems, as calculated in Step 4. Finally, the
yellow circle segment represents the percentage of error that can be
attributed to aleatory uncertainty. For both Theta and Cori, this is a rather
large amount, pointing to the fact that there exists a lot of innate noise in
the behavior of these systems, and setting a relatively high lower bound on
ideal model performance.

The similarity between the modeling error estimates (Steps 2.1 and 3.1) and
the actual updated model performance (Steps 2.2 and 3.2) is surprising and
serves as evidence for the quality of the error estimates. However, the
estimates of the five error classes do not add up to 100%. The first three
error estimates are just that - estimates, derived from a subset of data
(duplicate HPC jobs) that do not necessarily follow the same distribution as
the rest of the dataset and may be biased. If we add the estimates, we see
that on Theta 32.9% of the error is unexplained, and on Cori 13.5% of the
error is unexplained. Cori’s lower unexplained error may be due to the fact
that we have collected some 1.1M jobs compared to 100K on Theta.

##  XI Discussion

Developing production-ready machine learning models that analyze HPC jobs and
predict I/O throughput is difficult: the space of all application behaviors is
large, HPC jobs are competing for resources, and the system changes over time.
To efficiently improve these models, we present a taxonomy of HPC I/O modeling
errors that allows us to study the types of errors independently, quantify
their impact, and identify the most promising avenues for model improvement.
Our taxonomy breaks errors into five categories: (1) application and (2)
system modeling errors, (3) poor generalization, (4) resource contention, and
(5) I/O noise. We present litmus tests that quantify what percentage of
existing error belongs to the classes and show that models improved by using
the taxonomy are within several percentage points of an estimated best-case
I/O throughput modeling accuracy. We show that a large portion of I/O
throughput modeling error is irreducible and stems from I/O variability. We
provide tests that quantify the I/O variability and establish an upper bound
on how accurate models can become. Our test shows that jobs run on Theta and
Cori can expect an I/O throughput standard deviation of 5.7% and 7.2%,
respectively.

## References

  * [1] H. Luu, M. Winslett, W. Gropp, R. Ross, P. Carns, K. Harms, M. Prabhat, S. Byna, and Y. Yao, “A Multiplatform Study of I/O Behavior on Petascale Supercomputers,” in _Proceedings of the 24th International Symposium on High-Performance Parallel and Distributed Computing_. ACM, Jun. 2015, pp. 33–44.
  * [2] M. Isakov, E. d. Rosario, S. Madireddy, P. Balaprakash, P. Carns, R. B. Ross, and M. A. Kinsy, “HPC I/O Throughput Bottleneck Analysis with Explainable Local Models,” ser. SC ’20. IEEE Press, 2020.
  * [3] K. W. Cameron, A. Anwar, Y. Cheng, L. Xu, B. Li, U. Ananth, J. Bernard, C. Jearls, T. Lux, Y. Hong, L. T. Watson, and A. R. Butt, “MOANA: Modeling and Analyzing I/O Variability in Parallel System Experimental Design,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 30, no. 8, pp. 1843–1856, 2019.
  * [4] S. Madireddy, P. Balaprakash, P. Carns, R. Latham, R. Ross, S. Snyder, and S. M. Wild, “Machine Learning Based Parallel I/O Predictive Modeling: A Case Study on Lustre File Systems,” in _High Performance Computing_ , 2018, pp. 184–204.
  * [5] S. Madireddy, P. Balaprakash, P. Carns, R. Latham, G. K. Lockwood, R. Ross, S. Snyder, and S. M. Wild, “Adaptive Learning for Concept Drift in Application Performance Modeling,” in _Proceedings of the 48th International Conference on Parallel Processing_ , ser. ICPP 2019, 2019.
  * [6] M. Isakov, E. d. Rosario, S. Madireddy, P. Balaprakash, P. Carns, R. B. Ross, and M. A. Kinsy, “Toward Generalizable Models of I/O Throughput,” in _2020 IEEE/ACM International Workshop on Runtime and Operating Systems for Supercomputers (ROSS)_ , 2020, pp. 41–49.
  * [7] S. Kim, A. Sim, K. Wu, S. Byna, Y. Son, and H. Eom, “Towards HPC I/O Performance Prediction through Large-Scale Log Analysis,” in _Proceedings of the 29th International Symposium on High-Performance Parallel and Distributed Computing_ , ser. HPDC ’20, 2020.
  * [8] E. del Rosario, M. Currier, M. Isakov, S. Madireddy, P. Balaprakash, P. Carns, R. B. Ross, K. Harms, S. Snyder, and M. A. Kinsy, “Gauge: An Interactive Data-Driven Visualization Tool for HPC Application I/O Performance Analysis,” in _2020 IEEE/ACM Fifth International Parallel Data Systems Workshop (PDSW)_ , 2020, pp. 15–21.
  * [9] E. Ates, O. Tuncer, A. Turk, V. J. Leung, J. Brandt, M. Egele, and A. K. Coskun, “Taxonomist: Application Detection Through Rich Monitoring Data,” in _Euro-Par 2018: Parallel Processing_ , 2018.
  * [10] S. Kim, A. Sim, K. Wu, S. Byna, T. Wang, Y. Son, and H. Eom, “DCA-IO: A Dynamic I/O Control Scheme for Parallel and Distributed File Systems,” in _2019 19th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing (CCGRID)_ , 2019, pp. 351–360.
  * [11] Y. Li, K. Chang, O. Bel, E. L. Miller, and D. D. E. Long, “CAPES: Unsupervised Storage Performance Tuning Using Neural Network-Based Deep Reinforcement Learning,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , ser. SC ’17. New York, NY, USA: Association for Computing Machinery, 2017. [Online]. Available: https://doi.org/10.1145/3126908.3126951
  * [12] F. Isaila, P. Balaprakash, S. M. Wild, D. Kimpe, R. Latham, R. Ross, and P. Hovland, “Collective I/O Tuning Using Analytical and Machine Learning Models,” in _Proceedings of the 2015 IEEE International Conference on Cluster Computing_ , ser. CLUSTER ’15, 2015.
  * [13] T. Wang, S. Snyder, G. K. Lockwood, P. Carns, N. J. Wright, and S. Byna, “IOMiner: Large-Scale Analytics Framework for Gaining Knowledge from I/O Logs,” _2018 IEEE International Conference on Cluster Computing (CLUSTER)_ , pp. 466–476, 2018.
  * [14] D. Skinner and W. Kramer, “Understanding the causes of performance variability in HPC workloads,” in _Proceedings of the IEEE Workload Characterization Symposium_ , 2005, pp. 137–149.
  * [15] G. K. Lockwood, S. Snyder, T. Wang, S. Byna, P. Carns, and N. J. Wright, “A Year in the Life of a Parallel File System,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis_ , ser. SC ’18. IEEE Press, 2018.
  * [16] E. Costa, T. Patel, B. Schwaller, J. M. Brandt, and D. Tiwari, “Systematically Inferring I/O Performance Variability by Examining Repetitive Job Behavior,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , ser. SC ’21. New York, NY, USA: Association for Computing Machinery, 2021. [Online]. Available: https://doi.org/10.1145/3458817.3476186
  * [17] B. Li, S. Chunduri, K. Harms, Y. Fan, and Z. Lan, “The Effect of System Utilization on Application Performance Variability,” in _Proceedings of the 9th International Workshop on Runtime and Operating Systems for Supercomputers_ , ser. ROSS ’19, 2019.
  * [18] L. Wan, M. Wolf, F. Wang, J. Y. Choi, G. Ostrouchov, and S. Klasky, “Comprehensive Measurement and Analysis of the User-Perceived I/O Performance in a Production Leadership-Class Storage System,” in _2017 IEEE 37th International Conference on Distributed Computing Systems (ICDCS)_ , 2017, pp. 1022–1031.
  * [19] B. Xie, Y. Huang, J. S. Chase, J. Y. Choi, S. Klasky, J. Lofstead, and S. Oral, “Predicting Output Performance of a Petascale Supercomputer,” in _Proceedings of the 26th International Symposium on High-Performance Parallel and Distributed Computing_ , ser. HPDC ’17. New York, NY, USA: Association for Computing Machinery, 2017, p. 181–192. [Online]. Available: https://doi.org/10.1145/3078597.3078614
  * [20] S. Madireddy, P. Balaprakash, P. Carns, R. Latham, R. Ross, S. Snyder, and S. Wild, “Modeling I/O Performance Variability Using Conditional Variational Autoencoders,” in _2018 IEEE International Conference on Cluster Computing (CLUSTER)_ , 2018, pp. 109–113.
  * [21] R. Egele, R. Maulik, K. Raghavan, P. Balaprakash, and B. Lusch, “AutoDEUQ: Automated Deep Ensemble with Uncertainty Quantification,” _CoRR_ , vol. abs/2110.13511, 2021. [Online]. Available: https://arxiv.org/abs/2110.13511
  * [22] G. K. Lockwood, W. Yoo, S. Byna, N. J. Wright, S. Snyder, K. Harms, Z. Nault, and P. Carns, “UMAMI: A Recipe for Generating Meaningful Metrics through Holistic I/O Performance Analysis,” in _Proceedings of the 2nd Joint International Workshop on Parallel Data Storage & Data Intensive Scalable Computing Systems_, ser. PDSW-DISCS ’17, 2017.
  * [23] X. Yang, J. Jenkins, M. Mubarak, R. B. Ross, and Z. Lan, “Watch Out for the Bully! Job Interference Study on Dragonfly Network,” in _SC ’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2016, pp. 750–760.
  * [24] S. Snyder, P. Carns, K. Harms, R. Ross, G. K. Lockwood, and N. J. Wright, “Modular HPC I/O Characterization with Darshan,” in _5th Workshop on Extreme-Scale Programming Tools (ESPT)_ , 2016.
  * [25] H. Shan and J. Shalf, “Using IOR to analyze the I/O Performance for HPC Platforms,” 6 2007. [Online]. Available: https://www.osti.gov/biblio/923356
  * [26] C. M. Bishop, _Pattern Recognition and Machine Learning (Information Science and Statistics)_. Berlin, Heidelberg: Springer-Verlag, 2006.
  * [27] O. Tuncer, E. Ates, Y. Zhang, A. Turk, J. Brandt, V. J. Leung, M. Egele, and A. K. Coskun, “Diagnosing Performance Variations in HPC Applications Using Machine Learning,” in _High Performance Computing_ , 2017, pp. 355–373.
  * [28] B. Xie, Z. Tan, P. Carns, J. Chase, K. Harms, J. Lofstead, S. Oral, S. S. Vazhkudai, and F. Wang, “Interpreting Write Performance of Supercomputer I/O Systems with Regression Models,” in _2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_ , 2021, pp. 557–566.
  * [29] T. Chen and C. Guestrin, “XGBoost: A Scalable Tree Boosting System,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , ser. KDD ’16, 2016.
  * [30] R. Egele, P. Balaprakash, V. Vishwanath, I. Guyon, and Z. Liu, “AgEBO-Tabular: Joint Neural Architecture and Hyperparameter Search with Autotuned Data-Parallel Training for Tabular Data,” _CoRR_ , vol. abs/2010.16358, 2020\. [Online]. Available: https://arxiv.org/abs/2010.16358
  * [31] B. Lakshminarayanan, A. Pritzel, and C. Blundell, “Simple and Scalable Predictive Uncertainty Estimation Using Deep Ensembles,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , ser. NIPS’17. Red Hook, NY, USA: Curran Associates Inc., 2017, p. 6405–6416.
  * [32] F. Wenzel, J. Snoek, D. Tran, and R. Jenatton, “Hyperparameter Ensembles for Robustness and Uncertainty Quantification,” _ArXiv_ , vol. abs/2006.13570, 2020.
  * [33] S. Zaidi, A. Zela, T. Elsken, C. C. Holmes, F. Hutter, and Y. W. Teh, “Neural Ensemble Search for Uncertainty Estimation and Dataset Shift,” 2020.


11institutetext: Télécom Paris, Institut Polytechniques de Paris, France
11email<EMAIL_ADDRESS>
# A Self-Encoder for Learning Nearest Neighbors

Armand Boschin  Thomas Bonald  Marc Jeanmougin

###### Abstract

We present the self-encoder, a neural network trained to guess the identity of
each data sample. Despite its simplicity, it learns a very useful
representation of data, in a self-supervised way. Specifically, the self-
encoder learns to distribute the data samples in the embedding space so that
they are linearly separable from one another. This induces a geometry where
two samples are close in the embedding space when they are not easy to
differentiate. The self-encoder can then be combined with a nearest-neighbor
classifier or regressor for any subsequent supervised task. Unlike regular
nearest neighbors, the predictions resulting from this encoding of data are
invariant to any scaling of features, making any preprocessing like min-max
scaling not necessary. The experiments show the efficiency of the approach,
especially on heterogeneous data mixing numerical features and categorical
features.

###### Keywords:

Neural network, nearest neighbors, linear invariance.

##  1 Introduction

Despite the recent progress of machine learning, the question of the optimal
encoding of data remains open, especially for tabular data [2]. In this paper,
we present the self-encoder, a neural network trained to guess the identity of
each data sample. Given $n$ data samples $x_{1},\ldots,x_{n}\in{\bbbr}^{d}$,
the self-encoder maps any sample $x\in{\bbbr}^{d}$ to a probability
distribution $p$ over $\\{1,\ldots,n\\}$, in such a way that $p(x_{i})$ is
close to a Dirac in $i$ for each $i\in\\{1,\ldots,n\\}$. In other words, the
self-encoder is a classifier where each sample of the train set has its own
label (its index). As such, it belongs to the category of self-supervised
learning methods, like auto-encoders. The key difference is that, while auto-
encoders rely on a reconstruction task, with the output in the same space
${\bbbr}^{d}$ as the original sample, our self-encoder relies on an
identification task, with the output in the set of probability distributions
on the set of indices $\\{1,\ldots,n\\}$.

Despite its simplicity, the self-encoder learns a very useful representation
of data. It learns to distribute data in the embedding space in a way that
makes them linearly separable from one another. This induces a geometry where
two samples are close in the embedding space when they are not easy to
differentiate. In particular, the self-encoder can be used for any
classification or regression task, using the $k$ nearest neighbors in the
sense of this geometry, as given by the ranking of the predicted probabilities
$p_{1}(x),\ldots,p_{n}(x)$, for any sample $x$. Interestingly, these nearest
neighbors do not correspond to those given by the Euclidean distance in the
original space ${\bbbr}^{d}$ (nor by any other usual metric like a Minkowski
metric or cosine similarity for instance). The geometry is learned by the
model. In particular, the predictions resulting from this encoding of data are
invariant to any scaling of the features, making any preprocessing like min-
max scaling not necessary.

A drawback of the self-encoder is its complexity, as the dimension of the
output is equal to $n$, the size of the training set. This induces a time
complexity in $O(n^{2})$ for training. To overcome this problem, we present a
variant based on sampling where the model is trained to predict the identity
of samples in a random subset of the training set, reducing the training time.

The rest of the paper is organized as follows. We present the related work in
Section 2. The self-encoder is presented in Section 3. In Section 4, we prove
that the learned geometry is invariant to linear transformations, like any
scaling of the features. The behavior of the self-encoder in the presence of
categorical features is analyzed in Section 5. The variant of the model based
on sampling is described in Section 6. The experiments are presented Section
7. Section 8 concludes the paper.

##  2 Related work

###  2.0.1 Nearest Neighbors.

A simple and yet fundamental method to solve tasks in machine learning is the
proximity search. It relies on the intuition that close vectors in the feature
space should have close properties (similar labels in classification and close
values in regression). The most common classifier is the $k$ nearest neighbor
($k$-NN) algorithm [4], which assigns a label to a point by choosing the most
present label among its $k$ nearest neighbors. Usual similarity measures are
the Euclidean distance or the cosine similarity but $k$-NN can also be applied
to any similarity measure, such as that learned by the proposed self-encoder
model.

###  2.0.2 Kernel methods.

Many machine learning algorithms are able to work with kernels, i.e.,
similarity functions that differ from the usual Euclidean metric. This is the
case of the Support Vector Machine (SVM) [3], a classifier that tries to find
linear hyperplanes separating samples of different labels. When the dataset is
not linearly separable, the SVM can still be used with the help of a kernel
that moves the training vectors into another space, possibly of different
dimension, in which they are more likely to be separable. The key point is
that the kernel must be chosen by the user. In comparison, the self-encoder
learns a geometry that makes data samples linearly separable from another.
Moreover, it is a self-supervised learning technique, that does not rely on
labeled data.

###  2.0.3 Auto-encoders.

In its simplest form, an auto-encoder is a neural network model that learns a
compact latent representation of data points with an encoder part and that is
able to reconstruct them from this representation with a decoder part. The
rationale is that a performing encoder should be able to extract the core
discriminative features that define the samples of the dataset and this is
measured by the ability of the decoder to reconstruct the original data from
those features.

There is a large variety of application domains for auto-encoders and even
more different structures. Let us cite for example variational auto-encoders
that are used in a probabilistic framework for variational inference [6] or
denoising auto-encoders that can be used for example to improve image
resolutions by changing the reconstruction criterion [8] of manually noised
data samples.

The output of an auto-encoder being a reconstructed version of the input
vector, the measure of reconstruction loss is key in the training process.
Unlike auto-encoders, in the case of the Self-Encoder, there is no need to
engineer a good similarity measure between the original data sample and its
reconstruction because the objective is simply to predict the identity of each
data sample as an index in $\\{1,\ldots,n\\}$ and not the data point in
${\bbbr}^{d}$ itself.

##  3 The self-encoder

Let $x_{1},\ldots,x_{n}\in{\bbbr}^{d}$ be the set of training samples with $d$
the dimension of the feature space and $n$ the number of samples. The Self-
Encoder is a multi-layer perceptron with input dimension $d$ and output
dimension $n$, trained to predict the identity $i$ of each data sample
$x_{i}$.

###  3.0.1 Hidden layers

The encoder consists of $L$ hidden layers. Each layer $l=1,\ldots,L$ consists
in an affine transformation followed by an activation function:

$\displaystyle h^{(1)}=\phi^{(1)}\left(W^{(1)}x+b^{(1)}\right)$ $\displaystyle
h^{(2)}=\phi^{(2)}\left(W^{(2)}h^{(1)}+b^{(2)}\right)$ $\displaystyle\vdots$
$\displaystyle h^{(L)}=\phi^{(L)}\left(W^{(L)}h^{(L-1)}+b^{(L)}\right)$

where $\phi^{(1)},\ldots,\phi^{(L)}$ are the activation functions, typically
non-linear. The dimensions of the successive outputs, say
$d_{1},\ldots,d_{L}$, are hyper-parameters. The weight matrices
$W^{(1)},\ldots,W^{(L)}$ and the bias vectors $b^{(1)},\ldots,b^{(L)}$ must be
learned.

###  3.0.2 Output layer

The output layer is a fully connected layer with input dimension $d_{L}$ (the
output dimension of the last hidden layer) and output dimension $n$. This
affine transformation is followed by an activation function $\phi$ which is
either a coordinate-wise sigmoid function:

$\forall
u\in{\bbbr}^{n},\quad\phi(u)=\left(\frac{e^{u_{1}}}{1+e^{u_{1}}},\ldots,\frac{e^{u_{n}}}{1+e^{u_{n}}}\right)$

or a SoftMax function:

$\forall
u\in{\bbbr}^{n},\quad\phi(u)=\left(\frac{e^{u_{1}}}{\sum_{i=1}^{n}e^{u_{i}}},\ldots,\frac{e^{u_{n}}}{\sum_{i=1}^{n}e^{u_{i}}}\right)$

The output of the network is then a vector $p=\phi(Wh^{(L)}+b)\in[0,1]^{n}$
that can be interpreted as probabilities: the $i$th component $p_{i}$ is the
probability that the input $x$ corresponds to the training sample $x_{i}$.
Observe that the probabilities sum to 1 only with the SoftMax function (with
the sigmoid function, the probabilities are learned independently for each
training sample $i$ by the output layer). The weight matrix $W$ and the bias
vector $b$ must be learned, together with the other parameters.

###  3.0.3 Loss function

In the following, $f$ denotes the learned function of the network, mapping
sample vectors $x\in{\bbbr}^{d}$ to probability vectors $p\in[0,1]^{n}$. $f$
is a parametric function and its parameters are the weight matrices
$W^{(1)},\ldots,W^{(L)},W$ and the bias vectors $b^{(1)},\ldots,b^{(L)},b$.
Those are learned by minimizing the following Binary Cross Entropy (BCE) loss
by gradient descent:

${\cal L}=-\sum_{i=1}^{n}\left(\log f_{i}(x_{i})+\sum_{j\neq
i}\log(1-f_{j}(x_{i}))\right)$ (1)

###  3.0.4 Interpretation

The Self-Encoder learns a latent representation of data, given by the last
hidden layer, where the $n$ training samples are linearly separable. It is the
role of the output layer to find the hyperplanes (given by the weight matrix
$W$ and the bias vector $b$) separating each training sample.

###  3.0.5 No hidden layer

Observe that the Self-Encoder can also be trained without any hidden layer. It
then reduces to the output layer, i.e. a perceptron with input dimension $d$
and output dimension $n$. Equivalently, the Self-Encoder then consists of $n$
binary logistic regressions (for the sigmoid activation function) or a single
multinomial logistic regression (for the SoftMax activation function).

###  3.0.6 Geometry

In both cases (with or without hidden layers), the Self-Encoder learns a
specific similarity measure in the sense that it can predict the training
samples that are the most similar to any new data sample $x$. This measure
depends on the distribution of the training samples $x_{1},\ldots,x_{n}$ in
the original space ${\bbbr}^{d}$. Any sample $x\in{\bbbr}^{d}$ is said to be
close to the training sample $x_{i}$ if the corresponding predicted
probability $p_{i}=f_{i}(x)$ is close to 1.

##  4 Invariance property

The Self-Encoder is invariant to invertible affine transformations of the
training data, as stated below.

###### Proposition 1

Let $f$ be the mapping learned by the encoder with training samples
$x_{1},\ldots,x_{n}$. For any invertible matrix $M\in{\bbbr}^{d\times d}$ and
vector $v\in{\bbbr}^{d}$, let $\tilde{x}_{1},\ldots,\tilde{x}_{n}$ be the new
training samples obtained by affine transformation $x\mapsto Mx+v$. The new
mapping $\tilde{f}$ defined by:

$\forall\tilde{x}\in{\bbbr}^{d},\quad\tilde{f}(\tilde{x})=f(M^{-1}(\tilde{x}-v))$

minimizes the cross-entropy loss (1) for the training samples
$\tilde{x}_{1},\ldots,\tilde{x}_{n}$. Both encoders are related through the
affine transformation:

$\forall x\in{\bbbr}^{d},\quad f(x)=\tilde{f}(Mx+v).$

In other words, if the training data are the same up to some invertible affine
transformations, so are the mappings learned by the encoder.

###### Proof

The mapping $\tilde{f}$ is the same encoder as $f$ except for the first hidden
layer, whose output $\tilde{h}^{(1)}$ for the input $\tilde{x}$ is given by:

$\displaystyle\tilde{h}^{(1)}\left(\tilde{x}\right)$
$\displaystyle={h}^{(1)}\left(M^{-1}(\tilde{x}-v)\right)$
$\displaystyle=\phi^{(1)}\left(W^{(1)}(M^{-1}(\tilde{x}-v))+b^{(1)}\right)$
$\displaystyle=\phi^{(1)}\left(W^{(1)}M^{-1}\tilde{x}+b^{(1)}-W^{(1)}M^{-1}v\right)$
$\displaystyle=\phi^{(1)}(\tilde{W}^{(1)}\tilde{x}+\tilde{b}^{(1)}),$

for the weight matrix and bias vector:

$\displaystyle\tilde{W}^{(1)}$ $\displaystyle=W^{(1)}M^{-1}$
$\displaystyle\tilde{b}^{(1)}$ $\displaystyle=b^{(1)}-W^{(1)}M^{-1}v.$

The corresponding binary cross-entropy loss is minimized, given that:

$\displaystyle{\cal L}$
$\displaystyle=-\sum_{i=1}^{n}\left(\log\tilde{f}_{i}(\tilde{x}_{i})+\sum_{j\neq
i}\log(1-\tilde{f}_{j}(\tilde{x}_{i}))\right)$
$\displaystyle=-\sum_{i=1}^{n}\left(\log f_{i}(x_{i})+\sum_{j\neq
i}\log(1-f_{j}(x_{i}))\right)$

This invariance property is a clear difference with the usual Euclidean
distance. To illustrate this, Figure 1 shows the Voronoi diagram (regions
formed by the nearest neighbors) associated with $n=4$ points of ${\bbbr}^{2}$
for the Self-Encoder and for the Euclidean distance. On the left, the 4 points
form a square and the Voronoi diagrams coincide, by symmetry. On the right, a
linear transformation is applied; the Voronoi diagram is obtained by the same
linear transformation for the Self-Encoder, while it changes completely for
the Euclidean distance.

This invariance property is handy as it simplifies the pre-processing steps,
which can become tedious and that usually have an impact on the performances
of many classifiers. This will be showed later in the experiments for the
Euclidean $k$-NN classifier.

(a) Euclidean distance

(b) Self-Encoder

Figure 1: Impact of a linear transformation on the Voronoi diagram (regions of
nearest neighbors) formed by $n=4$ points in ${\bbbr}^{2}$.

##  5 Categorical features

We claim that the Self-Encoder robustly handles categorical features, in the
sense that it does not depend on the number of bits they are encoded on,
unlike Euclidean Nearest Neighbor (NN).

Given $x^{(1)},\dots,x^{(n)}$ training samples with categorical features in
$\\{0,1\\}^{d}$, data redundancy in the features does not modify the optimum
of the function. If there is a pair $(k_{1},k_{2})$ such that $\forall
i,x^{(i)}_{k_{1}}=x^{(i)}_{k_{2}}$ then the contributions to the loss
involving these features is duplicated and they should result in the same
corresponding weights in the input layer, which should not impact the learned
geometry. If there is a pair $(k_{1},k_{2})$ such that $\forall
i,x^{(i)}_{k_{1}}=1-x^{(i)}_{k_{2}}$ then the $k_{1}$-th feature is a
invertible affine transformation of the $k_{2}$-th feature and the invariance
property of Proposition 1 tells us that the first case of redundancy applies.

Let us illustrate on a simple example how the similarity measure learned by a
Self-Encoder can outperform the usual Euclidean metric on a nearest-neighbor
search by handling categorical features differently.

$\displaystyle X_{1}$ $\displaystyle=\begin{pmatrix}x_{1}^{(1)}\\\
x_{1}^{(2)}\\\ x_{1}^{(3)}\\\ x_{1}^{(4)}\\\
x_{1}^{(5)}\end{pmatrix}=\begin{pmatrix}0&0&0&1\\\ 0&0&1&0\\\ 0&1&0&0\\\
1&0&0&1\\\ 1&0&1&0\\\ \end{pmatrix}$ (2) $\displaystyle X_{2}$
$\displaystyle=\begin{pmatrix}x_{2}^{(1)}\\\ x_{2}^{(2)}\\\ x_{2}^{(3)}\\\
x_{2}^{(4)}\\\ x_{2}^{(5)}\end{pmatrix}=\begin{pmatrix}1&0&0&0&1\\\
1&0&0&1&0\\\ 1&0&1&0&0\\\ 0&1&0&0&1\\\ 0&1&0&1&0\\\ \end{pmatrix}$ (3)

Let $X_{1}$ and $X_{2}$ be defined in Equations 2 and 3. The only difference
between the two set of samples is that the binary feature of the first binary
column of $X_{1}$ is encoded on two bits in $X_{2}$ in the first two columns.
It is reasonable to expect that a classifier trained on $X_{1}$ and fed with
$\bar{x}_{1}=\begin{pmatrix}1&1&0&0\end{pmatrix}$ makes the same decision as a
classifier trained on $X_{2}$ and fed with
$\bar{x}_{2}=\begin{pmatrix}0&1&1&0&0\end{pmatrix}$. This is not the case for
a Euclidean nearest-neighbor as the closest vector from $\bar{x}_{1}$ in
$X_{1}$ is $\begin{pmatrix}0&1&0&0\end{pmatrix}$ and the closest ones from
$\bar{x}_{2}$ in $X_{2}$ are $\begin{pmatrix}1&0&1&0&0\end{pmatrix}$,
$\begin{pmatrix}0&1&0&0&1\end{pmatrix}$ and
$\begin{pmatrix}0&1&0&1&0\end{pmatrix}$. On the other hand, the Self-Encoder
returns $\begin{pmatrix}0&1&0&0\end{pmatrix}$ as the most similar to
$\bar{x}_{1}$ and $\begin{pmatrix}1&0&1&0&0\end{pmatrix}$ as the most similar
to $\bar{x}_{2}$.

##  6 Sampling

Under the reasonable assumption that the dimensions of the hidden layers do
not scale with $n$, the time complexity is that of the output layer:
$O(n^{2})$ for training and $O(n)$ for evaluation. The training time
complexity is higher than the $n$-linear time complexity of Euclidean $k$-NN.
However, training is only done once and then the trained model can be used in
$\mathop{}\mathopen{}O\mathopen{}\left(n\right)$ for similarity search.
Moreover, the natural way to implement a Multi-Layer Perceptron (MLP) model is
to use machine learning frameworks such as PyTorch [7] or TensorFlow [1],
which natively support GPU acceleration and highly reduce the computation
time.

To control the memory usage and time complexity of the approach and in order
to improve the performance, we propose to use a simple sampling strategy,
where a random subset of the training set is selected, thus reducing the space
and time complexity of the model. Given a set of training samples
$X=(x^{(i)})_{i\in\\{1,\dots,n\\}}$, sampling generates a new training subsets
of size $s$ by randomly sampling vectors from $X$. The model is then trained
on the new training subset. Experiments show that that the performance of the
sampled model is comparable to the complete one.

##  7 Experiments

To evaluate the quality of the similarity measure learned by the proposed
Self-Encoder, we apply it to a classification task with $k$-NN method . The
Self-Encoder $k$-NN is simply called Self-Encoder for simplicity. It is tested
with and without a hidden layer and in normal and sampling modes. Training
uses early stopping and learning rate decay.

The performances are compared to those of four classification baselines: the
Euclidean $k$-NN with $k=5$, the linear Support Vector Machine (SVM), the one-
vs-all logistic regression and the multi-layer perceptron with one hidden
layer of twenty neurons.

Ten datasets that are recurrent in the machine learning literature were
selected for the experiments. They are all available from the UCI
repository111https://archive.ics.uci.edu/ml/datasets.php and descriptive
figures can be found in Table 1. Among these, the German Credit dataset comes
in two versions: one with numerical features and another one with 13 out of 20
features being categorical. The only pre-processing applied to all the
datasets is the conversion of categorical features into a one-hot encoding.

| # samples |  | Feature
---
dimension
# classes
Breast Cancer | 699 | 9 | 2
Digits | 1,797 | 64 | 10
Ecoli | 336 | 7 | 8
German credit | 1,000 | 24 | 2
| German credit
---
(categorical)
1,000 | 20 | 2
Glass | 214 | 9 | 6
Ionosphere | 351 | 34 | 2
Iris | 150 | 4 | 3
Liver | 345 | 6 | 2
Wine | 178 | 13 | 3
Table 1: Dataset details

The Self-Encoder is implemented in PyTorch, the optimization is done using
Adam [5]. Some hyper-parameters are fixed: the learning rate decay is set to
0.995, the size of the hidden layer is fixed to 20 and for the sampling mode,
the number of visible samples is fixed to 100.

For each classification model and each dataset, the reported metric is the
accuracy. It is measured using 5-fold cross validation. The learning rate to
train each model on each fold is chosen according to a log-uniform
distribution between $0.001$ and $2$, using the Bayesian optimization library
hyperopt222https://hyperopt.github.io/hyperopt. The library is also used to
choose the normalization function between sigmoid and SoftMax. All experiments
can be reproduced using the code provided in the supplementary material.

| $k$-NN | MLP | Logistic | SVM
---|---|---|---|---
Breast cancer | 0.967 | 0.009 | 0.96 | 0.007 | 0.963 | 0.003 | 0.96 | 0.003
Digits | 0.983 | 0.004 | 0.934 | 0.013 | 0.965 | 0.004 | 0.983 | 0.002
Ecoli | 0.869 | 0.021 | 0.762 | 0.016 | 0.759 | 0.029 | 0.798 | 0.019
German credit | 0.684 | 0.012 | 0.741 | 0.026 | 0.738 | 0.012 | 0.754 | 0.012
German credit cat | 0.716 | 0.027 | 0.717 | 0.023 | 0.733 | 0.027 | 0.736 | 0.023
Glass | 0.616 | 0.066 | 0.439 | 0.099 | 0.579 | 0.046 | 0.616 | 0.06
Ionosphere | 0.835 | 0.032 | 0.903 | 0.033 | 0.855 | 0.039 | 0.838 | 0.03
Iris | 0.953 | 0.05 | 0.913 | 0.129 | 0.967 | 0.03 | 0.987 | 0.027
Liver | 0.693 | 0.025 | 0.609 | 0.098 | 0.684 | 0.063 | 0.687 | 0.057
Wine | 0.714 | 0.046 | 0.321 | 0.1 | 0.966 | 0.021 | 0.977 | 0.011
(a) Classification accuracy of the baseline methods.

| SE | SE hidden | Best Baseline
---|---|---|---
Breast cancer | 0.969 | 0.007 | 0.976 | 0.011 | 0.967 | 0.009
Digits | 0.983 | 0.012 | 0.655 | 0.117 | 0.983 | 0.004
Ecoli | 0.911 | 0.028 | 0.899 | 0.021 | 0.869 | 0.021
German credit | 0.768 | 0.011 | 0.749 | 0.007 | 0.754 | 0.012
German credit cat | 0.768 | 0.025 | 0.756 | 0.02 | 0.736 | 0.023
Glass | 0.794 | 0.033 | 0.766 | 0.041 | 0.616 | 0.06
Ionosphere | 0.94 | 0.023 | 0.963 | 0.031 | 0.903 | 0.033
Iris | 1.0 | 0.0 | 0.993 | 0.013 | 0.987 | 0.027
Liver | 0.759 | 0.02 | 0.768 | 0.016 | 0.693 | 0.025
Wine | 0.736 | 0.071 | 0.871 | 0.116 | 0.977 | 0.011
(b) Classification accuracy of the self-encoder. The score of the best
baseline is recalled for comparison.

| SE | SE hidden | Best Baseline
---|---|---|---
Breast cancer | 0.960 | 0.006 | 0.970 | 0.008 | 0.967 | 0.009
Digits | 0.853 | 0.019 | 0.528 | 0.161 | 0.983 | 0.004
Ecoli | 0.878 | 0.036 | 0.899 | 0.025 | 0.869 | 0.021
German credit | 0.741 | 0.012 | 0.750 | 0.008 | 0.754 | 0.012
German credit cat | 0.742 | 0.011 | 0.750 | 0.011 | 0.736 | 0.023
Glass | 0.766 | 0.03 | 0.771 | 0.027 | 0.616 | 0.06
Ionosphere | 0.897 | 0.028 | 0.946 | 0.017 | 0.903 | 0.033
Iris | 1.0 | 0.0 | 0.993 | 0.013 | 0.987 | 0.027
Liver | 0.762 | 0.024 | 0.760 | 0.007 | 0.693 | 0.025
Wine | 0.702 | 0.04 | 0.848 | 0.073 | 0.977 | 0.011
(c) Classification accuracies of the Self-Encoder in sampling mode (with 100
visible samples chosen uniformly at random). The score of the best baseline is
recalled for comparison.

Table 2: Classification accuracies for the Self-Encoder along with other
baselines on ten usual datasets. Performances are measured using a 5-fold
cross validation mechanism. Best results are in bold. German credit cat is the
version of the German credit dataset with categorical features.

All results are reported in Table 2. An obvious conclusion is that in both
normal and sampling settings, the Self-Encoder performs better than the
baselines. A second conclusion is that the Self-Encoder in the sampling
framework performs comparably or better than the other baselines in the normal
setting. This is reassuring because the Self-Encoder might need to be applied
in sampling mode while other lighter models have access to all the samples.

To measure the impact of the invariance feature, the performance of the Self-
Encoder has also been compared to the performance of the Euclidean $k$-NN
algorithm on normalized numerical datasets. The results are shown in Table 3.
As expected, normalization improves the score of the Euclidean $k$-NN in most
cases but not up to the score of the Self-Encoder.

|  |  $k$-NN
---
normalized
$k$-NN | SE
Breast cancer | 0.960 | 0.967 | 0.976
Digits | 0.970 | 0.983 | 0.983
Ecoli | 0.866 | 0.869 | 0.911
German | 0.693 | 0.684 | 0.768
Glass | 0.645 | 0.616 | 0.794
Ionosphere | 0.835 | 0.835 | 0.963
Iris | 0.960 | 0.953 | 1.0
Liver | 0.609 | 0.693 | 0.768
Wine | 0.966 | 0.714 | 0.977
Table 3: Classification accuracy of Euclidean $k$-NN on normalized datasets
and Euclidean $k$-NN and Self-Encoder on raw datasets. This comparison is only
reported on numerical datasets. Two best scores for each dataset are in bold.

##  8 Conclusion

In conclusion, the Self-Encoder is an unsupervised method that learns a
similarity measure specific to the training data that can be used for
downstream tasks such as classification, its primary objective being to
separate samples from one another.

## References

  * [1] Abadi, M., et. al.: TensorFlow: Large-scale machine learning on heterogeneous systems (2015), software available from tensorflow.org
  * [2] Borisov, V., Leemann, T., Seßler, K., Haug, J., Pawelczyk, M., Kasneci, G.: Deep neural networks and tabular data: A survey. IEEE Transactions on Neural Networks and Learning Systems (2022)
  * [3] Cortes, C., Vapnik, V.: Support-vector networks. Machine learning 20(3), 273–297 (1995)
  * [4] Cover, T., Hart, P.: Nearest neighbor pattern classification. IEEE Transactions on Information Theory (1967)
  * [5] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. In: Bengio, Y., LeCun, Y. (eds.) 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (2015), http://arxiv.org/abs/1412.6980
  * [6] Kingma, D.P., Welling, M.: Auto-encoding variational bayes. stat 1050, 1 (2014)
  * [7] Paszke, A., et. al.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc. (2019)
  * [8] Vincent, P., Larochelle, H., Lajoie, I., Bengio, Y., Manzagol, P.A.: Stacked Denoising Autoencoders: Learning Useful Representations in a Deep Network with a Local Denoising Criterion. The Journal of Machine Learning Research 11, 3371–3408 (Dec 2010)