# Adaptive Inference:
Theoretical Limits and Unexplored Opportunities

Soheil Hor, Ying Qian, Mert Pilanci, Amin Arbabian
Department of Electrical Engineering
Stanford University
Stanford, CA 94305
<EMAIL_ADDRESS>

###### Abstract

This paper introduces the first theoretical framework for quantifying the
efficiency and performance gain opportunity size of adaptive inference
algorithms. We provide new approximate and exact bounds for the achievable
efficiency and performance gains, supported by empirical evidence
demonstrating the potential for 10-100x efficiency improvements in both
Computer Vision and Natural Language Processing tasks without incurring any
performance penalties. Additionally, we offer insights on improving achievable
efficiency gains through the optimal selection and design of adaptive
inference state spaces.

##  1 Introduction

In recent years, neural networks have achieved human-level performance across
various domains, ranging from image classification using vision transformers
to intricate natural language processing tasks handled by Large Language
Models. However, this notable improvement in performance comes with the caveat
of training progressively larger models. The current high-performing vision
transformers and large language models can only be effectively deployed on
large cloud data-centers, leading to significant economical costs and
environmental implications in terms of carbon footprint and energy consumption
Anthony et al. (2020); McDonald et al. (2022); Desislavov et al. (2023).

As 80-90% of the current cloud workload constitutes of inference tasks
McDonald et al. (2022); Freund (2019), and this percentage is projected to
increase further as models achieve peak performance and maturity, the demand
for efficient inference has transitioned from a mere consideration to an
immediate necessity Samsi et al. (2023); Desislavov et al. (2023). This
urgency is particularly heightened in non-cloud (edge) applications, where
there is a demand for low latency execution of models on devices with limited
memory, compute, and power resourcesXu et al. (2018); Li et al. (2019);
Daghero et al. (2021).

The advent of network compression techniques, such as pruning and quantization
Han et al. (2023), marked a significant stride in efficient inference. The
initial achievements in this area paved the way for subsequent developments
such as resource-aware neural architecture search, model distillation, and
low-rank decomposition techniques, all aimed at enhancing the performance and
efficiency of neural networks, either during training or as post-processing
steps Xu & McAuley (2023); Li et al. (2023b). However, such efficiency
advancements have reached a plateau, necessitating fundamentally new
techniques that extend beyond the design space of conventional offline neural
network optimization methods.

One such technique is adaptive inference, founded on the intuition that
accurate inference for different instances naturally varies in difficulty. For
simpler instances, a smaller neural network might perform as accurately as a
larger one. Hence, an adaptive neural network (or an adaptive ensemble of
networks) capable of dynamically adjusting its size based on the difficulty of
the input instance can prove to be more efficient and, in some cases, even
more accurate than an equivalent “static” model.

Adaptive inference methods, as commonly explored in the literature, make use
of networks with dynamically tunable size and complexity Han et al. (2021).
Such networks are often adapted through techniques such as early exiting
Laskaridis et al. (2021); Ilhan et al. (2023); Yang et al. (2020) or through
adaptive selection and mixture of experts Meng et al. (2020); Li et al.
(2023a); Jawahar & Mukherjee (2023); Chen et al. (2023). One illustrative
examples in context of Computer Vision (CV) is AR-Net Meng et al. (2020),
which showcases the use of a simple policy network during the inference phase
to adaptively select between pre-trained classifiers with varying sizes and
resolutions, achieving efficient video-based activity classification. Another
example in context of Natural Language Processing (NLP) is a study by Rotem et
al. (2023) comparing the performance and efficiency of both multi-model and
early exiting approaches on large language models like BERT.

The relatively restricted adoption and application of adaptive inference
comparing to static neural network compression methods can be attributed to
the ad-hoc nature of existing methods, lacking a comprehensive framework for
designing adaptive inference data pipelines or gaining insight into the
potential and limitations of adaptive inference in specific tasks and
applications.

This paper is the first to establish a theoretical foundation for analyzing
adaptive inference methods and quantifying achievable efficiency and
performance gains for general inference tasks. The proposed framework aims to
bridge the gap between current ad-hoc methods and a more systematic approach
to understanding and leveraging adaptive inference.

Our contributions encompass:

  * •

A novel theoretical framework for the analysis of adaptive inference systems,
achieved through the definition of conceptual adaptive Oracles.

  * •

Introduction of both approximate and exact bounds on the performance and
efficiency gains achievable by an adaptive agent. This marks the inception of
new quantitative measures for adaptation potential.

  * •

Empirical findings showcasing adaptation potential and limits of models,
demonstrated in the realms of both Computer Vision (Image Classification) and
Natural Language Processing (Natural Language Inference).

  * •

Design considerations and recommendations for maximizing efficiency and
performance gain potential of neural networks.

  * •

Ground truth “adaptation labels” for optimal adaptive inference, presented for
two datasets and four neural networks in the context of image classification
on ImageNet and CommonSense NLI on HellaSwag.

##  2 A General Framework for Adaptive Inference

As previously discussed, adaptive inference serves as a broad term applicable
to a diverse range of systems, applications, and methodologies. Nevertheless,
the majority of adaptive systems can be effectively abstracted at the highest
level as a state machine. This abstraction enables the separation of
considerations and constraints imposed by the “adaptation state space” of a
system from the performance of a specific “agent” responsible for guiding the
system transitions between states.

In this section, we begin by defining an adaptation state space within the
context of a classification task. Subsequently, we present precise definitions
and equations that explore the theoretical limits of performance and
efficiency achievable by adaptive inference agents. This is achieved through
the analysis and definition of ideal “Adaptive Oracle”s.

###  2.1 Adaptation State Space for Classification

In classification-based adaptive inference, we often see the adaptation state
space defined either using an ensemble of backbone classifiers (like in AR-Net
Meng et al. (2020), switching between different classifiers from EfficientNet
family) or a single classifier with adaptable complexity (for instance, RA-Net
Yang et al. (2020) allowing different setups within a single backbone
network).

For simplicity, we conceptualize both scenarios by representing them as a
discrete set of $N$ backbone classifiers applied to a dataset $X$. These
classifiers constitute a state space, defined as:

$\\{S_{i}\\}\,\mathrm{for}\,i\in\\{1,2,3,...,N\\}.$ (1)

Suppose further that these classifiers are ranked based on the amount of
resources they consume in an increasing order. In other words, let the set
$\\{R_{i}\\}$ represent the resource consumption of each classifier in
$\\{S_{i}\\}$, then:

$R_{1}\leq R_{2}\leq...\leq R_{N}.\\\ $ (2)

Then, let the set $\\{A_{i}\\}$ represent the test accuracies of each
classifier within the set $\\{S_{i}\\}$. Hence, for each classifier or state
$S_{i}$, there exists a corresponding pair $(R_{i},A_{i})$, representing both
the resource consumption and accuracy of that specific classifier. In
practical systems, larger (and more resource-intensive) models are typically
employed only if they deliver better or equal performance compared to smaller
models. Consequently, we assume that the model accuracies follow an increasing
order

$A_{1}\leq A_{2}\leq...\leq A_{N}.\\\ $ (3)

Given the definition of the adaptation state space, an adaptive agent aims to
identify an optimal strategy that maximizes performance ($A$) and minimizes
resource consumption ($R$) by selecting the optimal state ($S$) for a given
input $x$.

We establish the performance and efficiency bounds attainable by any adaptive
agent through the concept of an “Adaptive Oracle,” as defined in the
subsequent section.

###  2.2 Ideal Adaptive Agents: Adaptive Oracles

An adaptive Oracle is characterized as an agent possessing simultaneous
knowledge of both resource consumption ($R$) and per-instance accuracy
(correctness) across all models for each instance $x$. Consequently, it can
select the model with the least resource consumption while still achieving the
highest possible accuracy (allowed by the adaptation state space) for every
classified instance.

It’s important to clarify that, in this definition, the adaptive Oracle is not
assumed to have any knowledge beyond what is provided within the defined
adaptation space (such as the ground truth label for each instance).
Consequently, it cannot in general ensure correct predictions (or $100\%$
accuracy) for every instance. This stands in contrast to the more commonly
used “label Oracle”s found in the literature.

As a conceptual example, consider a 2-state adaptation problem with two
backbone classifiers of different sizes:

Figure 1: Confusion matrix for a conceptual 2-state classification task.
Resource consumption of the Adaptive Oracle is only a function of the $P(IA)$

Consider a larger classifier characterized by $S_{L}=(R_{L},A_{L})$ and a
smaller model characterized by $S_{s}=(R_{s},A_{s})$. In Figure 1, there are
only four cases to be considered based on the per-instance accuracy of each
classifier. Given that opting for more resources (selecting a larger model) is
justified only if it leads to a better relative accuracy, it can be argued
that an adaptive Oracle would choose the larger model for a specific instance
only when the smaller model is inaccurate, i.e., incorrect, while the larger
model is accurate, i.e., correct, (as depicted in the IA case in Figure 1).

####  2.2.1 General formulation

Building upon the insights from this straightforward 2-state scenario, we
present the following general definition for an adaptive Oracle:

######  Definition 2.1.

Given a state space $\\{S_{i}\\}$, its corresponding $\\{R_{i}\\}$ and
$\\{A_{i}\\}$, and a dataset $X$, the adaptive Oracle is defined as an
adaptive agent that implements the following strategy:

$R_{oracle}(x)=\begin{cases}\min_{i}(R_{i})&s.t.\>Y_{i}(x)=Y_{GT}(x)\\\
R_{1}&O.W.\\\ \end{cases}$ (4)

for all $x\in X$,

Where $Y_{i}(x)$ is the predicted label from model $i$ on instance $x$, and
$Y_{GT}(x)$ is the ground truth label of instance $x$. The expected resource
consumption and accuracy achieved by such an Oracle are $R_{oracle}$ and
$A_{oracle}$ calculated as111For a detailed proof of each equation please see
the Appendix:

$\displaystyle R_{oracle}$ $\displaystyle=R_{1}(1-P(e_{1})+P(e_{N}))+$
$\displaystyle\sum_{i=2}^{N}R_{i}[P(e_{i-1})-P(e_{i})],$ $\displaystyle
A_{oracle}$ $\displaystyle=1-P(e_{N}),$ (5)

In which $P(e_{i})$ is the probability of event $e_{i}$ defined as:

$e_{i}=\\\ \\{Y_{1}\neq Y_{GT}\cap Y_{2}\neq Y_{GT}\cap\cdots Y_{i-1}\neq
Y_{GT}\cap Y_{i}\neq Y_{GT}\\},\\\ $

This can be interpreted as the event in which the $i$ smallest models fail to
classify an instance correctly.

To get a better intuition on the equations above one can use the Bayes rule,
and the fact that $A_{i}=1-P(Y_{i}\neq Y_{GT})$, to write each $P(e_{i})$ as:

$P(e_{i})=\begin{cases}\alpha_{i}(1-A_{i}),&i>1\\\ (1-A_{1})&i=1\\\
\end{cases}$ (6)

In which $\alpha_{i}$ is defined for $i>1$ and can be written as:

$\alpha_{i}=P(Y_{1}\neq Y_{GT}\cap Y_{2}\neq Y_{GT}\cap\cdots Y_{i-1}\neq
Y_{GT}|Y_{i}\neq Y_{GT}).\\\ $

Intuitively, larger $\alpha_{i}$ values (approaching 1) indicate states where
the errors of larger models are inherently challenging to resolve using any of
the smaller models. Conversely, a smaller $\alpha_{i}$ represents the scenario
in which an ensemble of smaller models are capable of resolving some or all of
the classification errors of a larger model.

Using this definition Equation (5) can be reformulated as:

$\displaystyle R_{oracle}=R_{1}+(R_{2}-R_{1})(1-A_{1})+$
$\displaystyle\sum_{i=3}^{N}[\alpha_{i-1}(R_{i}-R_{i-1})(1-A_{i-1})]-\alpha_{N}(R_{N}-R_{1})(1-A_{N})$
$\displaystyle A_{oracle}=1-\alpha_{N}[1-A_{N}],$ (7)

The resource consumption and accuracy of an adaptive Oracle calculated using
this equation can serve as an upper bound on the performance and efficiency
achievable by any adaptive agent applied on the same adaptation state space.
Calculating this upper bound however relies on knowledge about the backbone
state space, characterized by $R_{i}$’s and $A_{i}$’s, along with the hidden
term $\alpha_{i}$.

For pre-trained off-the-shelf backbone models, $R_{i}$ and $A_{i}$ values are
typically readily available since they can be calculated separately for each
classifier.

On the other hand, $\alpha_{i}$ captures the cross-dependencies among the
entire set of backbone models, a detail often not reported for static off-the-
shelf models. Obtaining an empirical estimate of $\alpha_{i}$ while not
impossible, necessitates access to both a representative validation set and a
comprehensive set of candidate backbone models. This poses a challenge,
especially when constructing an adaptive inference pipeline from scratch or
when the backbone models undergo frequent retraining to uphold performance
amidst real-world data distribution shifts.

Fortunately, for classifiers with similar structure that are trained using the
same training set, variations of $\alpha_{i}$s between states can be
relatively small (as shown in Figure 2). This allows for calculation of an
approximate performance and efficiency bound for adaptive Oracles without the
need for calculating $\alpha_{i}$s for the entire adaptation space. This is
the motivation for the constant-$\alpha$ formulas investigated in the next
section.

####  2.2.2 Constant-$\alpha$ Approximation

As previously discussed, $\alpha_{i}$ serves as a measure of the probability
that a large model making a classification error leads to errors in all of the
smaller models.

Intuitively, if the classifiers forming the adaptation state space were
statistically independent, one would expect $\alpha_{i}$ to quickly approach 0
as the number of states increases. This is because as the index $i$ (and
subsequently number of states included in calculation of $\alpha_{i}$)
increases, it becomes increasingly more likely that at least one of the
smaller models predicts a label correctly by chance.

However, the expected decrease in the rate of $\alpha_{i}$ is less pronounced
when models forming the state space are not completely independent. In such
cases, larger models are anticipated to correctly classify most, if not all,
of the samples that are correctly classified by the smaller models. This
tendency is commonly observed in backbone classifier families with similar
network structures, as demonstrated in Figure 2.

Building on this intuition, one straightforward approach is to assume
$\alpha_{i}$ to be constant and independent of the state index
($\alpha_{i}=\alpha$). In this scenario, Equation (7) can be modified as:

$\displaystyle R_{oracle}=R_{1}+(R_{2}-R_{1})(1-A_{1})+$
$\displaystyle\alpha\left[\sum_{i=3}^{N}[(R_{i}-R_{i-1})(1-A_{i-1})]-(R_{N}-R_{1})(1-A_{N})\right]$
$\displaystyle A_{oracle}=1-\alpha[1-A_{N}].$ (8)

This equation reveals that under the constant-$\alpha$ assumption, the
relationship between $R_{oracle}$ and $A_{oracle}$ is a line connecting a very
optimistic operating point with $A_{oracle}=1$ and
$R_{oracle}=R_{1}+(R_{2}-R_{1})(1-A_{1})$ (associated with $\alpha=0$) to a
more realistic “conservative” bound with an accuracy of $A_{oracle}=A_{N}$
corresponding to $\alpha=1$.

For $\alpha=1$ we have:

$\displaystyle{R_{oracle}}$
$\displaystyle=R_{1}+\sum_{i=2}^{N}(R_{i}-R_{1})(A_{i}-A_{i-1}),$
$\displaystyle A_{oracle}$ $\displaystyle=A_{N}.$ (9)

This equation serves as a conservative estimate for the performance and
efficiency gains achievable by an adaptive agent, requiring only knowledge
about the $R_{i}$ and $A_{i}$ values for each state that are typically readily
available for common well-known classifiers.

Figure 2: Empirical Measurements of $\alpha_{i}$ for different tasks and
models. $\alpha_{i}$ remains relatively constant for models with similar
architecture. Table 1: Estimated Adaptation Opportunity Bounds

|  |  | Baseline State Space |  | Conservative Bound
---
($\alpha=1$)
| Optimistic Bound
---
($\alpha=\alpha_{min}$)
|  |  |  | Performance
---
Baseline
| Efficiency
---
Baseline
| Efficiency Gain
---
Opportunity
| $\alpha$
---
Estimate
| Performance Gain
---
Opportunity
| Efficiency Gain
---
Opportunity
Context | Edge/Cloud | Model Family | $A_{N}$ |  | $R_{1}$
---
(GFLOPs)
| $R_{N}$
---
(GFLOPs)
| $R_{oracle}$
---
(GFLOPs)
| $\Delta R$
---
(GFLOPs)
$R_{ratio}$ | $\alpha_{min}$ |  $A_{oracle}$ | $\Delta A$ |  | $R_{oracle}$
---
(GFLOPs)
| $\Delta R$
---
(GFLOPs)
$R_{ratio}$
| Edge | EfficientNet | 83.95% | 0.39 | 37.75 | 0.60 | 37.15 | 63.43x | 0.58 | 90.67% | +6.72% | 0.54 | 37.21 | 70.26x
| Cloud | ViT | 88.60% | 4.41 |  1,016.72 | 23.21 | 993.51 | 43.80x | 0.52 | 94.00% | +5.66% | 15.56 | 1,001.16 | 65.33x
CV (ImageNet) | All | SOTA | 90.88% | 0.04 | 2,586.00 | 21.24 | 2,564.76 | 121.77x | -
| Edge | Pythia | 67.08% | 78.96 |  2,910.00 | 304.42 | 2,605.58 | 9.56x | 0.88 | 71.13% | +4.05% | 286.14 | 2,623.86 | 10.17x
| Cloud | Llama-2 | 83.79% | 1,670.00 | 17,570.00 | 2,423.36 | 15,146.64 | 7.25x | 0.90 | 85.43% | +1.64% | 2,385.65 | 15,184.35 | 7.36x
NLP (HellaSwag) | All | SOTA | 90.96% | 1.90 | 1,352,640.00 | 16,694.40 | 1,335,945.62 | 81.02x | -

(a) Imagenet

(b) HellaSwag

Figure 3: Benchmarks and corresponding $\alpha=1$ bounds. The shaded area
shows the space of operation points achievable using adaptive inference
techniques. The state of the art (SOTA) baseline is used as a proxy for the
inherent efficiency versus performance trade-off of each task.

##  3 Experiments

In the preceding section, we introduced both exact and approximate bounds
aimed at evaluating the achievable efficiency and performance gains of an
adaptive agent. This section delves into the practical applications of these
bounds on real-world off-the-shelf neural networks.

Our exploration of adaptation spaces focuses on two distinct inference tasks:
image classification on ImageNet Russakovsky et al. (2015) and Natural
Language Inference on HellaSwag Zellers et al. (2019). For each of these
tasks, we assessed models tailored for efficiency-critical applications (e.g.,
image classification at the edge using Efficient-Net) as well as performance-
critical applications (e.g., state-of-the-art Llama-2 LLM models deployed on
the cloud). Within each state space formed by these models, we present the
maximum accuracy achievable by an adaptive agent ($A_{oracle}$) along with the
minimum resource consumption required to attain such accuracy ($R_{oracle}$).

The estimated $R_{oracle}$ is then employed to derive two quantitative
measures of efficiency gain: $\Delta R=R_{N}-R_{oracle}$ and
$R_{ratio}=R_{N}/R_{oracle}$, and ($\Delta A=A_{oracle}-A_{N}$) as a measure
of performance gain as detailed in Table 1.

###  3.1 Image Classification Benchmark: ImageNet

ImageNet stands out as one of the most renowned and demanding datasets for
image classification, featuring high-resolution images spanning 1000 classes
of diverse objects. Off-the-shelf classifiers trained on ImageNet range from
compact and efficient models tailored for resource-limited at-the-edge
inference (such as Efficient-Net Tan & Le (2019)) to high-performance models
typically deployed on the cloud e.g., Vision Transformers Dosovitskiy et al.
(2020)).

In Figure 3(a), we present the Performance (Accuracy) versus Resource
Consumption (GFLOPs) profiles for two of the prominent pre-trained classifiers
on ImageNet, encompassing a broad spectrum of resource requirements and
performance capabilities. Additionally, we have calculated the GFLOPs versus
accuracy envelope of the state-of-the-art on ImageNet, serving as a proxy for
the global adaptation potential of ImageNet (datapoints sourced from the
papers with code leaderboard Paperswithcode (2024)).

Utilizing Equation (9), in conjunction with GFLOPs and accuracy metrics
reported in literature for each model, we derived a conservative estimate of
the achievable adaptation bounds for each model, as illustrated in Figure 3(a)
and summarized in Table 1.

For ImageNet models with a large number of states (e.g., EfficientNet, ViT),
even the conservative assumption of $\alpha=1$ suggest a substantial
efficiency improvement potential, in orders of 43-63x. Moreover, the analysis
indicates potential for efficiency gains exceeding 121x using the entire
state-of-the-art envelope of ImageNet.

###  3.2 Natural Language Inference Benchmark: HellaSwag

HellaSwag serves as a widely adopted benchmark in the domain of Commonsense
Natural Language Inference. Comprising over 10,000 sets of sentences, each
with four potential endings, this dataset tasks language models with selecting
the most probable conclusion for a given sentence. The dataset is crafted
specifically to necessitate commonsense reasoning based on contextual cues in
addition to the words within a sentence.

For this particular task, we chose a large language model (Llama-2 Touvron et
al. (2023)), typically deployed on cloud infrastructure, alongside a more
compact language model (Pythia Biderman et al. (2023)), crafted for deployment
on resource-limited systems. Additionally, we incorporated the GFLOPs vs
Accuracy envelope for state-of-the-art models from the HuggingFace LLM
leaderboard Hugginface (2024); Gao et al. (2023) as a representation of the
overarching Resource Consumption vs Accuracy trade-offs associated with
HellaSwag. The inference cost of each model (measured in GFLOPs) was computed
assuming a batch size of 1 and a maximum sequence length of 128, utilizing Ye
(2024) for calculations.

As depicted in Figure 3(b) and detailed in Table 1, state-of-the-art language
models show great potential for substantial efficiency improvements through
adaptive inference. Specifically, the smaller language model (Pythia) boasts a
conservative bound of over 9x in achievable adaptation efficiency gains.
Conversely, the larger language model demonstrates the potential for
efficiency gains exceeding 7x, a noteworthy accomplishment given the
considerable size and scale of such models, resulting in a relative resource
consumption reduction of more than 15 TFLOPs. Notably, these accomplishments
are surpassed only by the global adaptation potential of the state-of-the-art
models on HellaSwag, indicating over 81x potential improvements in efficiency
without sacrificing performance.

###  3.3 Empirical (Exact) Adaptation Bounds

The $\alpha=1$ bounds discussed in the preceding section serve as a
conservative estimate on the adaptation potential of various models and
datasets. More accurate estimates of the efficiency and performance achievable
by an adaptive agent can be obtained by considering the hidden dependencies
between models ($\alpha_{i}$s) within a specific adaptation state space.

Figure 2 presents the experimentally determined $\alpha_{i}$ values for four
state spaces highlighted in the earlier sections showing that $\alpha_{i}$
values are relatively constant for models with similar structures and
training. This observation motivates obtaining a more accurate estimate of
$\alpha$ using the constant-$\alpha$ equations provided in previous sections
(Equation (8)).

Table 1 showcases the minimum $\alpha$ values measured for each of the four
classifiers and their corresponding efficiency and performance gain
opportunity bounds. It’s crucial to highlight that, unlike the $\alpha=1$
bounds, which assumed that no accuracy gain are achievable through adaptation,
$constant-\alpha$ estimates can be employed to calculate both efficiency and
performance gains. For instance, within the EfficientNet family of
classifiers, usin $\alpha=\alpha_{min}$ to get an optimistic estimate on the
performance of an adaptive agent results in an estimated accuracy of
$90.67\%$—over $6.72\%$ more accurate comparing to the largest model in the
state space.

The reported resource consumption values represent the minimum resources
required to achieve the performance bounds, indicating that simultaneous
improvements in efficiency and performance are attainable for all models based
on the $\alpha$ values. For example, the EfficientNet and ViT families (with
smaller $\alpha$ values) can achieve accuracy gains of over $5-6\%$ while
realizing efficiency gains of over 70x and 65x, respectively. On the other
hand, tasks related to the HellaSwag dataset exhibit larger $\alpha$ values,
resulting in relatively smaller accuracy gains ($4.05\%$ and $1.64\%$ for
Pythia and LLama-2, respectively) but still showcasing efficiency gains of
over 7-10x for both studied language models.

(a) EfficientNet

(b) ViT

(c) Pythia

(d) Llama-2

Figure 4: Proposed $constant-\alpha$ bounds versus adaptation gain opportunity
space empirically measured for an adaptive Oracle. The shaded area shows the
space of operation points potentially achievable by an adaptive inference
algorithm.

Figure 4 illustrates the comparison between the approximate bounds and
empirical measurements of an ideal Oracle’s efficiency and performance across
the adaptation space. The visualization highlights that the proposed
conservative and optimistic estimates can serve as accurate bounds on the
actual operating points achievable by an adaptive Oracle (Opportunity Space),
particularly for classifiers with larger $\alpha$ values.

##  4 State Space Design Considerations

In the preceding section, we established the utility of Equation (9) as a
measure of adaptation potential for off-the-shelf state spaces. In this
section, we aim to provide intuitions and design guidelines for enhancing the
adaptation potential of a given adaptation state space.

###  4.1 Effect of State Dynamic Range

Equation (9) reveals a fundamental observation: $R_{oracle}$ is always bounded
below by $R_{1}$, while $A_{oracle}$ is directly tied to $A_{N}$.
Consequently, careful selection of the smallest and largest models in the
adaptation space, defining the state dynamic range, is of great importance.
Utilizing Equation (10) as a simplified form of Equation (9) for 2 states
underscores that adapting with state spaces featuring a broader dynamic range
incurs higher efficiency costs but yields more substantial performance
benefits.

$\displaystyle A_{oracle}$ $\displaystyle=A_{N},$ $\displaystyle R_{oracle}$
$\displaystyle=R_{1}+(A_{N}-A_{1})(R_{N}-R_{1}).$ (10)

####  4.1.1 Maximizing Efficiency: Optimum $S_{1}$

Given the equation above, the instinctive choice might be to optimize $S_{1}$
for maximum efficiency (rather than necessarily performance). However, there
exists a minimum accuracy threshold for even the smallest state. A practical
method to evaluate whether the accuracy of the smallest state justifies its
resource consumption is to confirm that the corresponding $R_{oracle}$ is
lower than a scenario in which the smallest state provides a random guess on
the target class without consuming any resources. Guided by this insight, a
design criterion for the resource consumption of the smallest state ($R_{1}$)
can be expressed as :

$R_{1}<\frac{A_{1}-\frac{1}{C}}{(1-A_{N}+A_{1})}R_{N}$ (11)

In which $C$ is is number of classes in the classification task.

####  4.1.2 Maximizing Performance: Optimum $S_{N}$

Unlike $S_{1}$, the accuracy of the Oracle directly hinges on the accuracy of
the most accurate state ($S_{N}$). Hence, it is logical to craft the largest
model in the adaptation space for performance, prioritizing it over
efficiency.

As demonstrated in Section 3, attaining higher accuracy for most real-world
models often incurs an exponential increase in resource consumption costs.
Drawing from Equation (10), the efficiency of the Oracle can exhibit a strong
correlation with $R_{N}$, particularly in state spaces with a larger dynamic
range. This implies that the Oracle’s resource consumption also grows
exponentially relative to its performance. In the subsequent section, we
illustrate that increasing the number of intermediate states (increasing state
granularity) can alleviate this issue, particularly in scenarios where $R_{N}$
is exceptionally large.

###  4.2 Effect of State Space Size and Granularity

Through a simple rearrangement of terms, Equation (9) can be expressed as:

$\displaystyle R_{oracle}$ $\displaystyle=R_{1}+(A_{N}-A_{1})(R_{N}-R_{1})-$
$\displaystyle\sum_{i=2}^{N}(A_{i}-A_{i-1})(R_{N}-R_{i})$ (12)

The first two terms in this equation represent $R_{oracle}$ calculated using
the state dynamic range. However, the third term
$\sum_{i=2}^{N}(A_{i}-A_{i-1})(R_{N}-R_{i})$ is a positive sum that reduces
$R_{oracle}$ as each additional state is added, constantly improving the
efficiency of the Oracle.

For real-world adaptive systems, expanding the number of states often
accompanies increased complexity. Therefore, it is crucial to discern the
minimum number of states that result in a sufficient adaptation gain. The
following section provides an example demonstrating how Equation (12) can be
applied to select a small but effective adaptation space.

Figure 5: Optimum discrete state spaces of different size and corresponding
$R_{ratio}$ for the ImageNet SOTA. The red dot shows the state with the most
utility relative to the immediately smaller state space.

####  4.2.1 Optimum Choice of States

Upon revisiting Equation (12), it becomes evident that the efficiency gain
resulting from iteratively growing a state space depends solely on the
resource consumption of each state ($R_{i}$) at each step and its accuracy
relative to the most similar states existing in the state space from the
previous steps ($A_{i}-A_{i-1}$). Therefore, the utility of each state for all
state space sizes can be calculated in polynomial time. Figure 5 illustrates
the results obtained using such an algorithm to determine the optimal
adaptation spaces of various sizes on the ImageNet SOTA.

Figure 5 is evidence that through an optimal design of the state space, a
remarkably high adaptation potential can be achieved even within relatively
small state spaces. The corresponding efficiency gains from optimally
expanding the adaptation state space are highlighted in Figure 6. Notably, for
ImageNet SOTA, it is possible to realize over a 100x efficiency gain
(equivalent to 90% of the efficiency gain of the largest discrete state space)
using only 7 states, and efficiency gains of over 80x using only 4 states.

Figure 6: Efficiency gains achievable for discrete state spaces of different
sizes. For the ImageNet SOTA it is possible to achieve 90% of the maximum
efficiency gain (over 100x $R_{ratio}$) utilizing only 7 states.

####  4.2.2 Optimum Number of States

The efficiency gain potential figures presented in the preceding sections
imply that expanding the state space size exhibits diminishing returns in
terms of efficiency gains. Nevertheless, given that Oracles with larger state
spaces consistently outperform those with smaller state spaces, considering
the concept of an Oracle with access to an infinite number of states
($N\rightarrow\infty$ forming a continuous state space) becomes valuable in
theoretically quantifying the achievable efficiency gains for a specific
dataset or benchmark.

Let $R_{h}=\lim_{N\rightarrow\infty}R_{N}$ and $A_{h}$ be defined as accuracy
of the largest state in the continuous adaptation space. The continuous
reformulation of Equation (12) can be then derived as:

$R_{oracle}=R_{1}+A_{h}(R_{h}-R_{1})-\int_{R_{1}}^{R_{h}}A(R)\,dR,$ (13)

where $A(R)$ represents the curve depicting the relationship between accuracy
and resource consumption in the continuous state space.

As an example, we utilized a continuous piece-wise linear approximation of the
SOTA adaptation spaces for ImageNet and HellaSwag to estimate the $\alpha=1$
bound for continuous adaptation. As presented in Table 2, the gains achievable
through continuous adaptation (160.94x and 122.18x for ImageNet and HellaSwag,
respectively) significantly exceed the corresponding figures reported in Table
1 for a discrete adaptation space (121.77x and 81.02x respectively).

Table 2: Conservative ($\alpha=1$) bounds for efficiency gain achievable for
each dataset assuming a continuous adaptation state space

Dataset | $A_{oracle}$ |  | $R_{oracle}$
---
(GFLOPs)
| $\Delta R$
---
(GFLOPs)
$R_{ratio}$
ImageNet | 90.88% | 16.07 | 2,569.93 | 160.94x
HellaSwag | 90.96% | 11,070.49 | 1,341,569.52 | 122.18x

##  5 Limitations and Future Work

The presented work has certain limitations, which will be explored in future
research. This includes extending the proposed framework to tasks beyond
classification, such as regression, and exploring more advanced models for
$\alpha_{i}$ (e.g., linear instead of constant-$\alpha$) to provide tighter
bounds on adaptation potential.

##  6 Conclusion

In this work we introduced a novel theoretical framework, quantifying the
efficiency and performance gains achievable by adaptive inference methods.

Empirical results demonstrated a substantial efficiency gain opportunity,
ranging from 40-70x for models like EfficientNet and ViT (ImageNet), and
exceeding 7x (equivalent to a relative computation saving of over 15 TFLOPs)
for large language models such as Llama-2. Theoretical estimates for datasets
like ImageNet (CV) and HellaSwag (NLP) suggest the potential for achieving
over 80-120x efficiency gain using adaptive inference techniques.

Furthermore, we provided insights and design considerations for further
enhancing the efficiency gain opportunity by carefully designing the
adaptation state space. Empirical results highlight that, for ImageNet,
efficiency improvements on the order of 100x can be attained with adaptation
space sizes of 7 or less. This paper establishes the theoretical foundation
for effective, quantifiable, and systematic design of adaptive inference
methods.

## References

  * Anthony et al. (2020) Anthony, L. F. W., Kanding, B., and Selvan, R.  Carbontracker: Tracking and predicting the carbon footprint of training deep learning models.  _arXiv preprint arXiv:2007.03051_ , 2020.
  * Biderman et al. (2023) Biderman, S., Schoelkopf, H., Anthony, Q. G., Bradley, H., O’Brien, K., Hallahan, E., Khan, M. A., Purohit, S., Prashanth, U. S., Raff, E., et al.  Pythia: A suite for analyzing large language models across training and scaling.  In _International Conference on Machine Learning_ , pp. 2397–2430. PMLR, 2023.
  * Chen et al. (2023) Chen, T., Chen, X., Du, X., Rashwan, A., Yang, F., Chen, H., Wang, Z., and Li, Y.  Adamv-moe: Adaptive multi-task vision mixture-of-experts.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 17346–17357, 2023.
  * Daghero et al. (2021) Daghero, F., Pagliari, D. J., and Poncino, M.  Chapter eight - energy-efficient deep learning inference on edge devices.  In Kim, S. and Deka, G. C. (eds.), _Hardware Accelerator Systems for Artificial Intelligence and Machine Learning_ , volume 122 of _Advances in Computers_ , pp. 247–301. Elsevier, 2021.  doi: https://doi.org/10.1016/bs.adcom.2020.07.002.  URL https://www.sciencedirect.com/science/article/pii/S0065245820300553.
  * Desislavov et al. (2023) Desislavov, R., Martínez-Plumed, F., and Hernández-Orallo, J.  Trends in ai inference energy consumption: Beyond the performance-vs-parameter laws of deep learning.  _Sustainable Computing: Informatics and Systems_ , 38:100857, 2023.
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Freund (2019) Freund, K.  Google cloud doubles down on nvidia gpus for inference, 2019.  _URL https://www. forbes. com/sites/moorinsights/2019/05/09/google-cloud-doubles-down-on-nvidia-gpus-for-inference_ , 2019.
  * Gao et al. (2023) Gao, L., Tow, J., Abbasi, B., Biderman, S., Black, S., DiPofi, A., Foster, C., Golding, L., Hsu, J., Le Noac’h, A., Li, H., McDonell, K., Muennighoff, N., Ociepa, C., Phang, J., Reynolds, L., Schoelkopf, H., Skowron, A., Sutawika, L., Tang, E., Thite, A., Wang, B., Wang, K., and Zou, A.  A framework for few-shot language model evaluation, 12 2023.  URL https://zenodo.org/records/10256836.
  * Han et al. (2023) Han, S., Liu, X., Mao, H., Pu, J., Pedram, A., Horowitz, M. A., and Dally, W. J.  Retrospective: Eie: Efficient inference engine on sparse and compressed neural network.  _arXiv preprint arXiv:2306.09552_ , 2023.
  * Han et al. (2021) Han, Y., Huang, G., Song, S., Yang, L., Wang, H., and Wang, Y.  Dynamic neural networks: A survey.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 44(11):7436–7456, 2021.
  * Hugginface (2024) Hugginface, 2024.  URL https://HuggingFaceH4/open_llm_leaderboard.  [Accessed 08-01-2024].
  * Ilhan et al. (2023) Ilhan, F., Liu, L., Chow, K.-H., Wei, W., Wu, Y., Lee, M., Kompella, R., Latapie, H., and Liu, G.  Eenet: Learning to early exit for adaptive inference.  _arXiv preprint arXiv:2301.07099_ , 2023.
  * Jawahar & Mukherjee (2023) Jawahar, G. and Mukherjee, S. S. e. a.  Automoe: Heterogeneous mixture-of-experts with adaptive computation for efficient neural machine translation.  In _ACL 2023_ , June 2023.
  * Laskaridis et al. (2021) Laskaridis, S., Kouris, A., and Lane, N. D.  Adaptive inference through early-exit networks: Design, challenges and directions.  In _Proceedings of the 5th International Workshop on Embedded and Mobile Deep Learning_ , pp. 1–6, 2021.
  * Li et al. (2019) Li, E., Zeng, L., Zhou, Z., and Chen, X.  Edge ai: On-demand accelerating deep neural network inference via edge computing.  _IEEE Transactions on Wireless Communications_ , 19(1):447–457, 2019.
  * Li et al. (2023a) Li, J., Su, Q., Yang, Y., Jiang, Y., Wang, C., and Xu, H.  Adaptive gating in mixture-of-experts based language models.  _arXiv preprint arXiv:2310.07188_ , 2023a.
  * Li et al. (2023b) Li, Z., Li, H., and Meng, L.  Model compression for deep neural networks: A survey.  _Computers_ , 12(3):60, 2023b.
  * McDonald et al. (2022) McDonald, J., Li, B., Frey, N., Tiwari, D., Gadepally, V., and Samsi, S.  Great power, great responsibility: Recommendations for reducing energy for training language models.  _arXiv preprint arXiv:2205.09646_ , 2022.
  * Meng et al. (2020) Meng, Y., Lin, C.-C., Panda, R., Sattigeri, P., Karlinsky, L., Oliva, A., Saenko, K., and Feris, R.  Ar-net: Adaptive frame resolution for efficient action recognition.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part VII 16_ , pp. 86–104. Springer, 2020.
  * Paperswithcode (2024) Paperswithcode, 2024.  URL https://paperswithcode.com/sota/image-classification-on-imagenet.  [Accessed 6-01-2024].
  * Rotem et al. (2023) Rotem, D., Hassid, M., Mamou, J., and Schwartz, R.  Finding the sweet spot: Analysis and improvement of adaptive inference in low resource settings.  _arXiv preprint arXiv:2306.02307_ , 2023.
  * Russakovsky et al. (2015) Russakovsky, O., Deng, J., Su, H., Krause, J., Satheesh, S., Ma, S., Huang, Z., Karpathy, A., Khosla, A., Bernstein, M., Berg, A. C., and Fei-Fei, L.  ImageNet Large Scale Visual Recognition Challenge.  _International Journal of Computer Vision (IJCV)_ , 115(3):211–252, 2015.  doi: 10.1007/s11263-015-0816-y.
  * Samsi et al. (2023) Samsi, S., Zhao, D., McDonald, J., Li, B., Michaleas, A., Jones, M., Bergeron, W., Kepner, J., Tiwari, D., and Gadepally, V.  From words to watts: Benchmarking the energy costs of large language model inference.  In _2023 IEEE High Performance Extreme Computing Conference (HPEC)_ , pp. 1–9. IEEE, 2023.
  * Tan & Le (2019) Tan, M. and Le, Q.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In _International conference on machine learning_ , pp. 6105–6114. PMLR, 2019.
  * Touvron et al. (2023) Touvron, H., Martin, L., Stone, K., Albert, P., Almahairi, A., Babaei, Y., Bashlykov, N., Batra, S., Bhargava, P., Bhosale, S., et al.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_ , 2023.
  * Xu & McAuley (2023) Xu, C. and McAuley, J.  A survey on model compression and acceleration for pretrained language models.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pp. 10566–10575, 2023.
  * Xu et al. (2018) Xu, X., Ding, Y., Hu, S. X., Niemier, M., Cong, J., Hu, Y., and Shi, Y.  Scaling for edge inference of deep neural networks.  _Nature Electronics_ , 1(4):216–222, 2018.
  * Yang et al. (2020) Yang, L., Han, Y., Chen, X., Song, S., Dai, J., and Huang, G.  Resolution adaptive networks for efficient inference.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 2369–2378, 2020.
  * Ye (2024) Ye, J. X., 2024.  URL https://huggingface.co/spaces/MrYXJ/calculate-model-flops.  [Accessed 08-01-2024].
  * Zellers et al. (2019) Zellers, R., Holtzman, A., Bisk, Y., Farhadi, A., and Choi, Y.  Hellaswag: Can a machine really finish your sentence?  _arXiv preprint arXiv:1905.07830_ , 2019.

##  Appendix A Appendix

###  A.1 Assumptions

General assumptions (applied to all cases):

  1. 1.

The adaptation state space $\\{S_{i}\\}$ is defined for an ensemble of
backbone classifiers.

  2. 2.

The states $S_{i}$’s are ranked in an ascending order based on their resource
consumption $R_{i}$.

  3. 3.

Model accuracies are ranked in an non-decreasing sequence, i.e. $A_{i}\leq
A_{j}$ if $j>i$.

  4. 4.

The adaptive Oracle has knowledge of both resource consumption ($R$) and
accuracy ($A$) across all models for each instance x, but does not have any
knowledge beyond what is provided within the defined adaptation space.

Conditional assumptions (applied only to certain situations):

  1. 1.

For Equation (8) and the lower bound Equation (9), assume
$\alpha_{i}=P(Y_{1}\neq Y_{GT}\cap Y_{2}\neq Y_{GT}\cap\cdots Y_{i-1}\neq
Y_{GT}|Y_{i}\neq Y_{GT})$ is constant across all $i$’s.

  2. 2.

The constant $\alpha_{i}$ assumption applies to all results in section 3 and
section 4 as well.

###  A.2 Proofs

  1. 1.

Proof of the general formula (5):

###### Proof.

Intuitively, we can estimate the value of $R_{oracle}$ using its expected
value, which can be written as:

$\displaystyle R_{oracle}$
$\displaystyle=\sum_{i=1}^{N}R_{i}P(x_{i})+R_{1}P(x_{f}),$ $\displaystyle
A_{oracle}$ $\displaystyle=1-P(x_{f}),$ (14)

where $x_{1}$, $x_{i}$’s, and $x_{f}$ are defined as the following events:

$\displaystyle x_{1}$ $\displaystyle:=\\{Y_{1}=Y_{GT}\\},$ $\displaystyle
x_{i}$ $\displaystyle:=\\{Y_{1}\neq Y_{GT}\cap Y_{2}\neq Y_{GT}\cap\cdots
Y_{i-1}\neq Y_{GT}\cap Y_{i}=Y_{GT}\\},$ $\displaystyle x_{f}$
$\displaystyle:=\\{Y_{1}\neq Y_{GT}\cap Y_{2}\neq Y_{GT}\cap\cdots\cap
Y_{N}\neq Y_{GT}\\}.$

Then, since $e_{i}=\\{Y_{1}\neq Y_{GT}\cap Y_{2}\neq Y_{GT}\cap\cdots
Y_{i-1}\neq Y_{GT}\cap Y_{i}\neq Y_{GT}\\}$, we see that for
$i=2,3,\cdots,N-1$:

$P(x_{i})=P(e_{i-1})-P(e_{i}).$

Moreover, it’s easy to see that $P(x_{1})=1-P(e_{1})$, and
$P(x_{f})=P(e_{N})$. We can then reformulate the expected value formula in
terms of $P(e_{i})$’s to get the general formula:

$\displaystyle R_{oracle}$
$\displaystyle=\sum_{i=2}^{N}R_{i}[P(e_{i-1})-P(e_{i})]+R_{1}(1-P(e_{1})+P(e_{N})),$
$\displaystyle A_{oracle}$ $\displaystyle=1-P(e_{N}).$ (15)

∎

  2. 2.

Calculations for rewriting the general formula in terms of $\alpha$ (proof of
Equation (7)):

###### Proof.

Given the general formula:

$\displaystyle R_{oracle}$
$\displaystyle=\sum_{i=2}^{N}R_{i}[P(e_{i-1})-P(e_{i})]+R_{1}(1-P(e_{1})+P(e_{N})),$
$\displaystyle A_{oracle}$ $\displaystyle=1-P(e_{N}),$ (16)

and Equation (6), we can do the following calculations:

$\displaystyle R_{oracle}$
$\displaystyle=\sum_{i=2}^{N}R_{i}[P(e_{i-1})-P(e_{i})]+R_{1}(1-P(e_{1})+P(e_{N}))$
$\displaystyle=\left(\sum_{i=2}^{N}R_{i}P(e_{i-1})-R_{i}P(e_{i})\right)+R_{1}(1-P(e_{1})+P(e_{N}))$
$\displaystyle=R_{2}P(e_{1})+\sum_{i=3}^{N}(R_{i}-R_{i-1})P(e_{i-1})-R_{N}P(e_{N})+R_{1}(1-P(e1)+P(e_{N}))$
$\displaystyle=R_{2}(1-A_{1})+$
$\displaystyle\sum_{i=3}^{N}(R_{i}-R_{i-1})\alpha_{i}(1-A_{i})-R_{N}\alpha_{N}(1-A_{N})+R_{1}(1-(1-A_{1})+\alpha_{N}(1-A_{N}))$
$\displaystyle=R_{1}-R_{1}(1-A_{1})+R_{1}\alpha_{N}(1-A_{N}))+R_{2}(1-A_{1})-R_{N}\alpha_{N}(1-A_{N})+$
$\displaystyle\sum_{i=3}^{N}(R_{i}-R_{i-1})\alpha_{i}(1-A_{i})$
$\displaystyle=R_{1}+(R_{2}-R_{1})(1-A_{1})+\sum_{i=3}^{N}(R_{i}-R_{i-1})\alpha_{i}(1-A_{i})-\alpha_{N}(R_{N}-R_{1})(1-A_{N}),$

which is what we have in Equation (7).

∎

  3. 3.

Proof of the criteria for choosing $R_{1}$ (Equation (11)):

###### Proof.

Assume a regular state space $\\{S_{i}\\}$ for a set of backbone classifiers.
The adaptive Oracle’s performance is estimated to be:

$\displaystyle{R_{oracle}}$
$\displaystyle=R_{1}+\sum_{i=2}^{N}(R_{i}-R_{1})(A_{i}-A_{i-1}),$
$\displaystyle A_{oracle}$ $\displaystyle=A_{N}$

according to (9) with the $\alpha_{i}=1$ assumption. Then, assume a special
state space $\\{S^{\prime}_{i}\\}$ where all states are identical as in
$\\{S_{i}\\}$, except that the first state $S_{1}$ is replaced by
$S^{\prime}_{1}$, which is a random agent with:

$\displaystyle R^{\prime}_{1}$ $\displaystyle=0,$ $\displaystyle
A^{\prime}_{1}$ $\displaystyle=\frac{1}{C},$

where $C$ is the number of classes in the classification task. Then, the
Oracle’s performance on this special state space is estimated to be:

$\displaystyle{R^{\prime}_{oracle}}$
$\displaystyle=\sum_{i=3}^{N}(R_{i})(A_{i}-A_{i-1})+R_{2}(A_{2}-\frac{1}{C}),$
$\displaystyle A^{\prime}_{oracle}$ $\displaystyle=A_{N}.$

Note here that $A_{oracle}=A^{\prime}_{oracle}$, so a specific choice of
$R_{1}$ only makes sense if $R_{oracle}<R^{\prime}_{oracle}$, or if
$R_{oracle}-R^{\prime}_{oracle}<0$. This then leads to the following
computation:

$\displaystyle R_{oracle}-R^{\prime}_{oracle}$
$\displaystyle=R_{1}+\sum_{i=2}^{N}(R_{i}-R_{1})(A_{i}-A_{i-1})-\sum_{i=3}^{N}(R_{i})(A_{i}-A_{i-1})-R_{2}(A_{2}-\frac{1}{C})$
$\displaystyle=R_{1}+\sum_{i=3}^{N}(R_{i}-R_{1})(A_{i}-A_{i-1})+(R_{2}-R_{1})(A_{2}-A_{1})-\sum_{i=3}^{N}(R_{i})(A_{i}-A_{i-1})-R_{2}(A_{2}-\frac{1}{C})$
$\displaystyle=R_{1}+\sum_{i=3}^{N}(R_{i})(A_{i}-A_{i-1})$
$\displaystyle-\sum_{i=3}^{N}(R_{1})(A_{i}-A_{i-1})-\sum_{i=3}^{N}(R_{i})(A_{i}-A_{i-1})+(R_{2}-R_{1})(A_{2}-A_{1})-R_{2}(A_{2}-\frac{1}{C})$
$\displaystyle=R_{1}-\sum_{i=3}^{N}(R_{1})(A_{i}-A_{i-1})+R_{2}(A_{2}-A_{1})-R_{2}(A_{2}-\frac{1}{C})-R_{1}(A_{2}-A_{1})$
$\displaystyle=R_{1}-R_{1}(A_{N}-A_{2})-R_{1}(A_{2}-A_{1})+R_{2}(A_{2}-A_{1})-R_{2}(A_{2}-\frac{1}{C})$
$\displaystyle=R_{1}-R_{1}(A_{N}-A_{1})-R_{2}(A_{1}-\frac{1}{C})$
$\displaystyle=R_{1}(1-A_{N}+A_{1})-R_{2}(A_{1}-\frac{1}{C})$

which gives $R_{oracle}-R^{\prime}_{oracle}<0$ if and only if
$R_{1}<\frac{(A_{1}-\frac{1}{C})}{(1-A_{N}+A_{1})}R_{2}$, as required in (11).
∎

###  A.3 Broader Societal Impact

As highlighted in the introduction, enhancing the efficiency of neural
networks is imperative, given the economic and environmental costs associated
with today’s large-scale models. This study establishes the foundation for
addressing these concerns by promoting more efficient model utilization during
the inference phase. However, optimizing system efficiency may impact
performance and robustness if not executed carefully. This motivates further
research into calculating robustness bounds, in addition to the efficiency and
performance bounds introduced in this paper.


# Understanding Representations by Exploring Galaxies in Chemical Space

Jan Weinreich  University of Vienna, Faculty of Physics, Kolingasse 14-16,
AT-1090 Wien, Austria  University of Vienna, Vienna Doctoral School in
Physics, Boltzmanngasse 5, 1090 Vienna, Austria  Konstantin Karandashev
University of Vienna, Faculty of Physics, Kolingasse 14-16, AT-1090 Wien,
Austria  Guido Falk von Rudorff  University Kassel, Department of Chemistry,
Heinrich-Plett-Str.40, 34132 Kassel, Germany  Center for Interdisciplinary
Nanostructure Science and Technology (CINSaT), Heinrich-Plett-Straße 40, 34132
Kassel<EMAIL_ADDRESS>
###### Abstract

We present a Monte Carlo approach for studying chemical feature distributions
of molecules without training a machine learning model or performing
exhaustive enumeration. The algorithm generates molecules with predefined
similarity to a given one for any representation. It serves as a diagnostic
tool to understand which molecules are grouped in feature space and to
identify shortcomings of representations and embeddings from unsupervised
learning. In this work, we first study clusters surrounding chosen molecules
and demonstrate that common representations do not yield a constant density of
molecules in feature space, with possible implications for learning behavior.
Next, we observe a connection between representations and properties: a linear
correlation between the property value of a central molecule and the average
radial slope of that property in chemical space. Molecules with extremal
property values have the largest property derivative values in chemical space,
which provides a route to improve the data efficiency of a representation by
tailoring it towards a given property. Finally, we demonstrate applications
for sampling molecules with specified metric-dependent distributions to
generate molecules biased toward graph spaces of interest.

##  I Introduction

Chemical compound space1 (CCS) encompasses all valid molecular graphs,
occurring naturally, synthesized, or existing only in theory. Investigating
this space is essential for enhancing predictive models of molecular
properties. Molecular representations,2; 3 whether based on graphs or
embedding vectors, qualify CCS as a metric space. While the concept of CCS is
well established there was less attention on studying its topology4 induced by
metric and representation – specifically, the set of molecules in the
neighborhood of any arbitrarily chosen molecule.5 Instead, typically a set of
molecules is predefined through a database and the diversity of molecules can
be visualized using dimensionality reduction techniques.6; 7; 8 This is to not
be confused with the underlying topology of CCS which is _independent_ of the
specific set of selected compounds.

Similarity-based learning methods assume if two representation vectors are
close so are the corresponding quantity values. Since grouping of molecules
depends on the representation vector we argue that understanding the emergent
topology is key to improving machine learning models that employ said
representations. In addition to uniqueness,9 a representation ideally contains
quantities that correlate with properties of interest. Furthermore, adapting
representations to reflect the physics at hand was shown to be beneficial for
model accuracy.10

This raises the question of whether a one fits all solution of representation
design is fundamentally flawed. Indeed, many state-of-the-art neural
networks,11; 12; 13; 14 use embeddings that are learned implicitly for each
property. Except for a few examples,15 this is not the case for physics-based
representations defined by analytical expressions. Supervised molecular
regression methods assume that a _good_ representation maps molecules with
similar properties close to each other and vice versa. Our method explores
chemical graph space obeying valence rules with operations inspired by genetic
algorithms,16; 17; 4; 18 and maintaining detailed balance to guarantee exact
sampling of the distributions of interest. Specifically, the implementation is
based on the Evolutionary Monte Carlo19 algorithm for generating molecular
graphs Molecular Optimization by Sampling Adaptively in Chemical Space20
(MOSAiCS). We use well-established molecular representation vectors2 and
generate molecules based on similarity in the original feature space, not in
projected space. This allows sampling distributions of molecular distances to
investigate the average properties of sets of compounds without the need for
exhaustive enumeration. Systematic control of the similarity between a
molecule of interest and generated molecules is achieved with potential
functions dependent on interpretable distances.

By identifying molecules with high similarity the proposed method can help
improve property-specific representations which remains a challenge _e.g._ for
learning HOMO-LUMO gap values21; 22 or addressing activity cliffs.23 Key
distinctions from generative models are: (i) no projection into latent space
or training is required. (ii) Unlike string-based methods such as EDBO24 or
REINVENT25 derived from SMILES26 or SELFIES,27; 28 the approach is feature
vector agnostic. (iii) Physics-based distance metrics and geometrically
invariant representations are used, offering a more thorough approach than
string-based encoding: generation of new molecules can be based on the
statistical ensemble of conformers by using an average representation.29

Here our focus is the topology generated by a representation and, to a lesser
degree, its correlation to properties. Nonetheless, such correlations of
representations with properties of interest are important for improving
representation efficiency for ML.30 Prevalent methods for probing chemical
space are generative models31 based on unsupervised learning mapping molecules
into a low-dimensional latent space, e.g. SELFIES using autoencoders,32; 33
recurrent neural networks,34; 35; 36 generative adversarial networks,37; 38
guided diffusion models,39 or stochastic processes on molecular graph space.40
While deep learning methods will likely continue leading molecule
generation,41 they may lack chemical interpretability or be less accurate when
applied far from the training set.42 We emphasize that no alternative to
generative models is provided, but an exact sampling of distributions of
chemical space in terms of distances. At best, this contribution aims at
bridging the gap between purely data-driven unsupervised generative models and
physics-inspired representations by combining exact feature vectors and Monte
Carlo (MC) on chemical graphs. However, conceptually the approach may have the
following advantages: (i) No reconstruction problem,42 _i.e._ finding a valid
molecular graph that corresponds to the optimized point in latent space; (ii)
introducing constraints, such as the number of heavy atoms or types of allowed
bonds is trivial; (iii) representation vectors with interpretable features can
be used.

First, we introduce potential functions defined in chemical space and present
details of the algorithm (Sec. II.1). Secondly, we discuss results based on
the sampling of local chemical environments of molecules (Sec. III). Finally,
we analyze the relationship between property values and representation
distances (Sec. III.4).

##  II Methods

###  II.1 Potentials in Chemical Space

We study the topology of chemical space, by examining sets of molecules of
varying “closeness” to each other by generating distributions of molecular
graphs satisfying distance conditions. The method applies to any metric,
including rather sophisticated ones.43; 44 Here we select the distance between
two molecular graphs $A$ and $B$ as

$\displaystyle d(A,B)=||\mathbf{X}(A)-\mathbf{X}(B)||_{2},$ (1)

where $\mathbf{X}$ is a representation vector and $||\ldots||_{2}$ is the
Euclidean norm. Smaller values of $d(A,B)$ mean greater similarity between $A$
and $B$. Although kernel functions also capture non-linear patterns, using
distances rather than kernel elements obviates the need for selecting kernel-
width parameters and enables immediate use of any accessible representation.
We will focus on sets of molecules drawn from a distribution that is uniform
over all such molecules $B$ that $\gamma\leq d(A_{\mathrm{c}},B)\leq\sigma$,
where $A_{\mathrm{c}}$ is a _central molecule_. This is achieved by a
potential function defined in terms of distance to the central molecule (from
now on $d$ for short)

$\displaystyle V(d)=\begin{cases}V_{0}(d-\gamma)^{2}&\text{if}\ d<\gamma\\\
0&\text{if}\ \gamma\leq d\leq\sigma\\\ V_{0}(d-\sigma)^{2}&\text{if}\
d>\sigma~{}.\end{cases}$ (2)

Similarly to parallel tempering,45; 46; 47 each replica is assigned an inverse
temperature parameter $\beta_{i}$ generating a trajectory whose probability
density converges to the Boltzmann distribution with unnormalized distribution
function $\rho_{i}$

$\rho_{i}(d)=\exp[-\beta_{i}V(d)].$ (3)

The potential $V(d)$ has a flat plateau depicted in Fig. 1, with the resulting
$\rho_{i}(d)$ distributions being biased towards $\gamma\leq d\leq\sigma$.
$\beta_{i}$ include infinitely large values (see _greedy replicas_ in Ref. 20)
that can only accept trial moves to graphs with potential functions smaller or
equal to the current one. Therefore, as the simulation progresses for each
greedy replica the potential decreases until the replica reaches the plateau.
Once it is inside the plateau it can never leave it by definition. The
presence of finite values of $\beta_{i}$ (see _exploration replicas_ in Ref.
20) ensures our MC simulations do not get obstructed in moving between subsets
of chemical space separated by barriers of $V$. The $\beta_{i}$ and curvature
values $V_{0}$ of the potential are listed in the Supplementary Information
(SI).

Figure 1: Two quadratic potential functions are defined by Eq. (2). Blue
potential with $\gamma=0$ and $\sigma=2$ defined to sample molecules up to a
distance of $\sigma=2$. Orange curve with $\gamma=2$ and $\sigma=4$ is adapted
to probe molecules in the distance interval $2\leq d\leq 4$.

Here we study the following choices for representations:

  1. 1.

ECFP448 fingerprint, a group-based representation with bit resolution
$n_{\mathrm{bits}}$, unless specified otherwise, $n_{\mathrm{bits}}=2048$. The
diameter parameter, here four, is the maximum distance of neighborhoods of
each atom.

  2. 2.

Smooth overlap of atomic positions49 (SOAP) with averaging over atomic
sites,50 for the remainder of the text referred to as SOAP. 3D conformations
are generated with the Morfeus51 package using MMFF94 forcefield,52;
*Halgren:1996_II; *Halgren:1996_III; *Halgren_Nachbar:1996_IV;
*Halgren:1996_V; *Halgren:1999_VI; *Halgren:1999_VII; 51 as implemented59 in
RdKit,60 and the recipe from Ref. 61. Any method that can generate conformers
from graphs62 can be used. Conformer-averaging of the representation accounts
for the variance of the distance within the classical Boltzmann ensemble for a
given temperature.29

We focus on physics-based representations,3 but representations from
unsupervised learning, such as SchNet,11 may also be used. It is important to
identify relevant distances within the chemical space of interest, as the same
metric evaluated on different representations can differ significantly. For
example, one SOAP distance unit may correspond to a conformation but could
represent a change of constitution for ECFP4. Since the potential is based on
an absolute measure of similarity and predefined representation functions, the
range of representations must be tailored to the specific chemical space of
interest or machine learning application. For instance, comparing molecules of
different sizes may yield larger distances in representation vectors than
comparing molecules of the same size. Another example is how minor changes of
the molecular structure corresponding to small differences in the
representation can result in considerable differences in the property of
interest, e.g. for the HOMO-LUMO gap, necessitating application-dependent
adjustments.22

###  II.2 ChemSpaceSampler code

The Python code for the ChemSpaceSampler is available at
https://github.com/chemspacelab/mosaics. Notebooks with examples can be found
here
https://github.com/chemspacelab/mosaics/tree/main/examples/05_chemspacesampler.
As illustrated in Fig. 2 usage steps are as follows: (1) define the central
molecule, (2) representation, (3) potential type, and (4) constraints for the
chemical space, such as permitted bonds, number of heavy atoms, along with
simulations details such as number of MC steps, moves, and set of $\beta_{i}$.
In practice, using the code is as simple as defining a set of parameters and
calling the main function as follows:

⬇

1params = {

2 ’min_d’: 2.0,

3 ’max_d’: 5.0,

4 ’Nsteps’: 500,

5 ’possible_elements’: ["C", "O", "N", "F"],

6 ’forbidden_bonds’: [(8, 9), (8, 8), (9, 9), (7, 7)],

7 ’nhatoms_range’: [1,9],

8 ’betas’: [8, 4, 2],

9 "nBits": 2048,

10 ’rep_name’: ’ECFP’

11 }

12

13MOLS, D =
chemspace_potentials.chemspacesampler_ECFP(smiles="CC(=O)OC1=CC=CC=C1C(=O)O",
params=params)

The parameter dictionary comprises setting the desired minimal and maximal
distance of the constant potential, the number of MC steps, and possible
elements. It includes the option of excluding certain bond combinations. The
chemical spaces used in this work were inspired by the QM9 dataset63 which
consists of 134k molecules containing up to 9 heavy atoms (C, O, N, and F).
Unless specified otherwise, we will be considering the ”QM9sub” chemical space
which is the set of molecules that contains up to 9 heavy atoms of the same
elements but with the bond constraints forbidding covalent bonds O-F, O-O,
F-F, and N-N, as the latter are associated with increased chemical reactivity.
In addition, the number of bits $n_{\mathrm{bits}}=$ nBits, inverse
temperature, and representation type are defined. Finally, the simulation is
initialized with a SMILES string returning the encountered molecules and their
distances from the central molecule.

Figure 2: ChemSpaceSampler: (1) central molecule, (2) 2D or 3D
representations, and (3) potential function, s. Eq. (2), are defined. The (4)
MC simulation is initialized for a given number of steps, temperatures for the
replica exchanges as well as permitted elements and bonds, and number of heavy
atoms.

##  III Results

###  III.1 Diversity of Chemical Space

We sample CCS of five randomly selected central molecules (s. Fig. 3a) of
which four are from the QM9sub space with ECFP4 at $n_{\mathrm{bits}}=2048$
for $\gamma=0$, $\sigma=6$, see Eq. (2). In addition, we study aspirin, from
QM13sub with 13 heavy atoms. A comparison between the QM963 dataset (grey) and
the ChemSpaceSampler QM9sub results is shown as a two-dimensional UMAP64
visualization in Fig. 3b. A diverse chemical space is explored spanning
regions not covered by QM9.63 Sampling the chemical neighborhood of only five
distinct molecules is sufficient to span a chemical space much larger than the
original QM9 dataset – even though the sampled QM9sub space is more restricted
due to additional bond constraints. This demonstrates that even the use of
central potentials imposing locality results allows generation of chemically
diverse sets of molecules, creating molecules that distinctly fall outside of
the training set.

Figure 3: UMAP projection of the QM9sub chemical space explored by the
ChemSpaceSampler. Sampling the local chemical space of five central molecules
results in a diverse range of molecular structures, as illustrated by the
dispersion of points throughout the plot: For instance, purple dots correspond
to molecules found by ChemSpaceSampler using benzene as a central molecule.
Grey crosses represent molecules from the QM9 dataset.

###  III.2 Galaxies

We examine the chemical spaces surrounding the molecules from the previous
section and perform a principal component analysis of the ECFP4 fingerprints
for $32\leq n_{\mathrm{bits}}\leq 8192$. The number of heavy atoms is fixed to
exactly six for benzene and thirteen for aspirin. Histograms for the density
of the first two principal components of representation vectors for benzene,
aspirin, and nonane at $n_{\mathrm{bits}}=512$ and $8192$ are shown in Fig. 4.

Figure 4:  Chemical spaces with a fixed number of heavy atoms: benzene six,
nonane nine, and aspirin thirteen. We show density plots of the first two
principal components (PCA1 and PCA2) obtained from the sampled chemical space.
The color gradient denotes the distance $\tilde{d}$ (normalized for maximum
distance for each $n_{\mathrm{bits}}$) to the initial molecule in each
fingerprint (FP) dimension. Central molecules shown in the bottom panes are
represented as white stars in the PCA plots. The upper panel for bit-length
$n_{\mathrm{bits}}=512$, lower panel for $n_{\mathrm{bits}}=8192$. Crosses
within the visualization indicate the positions of molecules closest to the
respective cluster centers, all listed in the SI.

After applying a $k$-nearest neighbor algorithm the optimal number of clusters
is estimated using the elbow method by studying the explained variation as a
function of the number of clusters. Afterward, the most representative
molecule for each sub-cluster (s. SI) was identified as the compound with the
smallest distance to the $k$-means centers, cf. Fig. 4. For
$n_{\mathrm{bits}}=8192$, the chemical space around benzene and nonane shows
distinct clusters. For benzene, these clusters are even completely
disconnected suggesting that the convex hull of all clusters contains regions
with vastly different densities of valid molecular graphs. No such seemingly
disconnected clusters were found for the larger aspirin molecule indicating
that the clustering is partially due to the upper bound on the number of heavy
atoms.

Next, we study the capability of ECFP4 to distinguish chemical graphs as a
function of $n_{\mathrm{bits}}$. Decreasing the resolution from
$n_{\mathrm{bits}}=8196$ to $n_{\mathrm{bits}}=512$ results in smaller
distances between cluster centers and the merging of previously distinct
clusters; thus $n_{\mathrm{bits}}$ affects the separation of different
chemical environments. The number of clusters depends on an arbitrary
threshold parameter for the elbow method and also on the chemical space under
consideration, i.e. diverse chemical spaces may result in a larger number of
clusters. Thus a more robust measure for the resolution of CCS should be
defined. To remove any dependence on the number of clusters, we evaluate as a
function of fingerprint length the average distance between molecules closest
to the center of each cluster

$\displaystyle\bar{d}(n_{\mathrm{bits}})$
$\displaystyle=\frac{1}{M}\sum_{l=1}^{M}\frac{1}{C^{l}}\sum_{{k\neq
s}}^{C^{l}}d(A^{l}_{k},A^{l}_{s})$ (4)

where $M$ is the number of considered central molecules, $l=1,...,M$ is the
index of a central molecule. $C^{l}$ is the number of clusters of molecules
observed after running an MC simulation, $k$ and $s$ are indices of such
clusters, and $A_{k}^{l}$ is the molecule closest to the center of such a
cluster with index $k$. Each two-dimensional PCA plot contains up to $C^{l}$
clusters. Finally, the distances are averaged over all $M$ considered
molecules to ensure $\bar{d}$ does not depend on the individual chemical
diversity around each of the central molecules. Intuitively,
$\bar{d}(n_{\mathrm{bits}})\rightarrow 0$ corresponds to the merging of
previously separated clusters whereas a higher resolution leads to larger
average distances between clusters. At $n_{\mathrm{bits}}\approx 1024$, the
average distance appears to saturate (cf. Fig. 5b) supporting the choice of
the default value $n_{\mathrm{bits}}=1024$ appearing in RdKit. The closest and
farthest molecules for benzene and aspirin are shown in Fig. 5a. Closest point
sampling for aspirin results in small modifications of the ester group while
the rest of the molecule is unaffected. For benzene, the closest points are
hydrocarbons obtained by adding extra bonds to the central molecule, whereas
at larger distances linear molecules and different compositions are found.
Note that ECFP4 may result in similar or even completely identical distance
values even if the molecular graphs differ substantially. Thus sampling at a
given distance can help identify molecules at comparable distances yet
possessing vastly distinct chemistries. Such cases are relevant if properties
of interest change strongly between any two similar considered molecules and
typically result in particularly poor ML predictions.30

Figure 5:  (a): closest and farthest molecules encountered for benzene and
aspirin respectively using ECFP4 at $n_{\mathrm{bits}}=8192$. (b): mean
distance $\bar{d}$ of the clusters averaged over five different molecules, cf.
Eq. (4).

###  III.3 Molecules with predefined similarity

The core assumption of any similarity-based model is that systems of similar
representations are of similar properties $P$. We exploit this observation to
generate molecules with similar representation vectors and properties. First,
we select 52 QM9sub molecules as central molecules and perform $N=5000$ MC
steps for ECFP4 and $N=2000$ for SOAP (smaller number of MC steps caused by
the relatively high cost of generating conformers), permitting up to nine
heavy atoms. Finally, we evaluate the free energy of solvation $P=G$ and band
gap $P=\varepsilon$ of the sampled compounds based on GFN2-xTB65 calculations
and an analytical linearized Poisson-Boltzmann model66 to simulate the
presence of water using the same procedure as described in Ref. 20. Next, we
perform a principal component analysis of the initial molecules and their
sampled neighbors, as depicted in Fig 6. As expected from the superior
performance of SOAP and geometry-based molecular representations in most
regression tasks, better clustering of sampled molecules around the initial
molecules compared to the ECFP4 fingerprint is observed. Conversely, ECFP4
does not require chemical graphs to be geometrically constrained. The number
of encountered molecules surrounding the central molecule is very susceptible
to its chemical graph. This supports the observation of regions with low and
high molecular densities and indicates distinct topologies depending on the
location in chemical space. For example, with ECFP4 and ethanol, we identify
574 unique molecules in $2.75\leq d_{\text{ECFP4}}\leq 3.25$, but no other
molecule is found for NC1=NC(=O)N=C(N)N1 in the same interval (cf.
$\textbf{C}_{6}$ in Fig. 7). By contrast, for SOAP and the interval $50\leq
d_{\text{SOAP}}\leq 60$, no other molecule is found for ethanol, but for
NC1=NC(=O)N=C(N)N1, we find 21 molecules suggesting that the density of
molecules around central molecules is highly dependent on the representation
type. A more quantitative comparison of the density is challenging, as there
seems to be no straightforward way to compare distances between
representations. Systematic generation of molecules far from an initial
molecule may find application for farthest point sampling often used for
active learning strategies.67; 68 Conversely, generating molecules close to a
query molecule may help identify the most relevant neighbors in chemical space
for query-specific ML predictions.69; 70

Figure 6:  ChemSpaceSampler with ECFP4 (top) SOAP49 (bottom) representations:
comparing sampled chemical space projected principal components (PC) into two
dimensions along axis PC1 and PC2 with maximal variance. Large marker symbols
represent central molecules used in the potential function, Eq. (2). Small
opaque markers with matching colors and styles denote molecules obtained after
sampling. (a) distance intervals at small representation distances: $2.75\leq
d_{\text{ECFP4}}\leq 3.25$ for ECFP4 and $50\leq d_{\text{SOAP}}\leq 60$ for
SOAP; (b) same at larger distances: $4.00\leq d_{\text{ECFP4}}\leq 4.25$ and
$100\leq d_{\text{SOAP}}\leq 110$. Figure 7: Identified central molecules
$\mathbf{C}_{l}$ with extremal radial average property slope values cf. Tab.
1.

###  III.4 Property derivatives in chemical space

In the following we use atomic units. Given a central molecule $l$ one could
consider average free energy of solvation $\bar{G}_{l}(d)$ and band gap
$\bar{\varepsilon}_{l}(d)$ of molecules encountered at a distance $d$, _i.e._
average property values at $d$. However, an average over an exact distance can
be ill-defined due to the discreteness of chemical graph space and would be
impractical to evaluate with MC methods even if it was. Instead, we consider
averages over shells because this allows numerical evaluation. Therefore we
introduce binning and finite distance intervals with outer bin positions given
by

$\displaystyle d_{\text{ECFP4}}^{b}=\\{1.25+b\cdot 0.5\\}~{},$ (5)

where $b\in 0...,10$. This results in discrete _shell_ intervals of

$\displaystyle
I_{\text{ECFP4}}=\\{[0.75,1.25],\underbrace{(1.25,1.75]}_{b=1},...,(5.75,6.25]\\}$
(6)

for ECFP4 and SOAP

$\displaystyle d_{\text{SOAP}}^{b}=\\{50+b\cdot 10\\},$ (7)

where $b\in 0,...,10$. The shell intervals $I_{\text{SOAP}}$ are defined
analogously based on preliminary simulations. The bin centers are chosen
according to meaningful distances for each representation respectively. For
instance, at a SOAP distance of $d=50$, only a single bond may differ from the
central molecule, whereas at $d=150$ encountered molecules will hardly
resemble the central molecule.

Note that the distances do not start at zero as we could not find any
molecules for smaller distances. Subsequently and for all 52 central molecules
$l$ we will consider, the average of both properties within each finite shell
interval is computed and denoted as $\bar{P}(d)$ (the property $P$ can be $G$
or $\varepsilon$), _e.g._ $\bar{G}(d=1.25)$ is the average value of free
energy of solvation over all molecules in the interval $[0.75,1.25]$ for
ECFP4. More generally, we define shell averages as

$\displaystyle\bar{P}_{l}(d^{b})=\frac{1}{M_{b}}\sum_{j=1}^{M_{b}}P_{k},$ (8)

where $M_{b}$ is the number of molecules in shell $b$. These are molecules
whose distance to the central molecule $l$ is between $d^{b-1}$ and $d^{b}$.
The indices $j\in 1,...,M_{b}$ enumerate these molecules. The property value
of the central molecule $l$ is $\bar{P}_{l}(0)$. The definition of shell
averages is also depicted in Fig. 8. For simplicity the index $b$ representing
the shell index is removed and $d$ always represents the upper shell boundary.

Figure 8:  Illustration of the local chemical environment of ethanol obtained
from simulations. The first shell is defined up to a SOAP distance of $d_{1}$,
and a second shell contains molecules $k$ in the interval $d_{1}\leq
d(0,l)\leq d_{2}$. Evaluating the encountered molecules within different
shells allows extracting shell averages $\bar{P}$, c.f. Eq. (8).

Considering the number of encountered molecules grows quickly with distance we
introduce a cutoff, evaluating only up to 300 randomly selected molecules per
shell to limit computational costs. As shown in Fig. 9a for small distances,
$d_{\text{SOAP}}\leq 60$, property values of new molecules are similar to
those of the two central molecules. Increasing the distance further leads to a
larger drift in the property values – consistent with the assumption of
regression models that small distances $d$ between training and test molecules
correlate to similar $P$.

Figure 9:  (a): radially averaged solvation energy $\bar{G}(d)$ and HOMO-LUMO
gap $\bar{\varepsilon}(d)$ as a function of SOAP distance $d$ for molecules
shown in the bottom panel. (b): Absolute values of radially averaged property
shifts $|\overline{P}(d)-P(0)|$ compared to the values $P(0)$ of several
initial molecules (opaque) for increasing ECFP4 distances $d$. Black dashed
lines show respective $95^{\text{th}}$ percentiles. The blue line with error
bars shows the average shift overall considered molecules at each $d$. top:
$P=G$, solvation energy, bottom $P=\varepsilon$, HOMO-LUMO gap.

Collecting the first two principal components of all encountered molecules
resulting from SOAP simulations shows how the HOMO-LUMO gap
$\bar{\varepsilon}$ (Fig. 10a) changes between molecules. The corresponding
figure for the free energy $\bar{G}$ is shown in SI Fig. 2. Both $\bar{G}(d)$
and $\bar{\varepsilon}(d)$ vary substantially from molecule to molecule, as
shown in Fig. 9b making it challenging to extract general trends for
$\bar{P}(d)$. Still, we fit linear functions to $\bar{P}_{l}(d)$ to
approximate average trends of properties through chemical space for each
molecule $l$ as follows,

$\bar{P}_{l}(d)\approx m^{P}_{l}\cdot d+c^{P}_{l},$ (9)

where $m^{P}_{l}$ and $c^{P}_{l}$ are the fitted slope and intercept. After
extracting $\bar{\varepsilon}_{l}(d)$ and $\bar{G}_{l}(d)$ for each molecule
and each $d$, we compute the slope values $m^{G}_{l}$ and
$m^{\varepsilon}_{l}$ by regressing Eq. (9) to the outer shell values as $d$
and the 10 shell averages $\bar{P}_{l}(d)$. Note that $c^{P}_{l}$ is
determined by the value of the central molecule. Scatter plots displaying the
slope values $m^{P}_{l}$, normalized by the average pairwise distance of all
central molecules $\langle d_{ks}\rangle$, against the property values of the
central molecules ($d=0$) are presented in Fig. 11a,b.

A near linear correlation between the radially averaged property $m_{l}^{P}$
derivative and the central property value is found,
$m_{l}^{P}\propto\bar{P}_{l}(0)$. For instance, the extracted slope values for
the solvation free energies $m_{l}^{G}$ correlate linearly with the initial
value of the respective central molecule $G_{l}(0)$. Average property values
increase (decrease) for molecules from the tails of the property distribution
corresponding to smaller (larger) values for increasing shell distances.
Radial average property values of sampled molecules from initial compounds
with more negative solvation energies than most other molecules drift towards
less soluble compounds and vice versa. The same trend applies to the HOMO-LUMO
gap $\varepsilon$. We also find that central molecules with average property
values have a near zero slope, indicated by vertical and horizontal lines in
Fig. 11ab respectively.

Figure 10:  (a): first two principal components (PC) of central molecules,
shown with markers, as well as sampled molecules, visited in all shells
$d\in\bigl{[}50,150\bigr{]}$ represented as dots colored by individual HOMO-
LUMO gap values $\varepsilon$. (b): scatter plot between solvation energy $G$
and band gap $\varepsilon$ for all molecules obtained by sampling. In both
cases SOAP was used. Figure 11:  Slopes $m^{P}_{l}/\langle d_{ks}\rangle$ of
average changes of a property [(a) free energy of solvation $G$ or (b) band
gap $\epsilon$] with respect to $d$ for molecules surrounding a central
molecule using ECFP4 (SOAP shown in the SI Fig. 1, and values given in Tab.
1). The slope values are scaled by the average pairwise distance $\langle
d_{ks}\rangle$ of all 52 molecules. Each dot corresponds to a molecule $l$.
y-axis displays the slope of a linear function fitted to Eq. (9) for each
property. x-axis represents the central molecule’s property. Molecules with
the largest and smallest slopes are shown in the inset. Top: solvation energy
$G$, bottom HOMO-LUMO gap. The vertical black dashed line indicates the
property average over all 52 central molecules and the horizontal line
indicates zero slope.

To investigate whether finding the linear trend for both band gap and
solvation energy is partially due to a correlation between both quantities, we
also show a scatter plot of both in Fig. 10b revealing that the linear slope
holds regardless since the properties themselves do not correlate. Given
$m^{P}_{l}\propto\bar{P}_{l}(0)$, an upper bound of the average property slope
$\bar{P}^{\prime}=m^{P}$ in a chemical space can be estimated by the compound
with the largest property value,

$\displaystyle\max\bar{P}^{\prime}(d)=\max
m^{P}=\max\limits_{l\in\text{CCS}}\left[\tilde{m}^{P}\bar{P}_{l}(0)+\tilde{c}^{P}\right],$
(10)

where $\tilde{m}^{P}$ and $\tilde{c}^{P}$ are linear slope and intercept
extracted from the correlation plots in Fig. 11ab. The most negative slope on
the other hand is defined by the molecule with the minimal property value.

The regression constants $\tilde{m}^{P}$, $\tilde{c}^{P}$ as well as the
average, minimal, and maximal slope values and Spearman’s rank correlation
coefficient $r_{s}$, measuring monotonic correlation between
$\bar{P}^{\prime}$ and $\Delta_{l}^{P}$, for both properties and
representations are listed in Tab. 1. The identified extremal molecules are
shown in Fig. 7. As expected $\tilde{m}^{P}$ and $\tilde{c}^{P}$ depend on the
property, representation, and permitted CCS. For both representations, the
magnitude of the maximal slope for $\bar{G}$ corresponds to the same compound
$\mathbf{C}_{6}$ (s. Fig. 7) which is also the most soluble in QM9. The
correlation coefficient $r_{s}$ for the SOAP representation is smaller than
that of the ECFP4 fingerprint. This is evident from the values recorded for
$G$, which are $-0.788$ and $-0.976$ respectively, as shown in Tab. 1. We
propose two potential explanations for the higher correlation coefficient of
ECFP4. First, we managed to perform more MC steps for ECFP4 simulations, thus
potentially making them more converged. Second, the conformer sampling as
required for the SOAP workflow introduces statistical noise due to a degree of
randomness in generating conformers used in SOAP representation.

Table 1: Minimal and maximal slope $m^{P}$ of the average radial property derivative in chemical space and corresponding molecules $\mathbf{C}_{i}$ identified among 52 QM9sub molecules as well as the Spearman correlation coefficient, shown in Fig. 7. Results for both molecular representations, ECFP448 as well as SOAP.49 | $P$ | Min mol | Max mol | $\text{min}~{}m^{P}$ | $\text{max}~{}m^{P}$ | $\tilde{m}^{P}$ | $\tilde{c}^{P}$ | $r_{s}$
---|---|---|---|---|---|---|---|---
ECFP4 | $G$ | $\textbf{C}_{1}$ | $\textbf{C}_{6}$ | $-2.02\times 10^{-3}$ | $6.42\times 10^{-3}$ | $-2.24\times 10^{-1}$ | $1.94\times 10^{-3}$ | -0.976
$\varepsilon$ | $\textbf{C}_{1}$ | $\textbf{C}_{5}$ | $-8.04\times 10^{-2}$ | $1.58\times 10^{-2}$ | $2.02\times 10^{-1}$ | $1.94\times 10^{-2}$ | -0.988
SOAP | $G$ | $\textbf{C}_{2}$ | $\textbf{C}_{6}$ | $-7.32\times 10^{-5}$ | $8.74\times 10^{-5}$ | $-5.24\times 10^{-3}$ | $4.66\times 10^{-5}$ | -0.788
$\varepsilon$ | $\textbf{C}_{3}$ | $\textbf{C}_{4}$ | $-3.41\times 10^{-3}$ | $4.99\times 10^{-4}$ | $-5.46\times 10^{-3}$ | $-7.47\times 10^{-4}$ | -0.874

Since Eq. (10) describes a relation over all of the permitted chemical space,
even approximate values for the maximal $m^{P}$ may help design molecular
representations or tuning model flexibility by taking into account distances
between data points and the expected variation in the target property. As an
illustration, we show the absolute property deviation from the initial values
in Fig. 9a,b: A too small slope of the curve under the $95^{\text{th}}$
percentile means weak correlation with the target property. Too large slopes
indicate discontinuities of the representation with respect to the property,
which are associated with larger errors.22 A good representation should find a
compromise between both. By calculating the derivative of the property machine
learning models could be fine-tuned for each property, enhancing a model’s
performance.

##  IV Conclusion

An MC-based method was introduced to identify shortcomings of molecular
representations designed for supervised learning tasks. It enables exact
sampling of distributions of molecules defined in terms of distances within
the original feature space with greater explainability than deep learning
approaches. Access to molecules with given similarity may help design new
feature-based representation vectors by tailoring the average property change
with respect to representation distance (s. Fig. 9a,b). The sampling
probabilities can be accessed by a simple comparison of molecular features.
This could be particularly useful for interpretable descriptors commonly used
in cheminformatics.

A near linear trend between the initial value of the central molecule and the
radial slope, $m^{P}_{l}\propto\bar{P}_{l}(0)$, was found and allows
estimating a lower bound for the maximal slope given by the molecule with the
largest property value [s. Eq. (10)]. We hypothesize this trend is due to the
difference between the initial molecule and the average property value, but
this observation warrants further investigation. Several aspects could be
investigated building on this work: (i) incorporating control for synthetic
accessibility71 into the generation of molecules will further enhance its
utility. (ii) Identifying pairs of molecules that are close in representation
but have significantly different properties, sometimes referred to as activity
cliffs,30 presents a challenge for existing predictive models – in particular
for ligand-based drug-design applications of ML. (iii) Investigating the
effect of different choices of metrics on the structure of CCS, for instance,
Tanimoto72 or Wasserstein.73 Finally, the method presented here may also serve
as a tool to invert molecular representations with applications to molecular
property optimization.

##  V Code Availability

The Python code and graphical user interface for the ChemSpaceSampler is
available at https://github.com/chemspacelab/mosaics. Notebooks with examples
can be found here
https://github.com/chemspacelab/mosaics/tree/main/examples/05_chemspacesampler.

##  VI Acknowledgments

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No 957189. Obtaining
the presented computational results has been facilitated using the queueing
system implemented at http://leruli.com. J.W. acknowledges support from the
Faculty of Physics and supervision by Anatole v. Lilienfeld and C. Dellago at
the University of Vienna. The authors thank Michael Sahre and Alexandre
Schoepfer for their helpful discussions and feedback.

##  VII Supplementary Information

The following Supplementary Information contains a list of the selected
inverse temperature values $\beta_{i}$, the SOAP representation figure
corresponding to Fig. 11 as well as the PCA corresponding to the solvation
free energies not shown in Fig. 10a.

## References

  * Kirkpatrick and Ellis (2004) P. Kirkpatrick and C. Ellis, Nature 432, 823 (2004).
  * Musil  _et al._ (2021) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Chemical Reviews 121, 9759 (2021).
  * Parsaeifard  _et al._ (2021) B. Parsaeifard, D. S. De, A. S. Christensen, F. A. Faber, E. Kocer, S. De, J. Behler, O. A. von Lilienfeld, and S. Goedecker, Machine Learning: Science and Technology 2, 015018 (2021).
  * Virshup  _et al._ (2013) A. M. Virshup, J. Contreras-García, P. Wipf, W. Yang, and D. N. Beratan, Journal of the American Chemical Society 135, 7296 (2013).
  * Fu  _et al._ (2021) T. Fu, C. Xiao, X. Li, L. M. Glass, and J. Sun, Proceedings of the AAAI Conference on Artificial Intelligence 35, 125 (2021).
  * Capecchi  _et al._ (2020) A. Capecchi, D. Probst, and J.-L. Reymond, Journal of Cheminformatics 12 (2020), 10.1186/s13321-020-00445-4.
  * Awale  _et al._ (2017) M. Awale, R. Visini, D. Probst, J. Arús-Pous, and J.-L. Reymond, CHIMIA 71, 661 (2017).
  * Willett (2006) P. Willett, Drug Discovery Today 11, 1046 (2006).
  * Huang and von Lilienfeld (2016) B. Huang and O. A. von Lilienfeld, The Journal of Chemical Physics 145, 161102 (2016), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4964627/13765901/161102_1_online.pdf .
  * Huang and von Lilienfeld (2021) B. Huang and O. A. von Lilienfeld, Chemical Reviews 121, 10001 (2021).
  * Schütt  _et al._ (2017) K. T. Schütt, P.-J. Kindermans, H. E. Sauceda, S. Chmiela, A. Tkatchenko, and K.-R. Müller, “Schnet: A continuous-filter convolutional neural network for modeling quantum interactions,” (2017), arXiv:1706.08566 [stat.ML] .
  * Westermayr and Marquetand (2020) J. Westermayr and P. Marquetand, Chemical Reviews ASAP (2020), 10.1021/acs.chemrev.0c00749.
  * Musaelian  _et al._ (2023) A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Nature Communications 14, 579 (2023).
  * von Lilienfeld  _et al._ (2015) O. A. von Lilienfeld, R. Ramakrishnan, M. Rupp, and A. Knoll, International Journal of Quantum Chemistry 115, 1084 (2015).
  * Helfrecht  _et al._ (2020) B. A. Helfrecht, R. K. Cersonsky, G. Fraux, and M. Ceriotti, Machine Learning: Science and Technology 1, 045021 (2020).
  * Globus  _et al._ (1999) A. Globus, J. Lawton, and T. Wipke, Nanotechnology 10, 290 (1999).
  * Brown  _et al._ (2004) N. Brown, B. McKay, F. Gilardoni, and J. Gasteiger, J. Chem. Inf. Comput. Sci. 44, 1079–1087 (2004).
  * Jensen (2019) J. H. Jensen, Chem. Sci. 10, 3567–3572 (2019).
  * Liang and Wong (2000) F. Liang and W. H. Wong, Stat. Sin. 10, 317 (2000).
  * Karandashev  _et al._ (2023) K. Karandashev, J. Weinreich, S. Heinen, D. J. A. Arrieta, G. F. von Rudorff, K. Hermansson, and O. A. von Lilienfeld, “Evolutionary monte carlo of qm properties in chemical space: Electrolyte design,” (2023), arXiv:2307.15563 [physics.chem-ph] .
  * Karandashev and von Lilienfeld (2022) K. Karandashev and O. A. von Lilienfeld, The Journal of Chemical Physics 156, 114101 (2022).
  * Mazouin  _et al._ (2022) B. Mazouin, A. A. Schöpfer, and O. A. von Lilienfeld, Mater. Adv. 3, 8306 (2022).
  * van Tilborg  _et al._ (2022a) D. van Tilborg, A. Alenicheva, and F. Grisoni, Journal of Chemical Information and Modeling 62, 5938 (2022a).
  * Shields  _et al._ (2021) B. J. Shields, J. M. Stevens, J. Li, M. Parasram, F. N. Damani, J. I. M. Alvarado, J. M. Janey, R. P. Adams, and A. G. Doyle, Nature 590, 89 (2021).
  * Olivecrona  _et al._ (2017a) M. Olivecrona, T. Blaschke, O. Engkvist, _et al._ , Journal of Cheminformatics 9, 48 (2017a).
  * Weininger (1988) D. Weininger, Journal of Chemical Information and Computer Sciences 28, 31 (1988).
  * Krenn  _et al._ (2020) M. Krenn, F. Häse, A. Nigam, P. Friederich, and A. Aspuru-Guzik, Machine Learning: Science and Technology 1, 045024 (2020).
  * Krenn  _et al._ (2022) M. Krenn, Q. Ai, S. Barthel, N. Carson, A. Frei, N. C. Frey, P. Friederich, T. Gaudin, A. A. Gayle, K. M. Jablonka, R. F. Lameiro, D. Lemm, A. Lo, S. M. Moosavi, J. M. Nápoles-Duarte, A. Nigam, R. Pollice, K. Rajan, U. Schatzschneider, P. Schwaller, M. Skreta, B. Smit, F. Strieth-Kalthoff, C. Sun, G. Tom, G. Falk von Rudorff, A. Wang, A. D. White, A. Young, R. Yu, and A. Aspuru-Guzik, Patterns 3, 100588 (2022).
  * Weinreich  _et al._ (2021) J. Weinreich, N. J. Browning, and O. A. von Lilienfeld, The Journal of Chemical Physics 154, 134113 (2021).
  * van Tilborg  _et al._ (2022b) D. van Tilborg, A. Alenicheva, and F. Grisoni, Journal of Chemical Information and Modeling 62, 5938 (2022b).
  * Vogt (2023) M. Vogt, Artificial Intelligence in the Life Sciences 3, 100064 (2023).
  * Nigam  _et al._ (2021) A. Nigam, R. Pollice, M. Krenn, G. d. P. Gomes, and A. Aspuru-Guzik, Chem. Sci. 12, 7079 (2021).
  * Schilter  _et al._ (2023) O. Schilter, A. Vaucher, P. Schwaller, and T. Laino, Digital Discovery , (2023).
  * Olivecrona  _et al._ (2017b) M. Olivecrona, T. Blaschke, O. Engkvist, and H. Chen, Journal of Cheminformatics 9 (2017b), 10.1186/s13321-017-0235-x.
  * Segler  _et al._ (2018) M. H. S. Segler, T. Kogej, C. Tyrchan, and M. P. Waller, ACS Central Science 4, 120 (2018).
  * Blaschke  _et al._ (2020) T. Blaschke, J. Arús-Pous, H. Chen, C. Margreitter, C. Tyrchan, O. Engkvist, K. Papadopoulos, and A. Patronov, Journal of Chemical Information and Modeling 60, 5918 (2020).
  * Prykhodko  _et al._ (2019) O. Prykhodko, S. Johansson, P.-C. Kotsias, J. Arús-Pous, E. Bjerrum, O. Engkvist, and H. Chen, Journal of Cheminformatics 11 (2019), 10.1186/s13321-019-0397-9.
  * Gómez-Bombarelli  _et al._ (2018) R. Gómez-Bombarelli, J. N. Wei, D. Duvenaud, J. M. Hernández-Lobato, B. Sánchez-Lengeling, D. Sheberla, J. Aguilera-Iparraguirre, T. D. Hirzel, R. P. Adams, and A. Aspuru-Guzik, ACS Central Science 4, 268 (2018).
  * Weiss  _et al._ (2023) T. Weiss, L. Cosmo, E. M. Yanes, S. Chakraborty, A. M. Bronstein, and R. Gershoni-Poranne, chemRxiv (2023), 10.26434/chemrxiv-2023-z8ltp.
  * Lee  _et al._ (2022) S. Lee, J. Jo, and S. J. Hwang, “Exploring chemical space with score-based out-of-distribution generation,” (2022), arXiv:2206.07632 [q-bio.BM] .
  * Westermayr  _et al._ (2023) J. Westermayr, J. Gilkes, R. Barrett, and R. J. Maurer, Nature Computational Science 3, 139 (2023).
  * Muenkler  _et al._ (2023) H. Muenkler, H. Misztela, M. Pikusa, M. Segler, N. Schneider, and K. Maziarz, “Are vaes bad at reconstructing molecular graphs?” (2023), arXiv:2305.03041 [cs.LG] .
  * Çaylak  _et al._ (2020a) O. Çaylak, O. A. von Lilienfeld, and B. Baumeier, Machine Learning: Science and Technology 1, 03LT01 (2020a).
  * Fabregat  _et al._ (2022) R. Fabregat, P. van Gerwen, M. Haeberle, F. Eisenbrand, and C. Corminboeuf, Machine Learning: Science and Technology 3, 035015 (2022).
  * Hukushima and Nemoto (1996) K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
  * Sambridge (2014) M. Sambridge, Geophys. J. Int. 196, 357–374 (2014).
  * Angelini and Ricci-Tersenghi (2019) M. C. Angelini and F. Ricci-Tersenghi, Phys. Rev. E 100, 013302 (2019).
  * Rogers and Hahn (2010) D. Rogers and M. Hahn, J. Chem. Inf. Mod. 50, 742 (2010).
  * Bartók  _et al._ (2013) A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev. B 87, 184115 (2013).
  * Himanen  _et al._ (2020) L. Himanen, M. O. J. Jäger, E. V. Morooka, F. Federici Canova, Y. S. Ranawat, D. Z. Gao, P. Rinke, and A. S. Foster, Computer Physics Communications 247, 106949 (2020).
  * Jorner (2022) K. Jorner, Zenodo (2022), 10.5281/zenodo.6685219.
  * Halgren (1996a) T. A. Halgren, J. Comput. Chem. 17, 490 (1996a).
  * Halgren (1996b) T. A. Halgren, J. Comput. Chem. 17, 520 (1996b).
  * Halgren (1996c) T. A. Halgren, J. Comput. Chem. 17, 553 (1996c).
  * Halgren and Nachbar (1996) T. A. Halgren and R. B. Nachbar, J. Comput. Chem. 17, 587 (1996).
  * Halgren (1996d) T. A. Halgren, J. Comput. Chem. 17, 616 (1996d).
  * Halgren (1999a) T. A. Halgren, J. Comput. Chem. 20, 720 (1999a).
  * Halgren (1999b) T. A. Halgren, J. Comput. Chem. 20, 730 (1999b).
  * Tosco  _et al._ (2014) P. Tosco, N. Stiefl, and G. Landrum, J. Cheminform. 6, 37 (2014).
  * sof (2023) “Rdkit: Open-source cheminformatics,”  (2023), python package.
  * Ebejer  _et al._ (2012) J.-P. Ebejer, G. M. Morris, and C. M. Deane, J. Chem. Inf. Model. 52, 1146–1158 (2012).
  * Lemm  _et al._ (2021) D. Lemm, G. Rudorff, and A. von Lilienfeld, Nature Communications 12 (2021), 10.1038/s41467-021-24525-7.
  * Ramakrishnan  _et al._ (2014) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. Anatole von Lilienfeld, Sci. Data 1, 140022 (2014).
  * McInnes  _et al._ (2018) L. McInnes, J. Healy, N. Saul, and L. Großberger, Journal of Open Source Software 3, 861 (2018).
  * Bannwarth  _et al._ (2019) C. Bannwarth, S. Ehlert, and S. Grimme, J. Chem. Theory Comput. 15, 1652 (2019).
  * Ehlert  _et al._ (2021) S. Ehlert, M. Stahn, S. Spicher, and S. Grimme, J. Chem. Theory Comput. 17, 4250 (2021).
  * Dral (2020) P. O. Dral, in  _Chemical Physics and Quantum Chemistry_, Advances in Quantum Chemistry, Vol. 81, edited by K. Ruud and E. J. Brändas (Academic Press, 2020) pp. 291–324.
  * Cordova  _et al._ (2022) M. Cordova, E. A. Engel, A. Stefaniuk, F. Paruzzo, A. Hofstetter, M. Ceriotti, and L. Emsley, The Journal of Physical Chemistry C 126, 16710 (2022).
  * Sheridan  _et al._ (2004) R. P. Sheridan, B. P. Feuston, V. N. Maiorov, and S. K. Kearsley, Journal of chemical information and computer sciences 44 6, 1912 (2004).
  * Lemm  _et al._ (2023) D. Lemm, G. F. von Rudorff, and O. A. von Lilienfeld, “Improved decision making with similarity based machine learning,” (2023), arXiv:2205.05633 [physics.chem-ph] .
  * Dolfus  _et al._ (2022) U. Dolfus, H. Briem, and M. Rarey, Journal of Chemical Information and Modeling 62, 3565 (2022).
  * Tanimoto (1958) T. T. Tanimoto,  _An elementary mathematical theory of classification and prediction_ , Tech. Rep. (IBM, 1958).
  * Çaylak  _et al._ (2020b) O. Çaylak, O. A. von Lilienfeld, and B. Baumeier, Machine Learning: Science and Technology 1, 03LT01 (2020b).