# Nonmyopic Multifidelity Active Search

Quan Nguyen  Arghavan Modiri  Roman Garnett

###### Abstract

Active search is a learning paradigm where we seek to identify as many members
of a rare, valuable class as possible given a labeling budget. Previous work
on active search has assumed access to a faithful (and expensive) oracle
reporting experimental results. However, some settings offer access to cheaper
surrogates such as computational simulation that may aid in the search. We
propose a model of _multifidelity_ active search, as well as a novel,
computationally efficient policy for this setting that is motivated by state-
of-the-art classical policies. Our policy is nonmyopic and budget aware,
allowing for a dynamic tradeoff between exploration and exploitation. We
evaluate the performance of our solution on real-world datasets and
demonstrate significantly better performance than natural benchmarks.

Multifidelity, Active Search, Nonmyopic, ICML

inputsinputs returnsreturns

##  1 Introduction

The goal of _active search_ is to identify members of a rare and valuable
class among a large pool of unlabeled points. This is a simple model of many
real-world discovery problems, such as drug discovery and fraud detection.
Active search proceeds by successively querying an oracle that returns a
binary label indicating whether or not a chosen data point exhibits the
desired properties. In many applications, this oracle is expensive, limiting
the number of queries that could be made. The challenge is to design a policy
to sequentially query the oracle in order to discover as many targets as
possible, subject to a given labeling budget.

Active search has been extensively studied under various settings (Garnett et
al., 2012; Jiang et al., 2017, 2018, 2019). Notably, Jiang et al. (2017)
proved a hardness of approximation result, showing that no polynomial-time
policy can approximate the performance of the Bayesian optimal policy within
any constant factor. Thus active search is surprisingly hard. However, the
authors also proposed an efficient approximation to the optimal policy that
delivers impressive empirical performance. The key feature of their policy is
its ability to account for the remaining budget and dynamically trade off
exploration and exploitation.

Most previous work on active search has assumed access to only a single
expensive oracle providing labels. In practice, however, we may have several
methods of probing the search space, including cheap, low-fidelity surrogates.
For example, a computational simulation may serve as a noisy approximation to
an experiment done in a laboratory, and we may reasonably seek to use a
cheaper computational search to help design the expensive experiments. This
motivates the problem of _multifidelity active search,_ where oracles of
different fidelities may be simultaneously accessed to accelerate the search
process. The central question in this task is how to effectively leverage
these oracles in order to maximize the rate of discovery.

We present a model for multifidelity active search and study the problem under
the framework of Bayesian decision theory. We propose a novel policy for this
setting inspired by the state-of-the-art single-fidelity policy mentioned
above. The result is an efficient approximation to the optimal multifidelity
policy that is specifically tailored to take advantage of low-fidelity
oracles. We also create aggressive branch-and-bound pruning strategies to
increase the efficiency of our proposed algorithm, enabling scaling to large
datasets. In a series of experiments, we investigate the performance of our
policy on several real-world datasets for scientific discovery. Our solution
outperforms various baselines from the literature by a large margin.

##  2 Problem Definition

We first introduce the general active search model and notations. Suppose we
are given a finite set of points $\mathcal{X}\triangleq\\{x_{i}\\}$, which
includes a rare, valuable subset $\mathcal{R}\subset\mathcal{X}$. The members
of $\mathcal{R}$, which we call _targets_ or _positives_ , are not known _a
priori_ , but whether a given point $x\in\mathcal{X}$ is a target can be
determined by making a query to an oracle that returns the binary label
$y\triangleq\mathds{1}\\{x\in\mathcal{R}\\}$. We assume that querying the
oracle is expensive and that we only have a limited budget of $t$ queries to
do so. We denote a given dataset of queried points and their corresponding
labels as $\mathcal{D}=\\{(x_{i},y_{i})\\}$. At times, we will use
$\mathcal{D}_{i}$ to denote the dataset collected after $i$ queries.

###  2.1 Multifidelity Active Search

fidelity $L$fidelity $H$time1472225282356 8 23242627

Figure 1:  An illustration of our multifidelity active search model. Each
short vertical line indicates when a query finishes and the next query is
made. The numbers indicate the order in which queries are made across the two
oracles. The low-fidelity oracle is $k=2$ times faster than the high-fidelity
oracle. The budget on $H$ is $t=10$; in total, $T=t+kt-k=28$ queries are made.

In addition to the expensive oracle that returns the exact label of a query,
we have access to other cheaper but more noisy oracles that are low-fidelity
approximations to the exact oracle. We limit our setting to two levels of
fidelity—one exact oracle and one noisy oracle, denoted as $H$ and $L$,
respectively—but note that our proposed algorithm can be extended to more than
two fidelities with minor modifications, as we will discuss later.

For a given point $x\in\mathcal{X}$, we denote its exact label on $H$ as
$y_{H}$ and its noisy label on $L$ as $y_{L}$. Under this setting, a dataset
$\mathcal{D}$ can be partitioned into $\mathcal{D}_{L}$, the observations made
on $L$, and $\mathcal{D}_{H}$, those that are made on $H$. Recall that our
objective is to query as many targets as possible; as such, to express our
preference over different datasets, we use the natural utility function

$u\left(\mathcal{D}=\mathcal{D}_{L}\cup\mathcal{D}_{H}\right)\triangleq\\!\sum_{y_{i}\in\mathcal{D}_{H}}\\!y_{i},$

the number of targets discovered on the $H$ fidelity.

We assume that queries to different oracles are run _in parallel_ , but a
high-fidelity query takes longer to complete than a low-fidelity one. In
particular, we will assume that each query on $H$ is $k$ times slower than a
query on $L$.111This assumption is for simplicity; _asynchronous_ queries
could be addressed with a slight modification of our proposed algorithm. That
is, each time we make a query on $H$, we may make $k$ sequential queries on
$L$ while waiting for the result. The process will then repeat until the
budget is depleted. This model aims to emulate real-life scientific discovery
procedures, the main motivation for active search, where multiple fidelities
(e.g., computational and experimental campaigns) are often run in parallel but
vary in response time.

If we are given a budget of $t$ queries on $H$, we can make $kt-k$ queries on
$L$; in total, $T=t+kt-k$ queries are made. (Although a total of $kt$ queries
can be made on $L$, after the final query on $H$ is made at iteration $T$,
further queries on $L$ do not affect the final utility. We thus terminate
after $T$ queries.) This query schedule illustrated in Figure 1 for $t=10$ and
$k=2$.

###  2.2 The Bayesian Optimal Policy

We now derive the optimal (expected-case) policy using Bayesian decision
theory. First, we assume access to a probabilistic classification model that
computes the posterior probability that a point $x\in\mathcal{X}$ is a target
given a dataset $\mathcal{D}$, $\Pr\left(y_{H}=1\mid x,\mathcal{D}\right)$, as
well as the posterior probability that the same $x$ is a positive on $L$,
$\Pr\left(y_{L}=1\mid x,\mathcal{D}\right)$.

Suppose we are currently at iteration $i+1\leq T$, having observed the dataset
$\mathcal{D}_{i}$, and now need to make the next query, requesting the label
of an unlabeled point $x_{i+1}$. The Bayesian optimal policy selects the point
that maximizes the expected utility of the terminal dataset $\mathcal{D}_{T}$,
assuming future queries will too be made optimally:

$x_{i+1}^{*}=\operatorname*{arg\,max}_{x_{i+1}\in\mathcal{X}\setminus\mathcal{D}_{i}}\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i}\right)\mid
x_{i+1},\mathcal{D}_{i}\bigr{]}.$

If the query is on $H$, the expectation is taken over the posterior
distribution of the label of the putative query $x_{i+1}$; if the query is on
$L$, it is over the _joint_ distribution of the $L$ label of $x_{i+1}$ and the
label of the pending $H$ query (recall that we sequentially query $k$ points
on $L$ while waiting for an $H$ query that runs in parallel to finish).

To compute this expected utility, we follow the backward induction procedure
described in Bellman (1957). In the base case, we are at iteration $T$ and
need to make the last $H$ query and

$\displaystyle\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{T-1}\mid
x_{T},\mathcal{D}_{T-1}\right)\bigr{]}\hfill$
$\displaystyle\qquad=\sum_{y_{H}}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{T-1}\right)\Pr\left(y_{H}\mid
x_{T},\mathcal{D}_{T-1}\right)$ $\displaystyle\qquad=\Pr\left(y_{H}=1\mid
x_{T},\mathcal{D}_{T-1}\right).$

The optimal decision at this final step is therefore to greedily query the
point most likely to be a target, maximizing $\Pr\left(y_{H}=1\mid
x_{T},\mathcal{D}_{T-1}\right)$. For the second-to-last query, which is on
$L$, the expected utility is

$\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{T-2}\mid
x_{T-1},\mathcal{D}_{T-2}\right)\bigr{]}\\\
=\mathbb{E}\left[\max_{x_{T}}\Pr\left(y_{H}=1\mid
x_{T},\mathcal{D}_{T-1}\right)\right].$ (1)

This is the expected future reward to be collected at the next and final step,
which we have shown to be optimally the greedy query. Again, at this
iteration, the second-to-last $H$ query is still pending, so the expectation
is taken over the joint distribution of the label of that query and that of
the putative $L$ query. In general, we compute this expected utility at
iteration $i+1\leq T$ for an $L$ query recursively as

$\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i}\mid
x_{i+1},\mathcal{D}_{i}\right)\bigr{]}\\\
=\mathbb{E}\left[\max_{x_{i+2}}\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i+1}\right)\mid
x_{i+2},\mathcal{D}_{i+1}\bigr{]}\right].$ (2)

For an $H$ query, this expected utility may still be recursively computed in
the same manner but has a different expansion:

$\displaystyle\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i}\mid
x_{i+1},\mathcal{D}_{i}\right)\bigr{]}$
$\displaystyle\qquad=\Pr\left(y_{H}=1\mid x_{i+1},\mathcal{D}_{i}\right)+$
$\displaystyle\qquad\quad~{}~{}\mathbb{E}\left[\max_{x_{i+2}}\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i+1}\right)\mid
x_{i+2},\mathcal{D}_{i+1}\bigr{]}\right].$ (3)

The new term in (2.2), $\Pr\left(y_{H}=1\mid x_{i+1},\mathcal{D}_{i}\right)$,
accounts for the possibility that our running reward increases if
$x_{i+1}\in\mathcal{R}$, as we are querying on fidelity $H$. This is simply
the probability that $x_{i+1}$ is indeed a target. Also different from (2),
the expectation is taken over the distribution of the label $y_{H}$ of the
putative point $x_{i+1}$ only, as there is no pending query at this time. An
intuitive interpretation of the sum in (2.2) is the balance between
exploitation, the immediate reward $\Pr\left(y_{H}=1\mid
x_{i+1},\mathcal{D}_{i}\right)$, and exploration, the expected future reward
$u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i+1}\right)$ conditioned on this
putative query. Again, the exploitation term is not present in (2) when an $L$
query is made, as the query cannot possibly increase our utility immediately.

Perhaps unsurprisingly, computing the optimal policy is a daunting task: the
time complexity of computing the expectation term in (2) and (2.2) is
$\mathcal{O}\bigl{(}(4n)^{\ell}\bigr{)}$, where $\ell=T-i$ is the number of
remaining queries and $n$ is the number of unlabeled points. This optimal
policy is computationally intractable, and suboptimal approximations are
needed in practice. One natural solution is to shorten the lookahead horizon,
pretending there are only $\ell^{\prime}<T-i$ iterations remaining. This idea
constitutes _myopic_ policies, the most straightforward of which is the greedy
strategy that queries the point with the highest probability of being a target
in the single-fidelity setting, setting $\ell^{\prime}=1$. However, the design
of a greedy policy for an $L$ query is not obvious, as no immediate reward can
be obtained by making the query.

###  2.3 Hardness of Approximation

In addition to demonstrating the intractability of the analogous Bayesian
optimal policy in the classical single-fidelity setting, Jiang et al. (2017)
proved that _no_ polynomial-time policy can approximate the optimal policy (in
terms of the expected terminal utility) by _any_ constant factor. This was
done by constructing an adversarial family of active search problems featuring
“hidden treasures” that are provably difficult to uncover without exponential
work. By modifying details of this construction, the performance of any
polynomial-time policy can be made arbitrarily worse in expectation than that
of the optimal policy.

This hardness result naturally extends to our multifidelity model, as even
access to a _perfectly faithful_ low-fidelity oracle cannot aid a polynomial-
time active search policy in one of these adversarial examples. If we assume
positives on $L$ also count towards our utility, one such active search
problem with a faithful low-fidelity oracle reduces to a single-fidelity
problem with a budget on $H$ increased by a factor of $(k+1)$. Unless $k$ is
_exponential_ in the initial budget, any polynomial-time policy remains
incapable of approximating the optimal policy. Under our model, only positives
on $H$ count towards our utility, so the performance of any policy is further
reduced. We therefore obtain the same hardness of approximation result.
However, we can still reasonably aim to design efficient approximations to the
optimal policy that perform well in practice.

##  3 Efficient Nonmyopic Approximation

Our proposed algorithm is inspired by the ens policy, introduced by Jiang et
al. (2017) for the single-fidelity active search setting. ens offers an
efficient and nonmyopic approximation to the optimal policy by assuming that
after the putative query, all remaining budget will be spent _simultaneously
in one batch_. Under this heuristic, the optimal decision following the
putative query is to greedily construct the batch with points having the
highest probabilities. The expected utility of this batch is simply the sum of
these highest probabilities due to linearity of expectation and therefore can
be computed efficiently. An interesting interpretation by Jiang et al. (2017)
about ens is that it matches the optimal policy given that after the putative
query, the labels of the remaining unlabeled points are conditionally
independent.

###  3.1 One-stage Approximation

As a stepping stone to our proposed policy, consider the following adoption of
ens. We reapply the heuristic by assuming that after the putative query, which
can be on either $L$ or $H$, all remaining $H$ queries will be made
simultaneously in one batch. Again, the optimal choice for this batch is the
set of most promising points, which we will refer to as the greedy $H$ batch.
Using the summation-prime symbol $\sum\nolimits_{s}^{\prime}$ to denote the
sum of the top $s$ terms, we approximate the maximum expected utility of the
remaining portion of the search in (2) and (2.2) as

$\max_{x_{i+2}}\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i+1}\right)\mid
x_{i+2},\mathcal{D}_{i+1}\bigr{]}\\\
\approx\sum\nolimits_{\ell_{H}}^{\prime}\Pr\left(y_{H}=1\mid
x_{i+2},\mathcal{D}_{i+1}\right),$ (4)

where $\ell_{H}=\lfloor(T-i-1)/(k+1)\rfloor$ is the number of remaining $H$
queries. The resulting policy queries the point that maximizes the expected
utility in (2) and (2.2), using (4) as an approximation. This policy is
cognizant of the remaining budget on $H$ and performs well in our experiments
(see appendix). However, by assuming that the greedy $H$ batch will be queried
immediately following the putative point, the strategy fails to consider
future queries that could be made to the low-fidelity oracle in its lookahead;
this motivates the design of our proposed policy described below.

###  3.2 Two-stage Approximation

Our main contribution is a policy we call mf–ens, an efficient approximation
to the optimal policy that additionally accounts for future $L$ queries. As
above, we assume in our lookahead that after the putative query, our $H$
budget will be spent simultaneously. However, prior to committing to that
final batch, we assume we may make $\overline{k}\leq k$ additional $L$ queries
(also simultaneously) with the goal of improving the expected utility of the
final $H$ batch. To faithfully emulate our search model, we set $\overline{k}$
to be the number of $L$ queries remaining before the next $H$ query is
made.222 When making an $H$ query, $\overline{k}=k$; when making an $L$ query,
we subtract the number of $L$ queries since the pending $H$ query. We denote
this batch of exploratory $L$ queries as
$X_{L}\subset\mathcal{X}\setminus\mathcal{D}_{L}$ and the corresponding labels
as $Y_{L}\in\\{0,1\\}^{\overline{k}}$. In the language of _approximate dynamic
programming_ (Bertsekas, 1995), the policy described in the previous
subsection is a _one-stage rollout_ policy where the base policy selects the
optimal batch on $H$ following the putative query. mf–ens on the other hand is
a _two-stage rollout_ policy whose base policy first queries the $L$ batch
$X_{L}$, and upon observing $Y_{L}$, adaptively queries the updated optimal
$H$ batch.

How should we construct $X_{L}$ to best improve the greedy $H$ batch at the
second stage? One option would be to appeal to nonmyopic policies for _batch
active search_ (Jiang et al., 2018), but the best-promising batch policies
become prohibitively expensive in this context as we would need to construct a
new batch for _each_ putative query and label. In general, the conditional
dependence among the labels in the set $Y_{L}$ poses a computational challenge
in approximating the expected utility gained from querying a given $X_{L}$
batch.

Recall that in the single-fidelity setting, the ens policy is optimal if
labels become conditionally independent after the chosen point. Let us make a
similar assumption to ease computation: we assume that labels become
conditionally independent after the putative query, _except_ for the pair of
labels (on $H$ and $L$) corresponding to each point. That is, revealing the
label $y_{L}$ of a point $x$ is allowed to affect our belief about the
corresponding label $y_{H}$, but not the belief about any other label of any
other point. This structure allows for efficiently sharing information between
the fidelities and enables our multistage lookahead approach.333This is only
used in policy construction and not in inference!

We now aim to quantify the value of querying an unlabeled point on $L$ in
improving the final $H$ batch. Our solution is motivated by a heuristic search
for alternatives to the members of the greedy $H$ batch that is assembled
under the one-stage approximation. Given a putative query, we still construct
that same greedy $H$ batch. Then, for each candidate $x$ not in the greedy $H$
batch, we consider the expected marginal gain in utility of querying it on $L$
and modifying the membership of the $H$ batch in light of its newly revealed
label $y_{L}$.

Recall that observing $y_{L}$ only changes the distribution of $y_{H}$ of the
same $x$ under our assumption. Let $p_{\ast}$ be the lowest success
probability among the current $H$ batch and consider two cases. If
$\Pr\left(y_{H}=1\mid x,\mathcal{D}\right)$ _exceeds_ $p_{\ast}$ as $y_{L}$ is
revealed, we swap out the corresponding least-promising member of the batch
with $x$ and thus increase the final batch’s expected utility. Otherwise, we
do not modify the current batch. The value of querying $x$ on $L$ and
observing $y_{L}$, denoted as $v\left(x\mid\mathcal{D}\right)$, is then

$v\left(x\mid\mathcal{D}\right)\triangleq\mathbb{E}_{y_{L}}\bigl{[}\max\left(\Pr\left(y_{H}=1\mid
x,y_{L},\mathcal{D}\right)-p_{*},0\right)\bigr{]}.$

This score can be rapidly computed for most models under our conditional
independence assumption. With this value function in hand, we then greedily
construct $X_{L}$ with the candidates having the highest values. Another
computational benefit of label independence is that
$v\left(x\mid\mathcal{D}\right)$ automatically vanishes for points already
labeled on $H$, as querying its $L$ label does not affect the belief of our
model and thus cannot improve the $H$ batch. This reduces our search space in
computing the exploratory batch $X_{L}$.

We now proceed to the last step of the rollout procedure in mf–ens:
marginalizing over $Y_{L}$, the labels of $X_{L}$. As previously described,
for each possible value of $Y_{L}$, we approximate the optimal sequence of
queries following the putative one and $X_{L}$ with the updated greedy $H$
batch given the newly revealed labels $Y_{L}$:

$\max_{x_{i+1}}\mathbb{E}\bigl{[}u\left(\mathcal{D}_{T}\setminus\mathcal{D}_{i}\right)\mid
x_{i+1},\mathcal{D}_{i}\bigr{]}\\\
\approx\mathbb{E}_{Y_{L}}\left[\sum\nolimits_{\ell_{H}}^{\prime}\Pr\left(y_{H}=1\mid
x,X_{L},Y_{L},\mathcal{D}_{i}\right)\right].$ (5)

At each iteration, mf–ens queries the candidate maximizing the expected
utility in (2) and (2.2), as approximated by (5). As an extension of ens, our
approach is nonmyopic and aware of the remaining budget; we will demonstrate
the impact of this reasoning in our experiments. The policy also factors in
future $L$ queries, actively taking advantage of the ability to query the low-
fidelity oracle. We give the pseudocode for the policy in the appendix.

###  3.3 Extension to Multiple Low Fidelities

So far, we have assumed our model only consists of one exact oracle $H$ and
one noisy oracle $L$. To extend mf–ens to settings where there are multiple
noisy oracles $\\{L_{i}\\}$ that approximate $H$, potentially at different
levels of fidelity, we still aim to design each query to maximize the expected
utility on $H$, marginalizing future experiments on the $\\{L_{i}\\}$ oracles.
Once again limiting the conditional dependence among labels to those of the
same point, now between _each_ lower-fidelity $L_{i}$ and $H$, we identify a
batch of appropriate size for each $L_{i}$, marginalize their labels, update
the probabilities on $H$, and compute the approximate expected utility. When
there is only one low fidelity, this strategy reduces to the base version of
mf–ens we have presented here.

###  3.4 Implementation and Pruning

Active search requires a classification model computing a given point’s
success probability with an oracle. We extend the $k$-nearest neighbor
introduced by Garnett et al. (2012) to our multifidelity setting by allowing
information observed on $L$ to propagate to $H$. Specifically, when
calculating the probability that an unlabeled point $x\in\mathcal{X}$ is a
positive on $H$, $\Pr\left(y_{H}=1\mid x,\mathcal{D}\right)$, we take into
account the revealed labels of its nearest neighbors on both fidelities, as
well as its own $L$ label, $y_{L}$. Effectively, we treat each given point
$x\in\mathcal{X}$ as having two separate copies: one corresponding to its $H$
label, denoted as $x_{H}$, and one corresponding to its $L$ label, denoted as
$x_{L}$. Compared to the single-fidelity $k$-nearest neighbor, the set of
nearest neighbors of $x_{H}$ is now doubled to include both copies of its
original neighbors and its own copy on $L$, $x_{L}$. Formally, denote
$\textsc{{nn}}_{\text{single}}(x)$ as the original nearest neighbor set of
$x$. The nearest neighbor set of $x_{H}$ under our multifidelity predictive
model is

$\displaystyle\textsc{{nn}}(x_{H})\triangleq\\{x_{L}\\}$
$\displaystyle\cup\bigl{\\{}x^{\prime}_{H}:x^{\prime}\in\textsc{{nn}}_{\text{single}}(x)\bigr{\\}}$
$\displaystyle\cup\bigl{\\{}x^{\prime}_{L}:x^{\prime}\in\textsc{{nn}}_{\text{single}}(x)\bigr{\\}}.$

To account for the unknown accuracy of the low-fidelity oracle, we apply a
damping factor $q\in(0,1)$ to the weights of the neighbors on $L$; $q$ is
dynamically set at each iteration via maximum likelihood estimation.

This model performs well in practice, is nonparametric, and can be efficiently
updated in light of new data. The last feature is essential in allowing for
fast lookahead, the central component of our method. The time complexity of a
naive implementation of mf–ens is $\mathcal{O}\big{(}2^{k}\,n^{2}\log
n\big{)}$,444Note that we are assuming $k$ is small enough that $2^{k}$ is a
small constant; if this is not the case, we may approximate (5) by sampling
instead, replacing $2^{k}$ by $s$, the number of samples used. where $k$ is
the number of $L$ queries made for each $H$ query and $n$ is the size of the
candidate pool. The corresponding time complexity of the single-fidelity ens
is $\mathcal{O}\big{(}n^{2}\log n\big{)}$ (Jiang et al., 2017), to which
mf–ens adds a factor of $2^{k}$ from the exhaustive marginalization of the
exploratory $L$ labels. We may also take advantage of the implementation trick
developed by Jiang et al. (2017), which reduces the time complexity to
$\mathcal{O}\big{(}2^{k}\,n\,(n+m\log m)\big{)}$, where $m\ll n$ is the
maximum number of unlabeled points whose probabilities are affected by a newly
revealed label. The interested reader may refer to §3.2 of Jiang et al. (2017)
for more detail.

We also extend existing branch-and-bound pruning strategies to further reduce
the computation time of our policy at each search iteration. First, following
previous work (Garnett et al., 2012; Jiang et al., 2017), we establish an
upper bound of the score function that is the approximate expected utility
defined by (5). This allows us to eliminate candidates whose score upper
bounds are lower than the current best score we have found, as their actual
scores cannot possibly be the final best score. These upper bounds are
computationally cheap to evaluate, so applying this pruning check only adds a
trivial overhead to each search iteration. Further, we develop an extension of
this pruning strategy by making use of the fact that computing the score of a
point involves marginalizing over its unknown label. We thus apply similar
pruning checks at every step during this marginalization, which allows us to
identify and prune suboptimal candidates “on the fly,” avoiding any
unnecessary computation.

Concretely, denote by $f(x)$ the score of an unlabeled point
$x\in\mathcal{X}$, defined by (5) for mf–ens. Suppose $x$ has a posterior
probability of $\pi=\Pr\left(y=1\mid x,\mathcal{D}\right)$, where $y$ is the
label to be returned by the oracle we are currently querying. As previously
described, we compute $f(x)$ by calculating its _partial values_ while
marginalizing over $y$:

$f(x)=\pi\,f\left(x\mid y=1\right)+(1-\pi)\,f\left(x\mid y=0\right),$

where $f(x\mid y)$ is the partial score of $x$ according to (5) when
conditioned on a value of $y$. Suppose before computing either $f(x\mid y=0)$
or $f(x\mid y=1)$, we know these partial scores are upper bounded by
$\overline{u}(x)$ given a new positive label and $\underline{u}(x)$ given a
new negative label:

$f(x\mid y=1)<\overline{u}(x);\quad f(x\mid y=0)<\underline{u}(x).$

As such, $f(x)$ need not be evaluated if

$\pi\,\overline{u}(x)+(1-\pi)\,\underline{u}(x)<f^{*},$ (6)

where $f^{*}$ is the current best score we have found. This pruning strategy
has been found to offer a significant speedup in previous work (Garnett et
al., 2012; Jiang et al., 2017, 2018).

We extend this strategy by considering the case in which a given candidate $x$
is not pruned because (6) is not satisfied, and we proceed with the
calculation of $f(x)$. Now, suppose we have computed only $f\left(x\mid
y=1\right)$ and not yet $f\left(x\mid y=0\right)$ and observe that

$\pi\,f\left(x\mid y=1\right)+(1-\pi)\,\underline{u}(x)<f^{*},$

then we may also safely conclude that $f(x)$ cannot possibly exceed $f^{*}$
without needing to go further and compute $f(x\mid y=0)$. If this condition is
met, we simply abort the computation of the current score $f(x)$ and move on
to the next unpruned candidate. For each $H$ query, we apply this partial
pruning check once (either after conditioning on $y_{H}=1$ and before on
$y_{H}=0$ or vice versa) for each candidate that is not eliminated by the full
pruning check (6). For an $L$ query, we may do this at most three times for
each candidate, as the marginalization over the joint distribution of the
putative label and the pending $H$ label when computing $f(x)$ involves four
different possible label combinations. We quantify the effectiveness of these
pruning strategies and show that they can significantly reduce the computation
time of our methods in the appendix.

At an iteration where the current best score $f^{*}$ does not exceed the
majority of the score upper bounds, many candidates may be left unpruned. In
order to help our policy avoid having to calculate the scores of a large
number of candidates and consequently spending too much time on a single
iteration, we place an upper limit on how many candidates are to be considered
before we terminate the search. Our approach is to first follow the lazy-
evaluation strategy introduced by Jiang et al. (2018) and sort the candidates
by their score upper bounds. With this sorting, candidates with higher upper
bounds will be considered first, and a candidate will never be considered if
it will be pruned later on. Now, at each iteration, if after having considered
$u$ candidates and noticing that there are still unpruned points remaining, we
simply consider a randomly selected subset of size at most $s$ of the unpruned
set, before terminating the search and returning the current best candidate.
In our experiments, we set $u=s=500$ for mf–ens.

We note that this strategy is only used when pruning fails to reduce our
search space to be below $u+s$, and to allow us to collect results over a long
horizon over many repeated experiments. When applied in a real-life planning
setting, mf–ens can still cover the entirety of a large pool of candidates if
desired, even with a significant portion of the pool unpruned. In our
experiments, mf–ens takes approximately 30 seconds to reach its quota of 1000
candidates when pruning is unsuccessful; the time for it to fully cover a pool
of 100 000 points (about the size of the real-world datasets used in our
experiments) is thus well under one hour. In short, our policy remains
tractable in real-life settings, even without successful pruning.

##  4 Related Work

Table 1:  Experiment results with an $H$ budget of $t=300$, averaged across
all repeated experiments for each setting. $\theta$ is the simulated false
positive rate of the low-fidelity oracle; $k$ is the number of $L$ queries
made between two $H$ queries (i.e., the speed ratio between the two
fidelities). Each entry denotes the average number of targets found across the
repeated experiments and the corresponding standard error in parentheses. The
best performance in each column is highlighted bold.

|  ecfp4 |  gpidaph3 | bmg
---|---|---|---
| $\theta=0.1$ | $\theta=0.3$ | $\theta=0.1$ | $\theta=0.3$ | $\theta=0.1$ | $\theta=0.3$
| $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$
ens | 218 (4.5) | 218 (4.5) | 213 (4.2) | 213 (4.2) | 197 (4.7) | 197 (4.7) | 186 (5.3) | 186 (5.3) | 279 (1.6) | 279 (1.6) | 278 (1.7) | 278 (1.7)
mf–ucb | 227 (4.8) | 238 (4.4) | 200 (5.7) | 206 (5.8) | 200 (5.7) | 212 (5.7) | 187 (5.5) | 200 (5.3) | 284 (1.0) | 289 (1.0) | 274 (2.0) | 277 (1.7)
ug | 228 (4.9) | 237 (4.4) | 207 (5.8) | 209 (5.8) | 204 (5.6) | 211 (5.7) | 191 (5.6) | 202 (5.6) | 283 (1.2) | 287 (1.1) | 278 (2.0) | 280 (1.6)
mf–ens | 244 (3.5) | 250 (3.2) | 226 (4.5) | 230 (4.4) | 230 (3.8) | 243 (3.4) | 208 (5.0) | 221 (4.3) | 290 (0.7) | 294 (0.6) | 286 (1.1) | 286 (1.1)

This work is an extension of the larger active search paradigm, first
introduced by Garnett et al. (2012). Active search is a variant of _active
learning_ (Settles, 2009) where the goal is not to learn an accurate model but
to find and query positive labels. Previous work has studied active search
under different settings such as finding a given number of targets as
fast/cheaply as possible (Warmuth et al., 2002, 2003; Jiang et al., 2019),
making queries in batches (Jiang et al., 2018), or when points have real-
valued utility (Vanchinathan et al., 2015). Our work generalizes ens, the
policy proposed by Jiang et al. (2017) from the single-fidelity setting. The
authors of that work demonstrated that their policy is nonmyopic and aware of
its exact budget, allowing it to automatically balance between exploration and
exploitation during search and outperform various baselines by a large margin.
We will make the same observations about our policy.

Multifidelity active search was first examined by Klyuchnikov et al. (2019),
who specifically considered the search problem of a recommender system:
identifying items that users of a given application are interested in. They
modeled predictions made by a trained preference model as output of a low-
fidelity oracle and proposed a co-kriging predictive model (Álvarez et al.,
2012) to perform inference on the users’ true preferences. Under their
setting, queries to the oracles are made sequentially, one after another. Our
multifidelity setting is different, modeling situations where oracles of
different fidelities are available to run in parallel and vary in their
response times, common in scientific experiments and testing. Regardless, our
proposed policy could be naturally adopted to their sequential model; the only
difference in the computation is that the marginalization over a pending $H$
label is no longer necessary. Further, as we will show in later experiments,
our algorithm outperforms the adoption of their upper confidence bound (ucb)
policy for various datasets. To our knowledge, our work is the first to tackle
multifidelity active search using Bayesian decision theory.

Active search is equivalent to _Bayesian optimization_ (bo) (Brochu et al.,
2010; Snoek et al., 2012) with binary observations and cumulative reward.
Multifidelity bo itself has been studied considerably, and policies
corresponding to common acquisition functions have been adopted to
multifidelity settings, including expected improvement (Huang et al., 2006;
Picheny et al., 2013), knowledge gradient (Poloczek et al., 2017; Wu et al.,
2020), and ucb (Kandasamy et al., 2017). However, most of these policies are
derived from or motivated by greedy approximations to the optimal policy under
different utility functions, and to our knowledge, no _nonmyopic_
multifidelity bo policies have been proposed.

Active search is related to the _multi-armed bandit_ (mab) problem (Lai &
Robbins, 1985). In particular, querying the label of a point can be viewed as
“pulling an arm” in mab; in active search, an arm cannot be pulled twice but
is correlated to its neighbors. Kandasamy et al. (2016) studied a formulation
of multifidelity mab in which each arm may be pulled on different fidelities
at different costs, and proposed a policy that is a variant of ucb. By
assuming the accuracy of the lower fidelities is known, the authors derived
strong theoretical guarantees for their proposed policy.

##  5 Experiments

$0$$50$$100$$150$$200$$250$$300$$0$$10$$20$number of $H$ queriesdifference in
utilitymean difference95% ci Figure 2: The difference in cumulative targets
found between mf–ens and ens, averaged across all experiments and datasets.

$0$$50$$100$$150$$200$$250$$300$$0.6$$0.8$number of $H$ queriesqueried
probabilitymf-ensmf-ucbug

$0$$50$$100$$150$$200$$250$$300$$0$$10$$20$number of $H$ queriesdifference in
utilitymean difference95% ci

Figure 3:  An illustration of budget-awareness exhibited by mf–ens. Left: The
average progressive probabilities of points queried on $H$ by different active
search policies. Right: The difference in cumulative targets found between
mf–ens and mf–ucb. The results are averaged across all experiments and
datasets.

We now compare the empirical performance of mf–ens against several
benchmarks.555 Matlab implementations of our policies are available at:
https://github.com/KrisNguyen135/multifidelity-active-search . As previously
described, ens is a state-of-the-art, nonmyopic active search policy in the
single-fidelity setting. In our experiments, this policy simply ignores the
low-fidelity oracle, serving as a single-fidelity baseline to illustrate the
benefit of having access to low-fidelity queries. We also test against the
mf–ucb policy for multifidelity active search, recently proposed by
Klyuchnikov et al. (2019). Under this policy, each candidate $x$ has a
ucb–style score of
$\alpha\left(x,\mathcal{D}\right)=\pi+\beta\sqrt{\pi(1-\pi)}$, where $\pi$ is
the probability of $x$ having a positive label on the fidelity being queried
and $\beta$ is the exploration/exploitation tradeoff parameter.666The authors
also considered an odd scenario in which the correlation between the two
oracles is _negative._ We assume that this correlation is always positive, and
it is constructed to be so. We set $\beta=0.01$ for $L$ queries and
$\beta=0.001$ for $H$ queries, as suggested in the same work. For another
benchmark, we consider a simple but natural heuristic for multifidelity
optimization, that low-fidelity queries should serve to narrow down the most-
promising search regions (exploration), so that more informed queries could be
made on the higher-fidelity oracle (exploitation). Inspired by this heuristic,
we design a policy we call ug (for _uncertainty_ and _greedy_ sampling) that
always queries the most uncertain points on $L$ and the points most likely to
be positive on $H$. Uncertainty sampling is chosen for the role of exploration
due to its popularity as an active learning technique (Lewis & Gale, 1994). ug
may also be viewed as the limit of ucb when $\beta$ approaches infinity on $L$
for maximum exploration and 0 on $H$ for maximum exploitation.

Datasets for multifidelity active search are not readily available due to the
relative novelty of the problem setting, at least not in our motivating area
of scientific discovery. However, numerous high-quality datasets are available
in the single-fidelity setting, which we will adopt and use to simulate
multifidelity search. Namely, for each dataset, we simulate the noisy labels
returned by the low-fidelity oracle using the following procedure. We first
create a duplicate of the true labels. Given $\theta\in(0,1)$, we randomly
select a fraction $\theta$ of the positives from this duplicate set and “flip”
their labels to negative. We also randomly select the same number of negatives
and flip their values to positive. This perturbed set of labels is then used
as the $L$ labels. This construction yields simulated $L$ labels with a false
positive rate of $\theta$ and a false negative rate of $\theta\,r/(1-r)$,
where $r\in(0,1)$ is the prevalence rate of the positive set $\mathcal{R}$ in
$\mathcal{X}$. In a typical active search problem, $\mathcal{R}$ is rare and
$r\ll 1$, making the false negative rate much lower than the false positive
rate, a common characteristic of many real-world scientific discovery and
testing procedures.

We set $\theta\in\\{0.1,0.3\\}$. We set $k$, the number of $L$ queries that
are made for each $H$ query, to be either 2 or 5, and set the budget on $H$ to
be $300$. In each experiment, a policy starts with an initial training dataset
of one randomly selected target whose $L$ label is also positive.

###  5.1 Datasets

We conducted experiments on three real-world scientific discovery datasets
used in previous studies (Jiang et al., 2017, 2018, 2019). The first two come
from drug discovery, where the goal is to discover chemical compounds
exhibiting binding activity with a given protein. Each protein defines the
target for an active search problem. Here we used the first 50 proteins from
the bindingdb database (Liu et al., 2007) described by Jiang et al. (2017). A
set of 100 000 compounds sampled from the zinc database (Sterling & Irwin,
2015) served as a shared negative set. Features for the compounds are binary
vectors encoding chemical characteristics, also known as a chemoinformatic
fingerprint. We considered two fingerprints delivering good performance in
previous studies: ecfp4 and gpidaph3\. The size of the positive set for these
50 targets ranged from 205 to 1488 (mean 538), having an average prevalence
rate of $r\approx 0.5\%$. For each of these two datasets and 50 targets, we
repeat each experiment five times, for a total of 500 search simulations.

The other dataset is related to a materials science application. The targets
in this case are alloys that can form bulk metallic glasses (bmgs), which have
higher toughness and better wear resistance than crystalline alloys. This
dataset comprises 106 810 alloys from the materials literature (Kawazoe et
al., 1997; Ward et al., 2016), 4275 of which exhibit glass-forming ability
($r\approx 4\%$). We repeated each experiment 50 times for this dataset.

We report the average number of targets found by each policy across the
experiments with standard errors in parentheses in Table 2; each column
corresponds to a specific setting of $\theta$ and $k$ under a dataset. We
observe that our policy mf–ens outperforms all baselines by a large margin. In
each column, a two-sided paired $t$–test rejects the hypothesis that the
average difference in the number of targets found between mf–ens and any
baseline is zero with overwhelming confidence, returning a $p$–value of at
most $4\times 10^{-5}$. We report these $p$–values in the appendix. Finally,
looking across the columns, we notice the expected trends: performance of all
algorithms except for ens, which does not utilize fidelity $L$, improves with
higher $k$ (when $L$ is cheaper) or with lower $\theta$ (when $L$ is more
accurate).

###  5.2 Performance Gain from Multifidelity Search

To further examine the benefit of having access to more than just the exact
oracle, we visualize the difference in the cumulative number of targets found
between mf–ens and the single-fidelity policy ens, averaged across all
experiments, in Figure 2. We observe that mf–ens completely dominates ens,
finding roughly linearly more targets throughout the search. This large
difference in empirical performance illustrates the usefulness of the
simulated low-fidelity oracle in our experiments, and suggests that our
approach is likely to benefit search with any budget.

###  5.3 Nonmyopic Behavior

We have claimed that mf–ens is nonmyopic and aware of its remaining budget at
any given time during a search. We demonstrate this nonmyopia by first
comparing the progressive probabilities of the points queried on $H$ (at the
time of the queries being made) by the policy against the myopic baselines
mf–ucb and ug, averaged across all experiments, in the left panel of Figure 3.
Initially, mf–ens chooses points with lower probabilities, exploring the
space. As the search progresses, mf–ens queries more promising points,
smoothly transitioning to exploitation. The opposite trend can be observed for
ucb and ug, whose probabilities decrease over time due to greedy behavior.
This difference translates to distinct patterns in the cumulative reward
achieved by these policies. The right panel of Figure 3 shows the difference
in the cumulative number of targets found between mf–ens and mf–ucb, also
averaged across all experiments. During the first half of the search, mf–ens
appears to perform _worse_ than mf–ucb, but quickly recovers and outperforms
the latter in the end. The corresponding plot comparing mf–ens and ug shows a
similar trend and is deferred to the appendix.

Overall, this phenomenon perfectly highlights the automatic tradeoff between
exploration and exploitation exhibited by mf–ens: the policy makes its initial
queries to explore the search space, often requesting labels that are not the
most likely to be positive and failing to collect substantial immediate
reward; however, as the budget decreases, its queries grow more exploitative
and are ultimately more successful than those from myopic policies by
leveraging what it has learned.

##  6 Conclusion

We have proposed a multifidelity active search model in which an exact oracle
and a cheaper surrogate are queried in parallel. We presented a novel
nonmyopic policy (based on two-stage rollout) for this setting that reasons
about the remaining queries on both fidelities seeking to maximize the total
number of discoveries. Our policy is aware of the remaining budget and
dynamically balances exploration and exploitation. Experiments on real-world
data demonstrate that the policy significantly outperforms myopic benchmarks.

## Acknowledgements

We would like to thank the anonymous reviewers for their feedback. qn and rg
were supported by the National Science Foundation (nsf) under award numbers
oac–1940224 and iis–1845434.

## References

  * Álvarez et al. (2012) Álvarez, M. A., Rosasco, L., and Lawrence, N. D.  Kernels for Vector-Valued Functions: A Review .  _Foundations and Trends in Machine Learning_ , 4(3):195–266, 2012.
  * Bellman (1957) Bellman, R.  _Dynamic Programming_.  Princeton University Press, 1957.
  * Bertsekas (1995) Bertsekas, D. P.  _Dynamic Programming and Optimal Control_ , volume 1.  Athena Scientific Belmont, MA, 1995.
  * Brochu et al. (2010) Brochu, E., Cora, V. M., and De Freitas, N.  A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning.  2010\.  arXiv preprint arXiv:1012.2599 [cs.LG].
  * Garnett et al. (2012) Garnett, R., Krishnamurthy, Y., Xiong, X., Schneider, J., and Mann, R.  Bayesian Optimal Active Search and Surveying.  In _Proceedings of the 29th International Conference on Machine Learning_ , 2012.
  * Huang et al. (2006) Huang, D., Allen, T. T., Notz, W. I., and Miller, R. A.  Sequential kriging optimization using multiple-fidelity evaluations.  _Structural and Multidisciplinary Optimization_ , 32(5):369–382, 2006.
  * Jiang et al. (2017) Jiang, S., Malkomes, G., Converse, G., Shofner, A., Moseley, B., and Garnett, R.  Efficient Nonmyopic Active Search.  In _Proceedings of the 34th International Conference on Machine Learning_ , pp. 1714–1723, 2017.
  * Jiang et al. (2018) Jiang, S., Malkomes, G., Abbott, M., Moseley, B., and Garnett, R.  Efficient nonmyopic batch active search.  In _Advances in Neural Information Processing Systems 31_ , pp. 1099–1109, 2018.
  * Jiang et al. (2019) Jiang, S., Garnett, R., and Moseley, B.  Cost effective active search.  In _Advances in Neural Information Processing Systems 32_ , pp. 4880–4889, 2019.
  * Kandasamy et al. (2016) Kandasamy, K., Dasarathy, G., Schneider, J., and Póczos, B.  The Multi-fidelity Multi-armed Bandit.  In _Advances in Neural Information Processing Systems 29_ , pp. 1785–1793, 2016.
  * Kandasamy et al. (2017) Kandasamy, K., Dasarathy, G., Schneider, J., and Póczos, B.  Multi-fidelity Bayesian Optimisation with Continuous Approximations.  In _Proceedings of the 34th International Conference on Machine Learning_ , pp. 2861–2878, 2017.
  * Kawazoe et al. (1997) Kawazoe, Y., Yu, J.-Z., Tsai, A.-P., and Masumoto, T. (eds.).  _Nonequilibrium Phase Diagrams of Ternary Amorphous Alloys_ , volume 37A of _Condensed Matters_.  Springer-Verlag, 1997.
  * Klyuchnikov et al. (2019) Klyuchnikov, N., Mottin, D., Koutrika, G., Müller, E., and Karras, P.  Figuring out the User in a Few Steps: Bayesian Multifidelity Active Search with Cokriging.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pp. 686–695, 2019.
  * Lai & Robbins (1985) Lai, T. L. and Robbins, H.  Asymptotically Efficient Adaptive Allocation Rules.  _Advances in Applied Mathematics_ , 6(1):4–22, 1985.
  * Lewis & Gale (1994) Lewis, D. D. and Gale, W. A.  A Sequential Algorithm for Training Text Classifiers.  In _Proceedings of the 17th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pp. 3–12, 1994.
  * Liu et al. (2007) Liu, T., Lin, Y., Wen, X., Jorissen, R. N., and Gilson, M. K.  BindingDB: A web-accessible database of experimentally determined protein–ligand binding affinities.  _Nucleic acids research_ , 35(suppl_1):D198–D201, 2007.
  * Picheny et al. (2013) Picheny, V., Ginsbourger, D., Richet, Y., and Caplin, G.  Quantile-Based Optimization of Noisy Computer Experiments with Tunable Precision.  _Technometrics_ , 55(1):2–13, 2013.
  * Poloczek et al. (2017) Poloczek, M., Wang, J., and Frazier, P.  Multi-Information Source Optimization.  In _Advances in Neural Information Processing Systems 30_ , pp. 4288–4298, 2017.
  * Settles (2009) Settles, B.  Active Learning Literature Survey.  Technical report, University of Wisconsin-Madison Department of Computer Sciences, 2009.
  * Snoek et al. (2012) Snoek, J., Larochelle, H., and Adams, R. P.  Practical Bayesian Optimization of Machine Learning Algorithms.  In _Advances in Neural Information Processing Systems 25_ , pp. 2951–2959, 2012.
  * Sterling & Irwin (2015) Sterling, T. and Irwin, J. J.  Zinc 15–Ligand Discovery for Everyone.  _Journal of Chemical Information and Modeling_ , 55(11):2324–2337, 2015.
  * Vanchinathan et al. (2015) Vanchinathan, H. P., Marfurt, A., Robelin, C.-A., Kossmann, D., and Krause, A.  Discovering Valuable Items from Massive Data.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pp. 1195–1204, 2015.
  * Ward et al. (2016) Ward, L., Agrawal, A., Choudhary, A., and Wolverton, C.  A general-purpose machine learning framework for predicting properties of inorganic materials.  _npj Computational Materials_ , 2(1):1–7, 2016\.
  * Warmuth et al. (2002) Warmuth, M. K., Rätsch, G., Mathieson, M., Liao, J., and Lemmen, C.  Active learning in the drug discovery process.  In _Advances in Neural Information Processing Systems 15_ , pp. 1449–1456, 2002.
  * Warmuth et al. (2003) Warmuth, M. K., Liao, J., Rätsch, G., Mathieson, M., Putta, S., and Lemmen, C.  Active Learning with Support Vector Machines in the Drug Discovery Process.  _Journal of Chemical Information and Computer Sciences_ , 43(2):667–673, 2003.
  * Wu et al. (2020) Wu, J., Toscano-Palmerin, S., Frazier, P. I., and Wilson, A. G.  Practical Multi-Fidelity Bayesian Optimization for Hyperparameter Tuning.  In _Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence_ , pp. 788–798, 2020.

##  Appendix A Pseudocode

We give the pseudocode for mf–ens in Algorithm 1 (for $H$ queries) and
Algorithm 2 (for $L$ queries). The two versions are nearly identical; however,
in Algorithm 2, we marginalize both the putative label $y_{L}$ and the pending
label $y_{H}$ for $L$ queries, and the $\Pr(y_{H}=1\mid x,\mathcal{D})$ term
is dropped from the final line according to (1) in the main paper.

##  Appendix B The One-Stage Approximation Policy

We name the one-stage approximation policy described in §3.1 of the main paper
$H$–ens. Recall $H$–ens applies the ens heuristic directly by assuming in its
lookahead that immediately after the putative query, all budget on $H$ is
spent simultaneously on a batch of queries, disregarding any $L$ queries that
may be made before that batch.

The time complexity of $H$–ens is the same as that of ens,
$\mathcal{O}\bigl{(}n\,(\log n+m\log m)\bigr{)}$, where $n$ is the size of the
candidate pool and $m$ is the maximum number of points whose probablities are
affected by a newly revealed label (Jiang et al., 2017). To avoid having to
consider a large number of unpruned candidates on a single iteration, we set
$u=s=1000$, where $u$ is the upper limit of the number of candidates are to be
considered before a random subset of size $s$ of the remaining candidates is
selected. This strategy is described in detail in §3.3 of the main paper.

Algorithm 1 mf–ens for $H$ queries

$x$, $\mathcal{D}$approximate expected utility for querying $x$on $H$,
$f(x)$design query by maximizing $y_{H}\in\\{-,+\\}$probabilities:
$\Pr(y_{H}\mid x,\mathcal{D})$determine $X_{L}$by finding top-$k$$v$scores
among unlabeled $L$points$Y_{L}\in\\{-,+\\}^{k}$probabilities: $\Pr(Y_{L}\mid
x,y_{H},X_{L},\mathcal{D})$$s\left(x,y_{H},X_{L},Y_{L}\right)\leftarrow\sum\nolimits_{\ell_{H}}^{\prime}\Pr\left(y^{\prime}_{H}=1\mid
x^{\prime},x,y_{H},X_{L},Y_{L}\right)$(4) $f\left(x\mid
y_{H}\right)\leftarrow\mathbb{E}_{Y_{L}}\big{[}s\left(x,y_{H},X_{L},Y_{L}\right)\mid
x,y_{H},X_{L}\big{]}$(4) $f(x)\leftarrow\Pr\left(y_{H}=1\mid
x,\mathcal{D}\right)+\mathbb{E}_{y_{H}}\big{[}f\left(x\mid y_{H}\right)\mid
x,\mathcal{D}\big{]}$(2)

Algorithm 2 mf–ens for $L$ queries

$x$, $\mathcal{D}$approximate expected utility for querying $x$on $L$,
$f(x)$design query by maximizing $y_{L},y_{H}\in\\{-,+\\}^{2}$probabilities:
$\Pr(y\mid x,\mathcal{D})$determine $X_{L}$by finding
top-$\overline{k}$$v$scores among unlabeled
$L$points$Y_{L}\in\\{-,+\\}^{\overline{k}}$probabilities: $\Pr(Y_{L}\mid
x,y_{L},X_{L},\mathcal{D})$$s\left(x,y_{L},X_{L},Y_{L}\right)\leftarrow\sum\nolimits_{\ell_{H}}^{\prime}\Pr\left(y^{\prime}_{H}=1\mid
x^{\prime},x,y_{L},X_{L},Y_{L}\right)$(4) $f\left(x\mid
y_{L}\right)\leftarrow\mathbb{E}_{Y_{L}}\big{[}s\left(x,y_{L},X_{L},Y_{L}\right)\mid
x,y_{L},X_{L}\big{]}$(4) $f(x)\leftarrow\mathbb{E}_{y_{L}}\big{[}f\left(x\mid
y_{L}\right)\mid x,\mathcal{D}\big{]}$(1)

##  Appendix C Further Experiment Results

Table 2:  Experiment results with an $H$ budget of $t=300$, averaged across
all repeated experiments for each setting. $\theta$ is the simulated false
positive rate of the low-fidelity oracle; $k$ is the number of $L$ queries
made between two $H$ queries (i.e., the speed ratio between the two
fidelities). Each entry denotes the average number of targets found across the
repeated experiments and the corresponding standard error in parentheses. The
best performance in each column is highlighted bold; those that are not
significantly worse than the best (using a one-sided paired $t$-test with
significance level $\alpha=0.05$) are in blue italic.

|  ecfp4 |  gpidaph3 | bmg
---|---|---|---
| $\theta=0.1$ | $\theta=0.3$ | $\theta=0.1$ | $\theta=0.3$ | $\theta=0.1$ | $\theta=0.3$
| $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$
ens | 218 (4.5) | 218 (4.5) | 213 (4.2) | 213 (4.2) | 197 (4.7) | 197 (4.7) | 186 (5.3) | 186 (5.3) | 279 (1.6) | 279 (1.6) | 278 (1.7) | 278 (1.7)
mf–ucb | 227 (4.8) | 238 (4.4) | 200 (5.7) | 206 (5.8) | 200 (5.7) | 212 (5.7) | 187 (5.5) | 200 (5.3) | 284 (1.0) | 289 (1.0) | 274 (2.0) | 277 (1.7)
ug | 228 (4.9) | 237 (4.4) | 207 (5.8) | 209 (5.8) | 204 (5.6) | 211 (5.7) | 191 (5.6) | 202 (5.6) | 283 (1.2) | 287 (1.1) | 278 (2.0) | 280 (1.6)
$H$–ens | _241 (3.6)_ | _249 (3.5)_ | _223 (4.3)_ | 231 (4.1) | 221 (4.5) | 245 (3.0) | 202 (4.8) | 222 (3.8) | 290 (0.7) | 294 (0.5) | _284 (1.2)_ | _285 (1.2)_
mf–ens | 244 (3.5) | 250 (3.2) | 226 (4.5) | _230 (4.4)_ | 230 (3.8) | _243 (3.4)_ | 208 (5.0) | _221 (4.3)_ | _290 (0.7)_ | _294 (0.6)_ | 286 (1.1) | 286 (1.1)

Table 2 shows the average number of targets found by all considered policies
including $H$–ens. $H$–ens performs well in comparison with myopic baselines
and the single-fidelity ens policy. Against $H$–ens, our main policy mf–ens
offers a slight improvement. Across all experiments, mf–ens finds an average
of 237 targets, whereas $H$–ens finds 235; a two-sided paired $t$–test rejects
the hypothesis that the average difference in the number of targets found
between the two policies is zero with a $p$–value of $p=0.044$.

Table 3:  $p$–values from two-paired $t$–tests, each testing for the
hypothesis that the average number of targets found between mf–ens and another
policy is zero.

|  ecfp4 |  gpidaph3 | bmg
---|---|---|---
| $\theta=0.1$ | $\theta=0.3$ | $\theta=0.1$ | $\theta=0.3$ | $\theta=0.1$ | $\theta=0.3$
| $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$ | $k=2$ | $k=5$
ens | $4\times 10^{-27}$ | $3\times 10^{-30}$ | $7\times 10^{-12}$ | $2\times 10^{-18}$ | $1\times 10^{-26}$ | $1\times 10^{-40}$ | $3\times 10^{-14}$ | $3\times 10^{-22}$ | $4\times 10^{-10}$ | $5\times 10^{-14}$ | $7\times 10^{-7}$ | $6\times 10^{-5}$
mf–ucb | $1\times 10^{-10}$ | $1\times 10^{-6}$ | $7\times 10^{-18}$ | $6\times 10^{-14}$ | $1\times 10^{-18}$ | $6\times 10^{-18}$ | $1\times 10^{-13}$ | $7\times 10^{-18}$ | $5\times 10^{-8}$ | $8\times 10^{-6}$ | $1\times 10^{-8}$ | $8\times 10^{-8}$
ug | $5\times 10^{-8}$ | $2\times 10^{-7}$ | $2\times 10^{-10}$ | $2\times 10^{-10}$ | $2\times 10^{-15}$ | $5\times 10^{-19}$ | $6\times 10^{-11}$ | $1\times 10^{-14}$ | $9\times 10^{-8}$ | $4\times 10^{-8}$ | $1\times 10^{-5}$ | $4\times 10^{-5}$

Table 3 shows the $p$–values from two-sided paired $t$–tests, each testing for
the hypothesis that the average number of targets found between mf–ens and
each of the baselines considered in §5 of the main paper is zero. All
hypotheses are rejected with overwhelming confidence.

We also show the difference in the cumulative number of targets found between
mf–ens and ug in Figure 4. The plot shows the same trend as that for mf–ens
and mf-ucb in the main paper: mf–ens underperforms the myopic policy at the
beginning of the search during its exploration phase, but quickly recovers and
is able to find significantly more targets in the end.

$0$$50$$100$$150$$200$$250$$300$$0$$10$$20$number of $H$ queriesdifference in
utilitymean difference95% ci Figure 4: The difference in cumulative targets
found between mf–ens and ug, averaged across all experiments and datasets.
Table 4:  Average pruning rate and the contribution of each of the two pruning
methods.

policies | full coverage rate | given full coverage
---|---|---
total prune % | partial prune %
$H$–ens | 45.2% | 85.2% | 14.0%
mf–ens | 27.6% | 89.4% | 9.2%

###  C.1 Effects of Pruning

Finally, we examine the effectiveness of our pruning strategies in helping
$H$–ens and mf–ens cover the entire search space. Each of the two policies has
an upper limit on how many candidates are to have their scores fully computed.
If this limit is reached, we only consider a random subset of the remaining
candidate pool. We classify each time a policy returns while there are
unpruned candidates remaining as an instance of failure to cover the entire
space.

First, we compute the fraction of iterations, across all experiments, in which
this limit is not exceeded, or in other words, how often each policy can
exactly consider all candidates with the help of pruning. When a point is
pruned, it may be before any actual computation (total pruning) if its score
upper bound is lower than the current highest score $f^{*}$, satisfying (5) in
the main paper. Otherwise, it may be pruned during the calculation of its
scores (partial pruning) if its partial score, combined with partial upper
bounds, is lower than $f^{*}$. We keep track of the fraction of points pruned
by each of the two methods; these results are averaged across the iterations
in which pruning helps cover the entire candidate pool. We report these
statistics in Table 4.

We observe that while we do not cover the entire candidate pool in many
iterations, the effect of pruning is dramatic when it is successful. The
existing pruning method (total pruning) helps eliminate most of the
candidates, and our extension (partial pruning) raises the combined pruning
rate to roughly 99% on average. In our experiments, successful pruning could
reduce the time a policy takes to produce a decision from under an hour to
mere seconds.


# Machine Learning the Dimension of a Polytope

Tom Coates  Department of Mathematics
Imperial College London
180 Queen’s Gate
London
SW7 2AZ
UK<EMAIL_ADDRESS>,  Johannes Hofscheier  School of Mathematical
Sciences
University of Nottingham
Nottingham
NG7 2RD
UK<EMAIL_ADDRESS>and  Alexander M. Kasprzyk  School of
Mathematical Sciences
University of Nottingham
Nottingham
NG7 2RD
UK<EMAIL_ADDRESS>
###### Abstract.

We use machine learning to predict the dimension of a lattice polytope
directly from its Ehrhart series. This is highly effective, achieving almost
$100\%$ accuracy. We also use machine learning to recover the volume of a
lattice polytope from its Ehrhart series, and to recover the dimension,
volume, and quasi-period of a rational polytope from its Ehrhart series. In
each case we achieve very high accuracy, and we propose mathematical
explanations for why this should be so.

###### Key words and phrases:

Machine learning, convex polytope, Ehrhart series, quasi-period.

###### 2020 Mathematics Subject Classification:

52B20, 68T05 (Primary); 14M25 (Secondary)

##  1\. Introduction

Let $P\subset\mathbb{Z}^{d}\otimes_{\mathbb{Z}}\mathbb{Q}$ be a convex lattice
polytope of dimension $d$, that is, let $P$ be the convex hull of finitely
many points in $\mathbb{Z}^{d}$ whose $\mathbb{Q}$-affine combinations
generate $\mathbb{Q}^{d}$. A fundamental invariant of $P$ is the number of
lattice points that it contains, $|P\cap\mathbb{Z}^{d}|$. More generally, let
$L_{P}(m)\coloneq|mP\cap\mathbb{Z}^{d}|$ count the lattice points in the $m$th
dilation $mP$ of $P$, where $m\in\mathbb{Z}_{\geq 0}$. Then $L_{P}$ is given
by a polynomial of degree $d$ called the  _Ehrhart polynomial_ [15]. The
corresponding generating series, called the  _Ehrhart series_ and denoted by
$\operatorname{Ehr}_{P}$, can be expressed as a rational function with
numerator a polynomial of degree at most $d$ [29]:

$\displaystyle\operatorname{Ehr}_{P}(t)\coloneq\sum_{m\geq
0}L_{P}(m)t^{m}=\frac{\delta_{0}+\delta_{1}t+\cdots+\delta_{d}t^{d}}{(1-t)^{d+1}},$
$\displaystyle\delta_{i}\in\mathbb{Z}.$

The coefficients $(\delta_{0},\delta_{1},\ldots,\delta_{d})$ of this
numerator, called the  _$\delta$ -vector_ or  _$h^{*}$ -vector_ of $P$, have
combinatorial meaning [14]:

  1. (i)

$\delta_{0}=1$;

  2. (ii)

$\delta_{1}=|P\cap\mathbb{Z}^{d}|-d-1$;

  3. (iii)

$\delta_{d}=|P^{\circ}\cap\mathbb{Z}^{d}|$, where
$P^{\circ}=P\setminus\partial P$ is the strict interior of $P$;

  4. (iv)

$\delta_{0}+\cdots+\delta_{d}=\operatorname{Vol}(P)$, where
$\operatorname{Vol}(P)=d!\operatorname{vol}(P)$ is the lattice-normalised
volume of $P$.

The polynomial $L_{P}$ can be expressed in terms of the $\delta$-vector:

$L_{P}(m)=\sum_{i=0}^{d}\delta_{i}\binom{d+m-i}{d}.$

From this we can see that the leading coefficient of $L_{P}$ is
$\operatorname{vol}(P)$, the Euclidean volume of $P$.

Given $d+1$ terms of the Ehrhart series, one can recover the $\delta$-vector
and hence the invariants $|P\cap\mathbb{Z}^{d}|$,
$|P^{\circ}\cap\mathbb{Z}^{d}|$, and $\operatorname{Vol}(P)$. This, however,
assumes knowledge of the dimension $d$.

######  Question 1.

Given a lattice polytope $P$, can machine learning recover the dimension $d$
of $P$ from sufficiently many terms of the Ehrhart series
$\operatorname{Ehr}_{P}$?

There has been recent success using machine learning (ML) to predict
invariants such as $\operatorname{Vol}(P)$ directly from the vertices of $P$
[4], and to predict numerical invariants from a geometric analogue of the
Ehrhart series called the Hilbert series [5]. As we will see in §§2.1–2.3, ML
is also extremely effective at answering Question 1. In §2.4 we propose a
possible explanation for this.

###  1.1. The quasi-period of a rational convex polytope

Now let $P$ be a convex polytope with  _rational_ vertices, and let
$k\in\mathbb{Z}_{>0}$ be the smallest positive dilation of $P$ such that $kP$
is a lattice polytope. One can define the Ehrhart series of $P$ exactly as
before:

$L_{P}(m)\coloneq|mP\cap\mathbb{Z}^{d}|.$

In general $L_{P}$ will no longer be a polynomial, but it is a  _quasi-
polynomial_ of degree $d$ and period $k$ [15, 27]. That is, there exist
polynomials $f_{0},\dots,f_{k-1}$, each of degree $d$, such that:

$\displaystyle L_{P}(qk+r)=f_{r}(q),$ whenever $0\leq r<k$.

The leading coefficient of each $f_{r}$ is $\operatorname{vol}(kP)$, the
Euclidean volume of the $k$th dilation of $P$.

It is sometimes possible to express $L_{P}$ using a smaller number of
polynomials. Let $\rho$ be the minimum number of polynomials needed to express
$L_{P}$ as a quasi-polynomial; $\rho$ is called the  _quasi-period_ of $P$,
and is a divisor of $k$. Each of these $\rho$ polynomials is of degree $d$
with leading coefficient $\operatorname{vol}(\rho P)$. When $\rho<k$ we say
that $P$ exhibits  _quasi-period collapse_ [26, 8, 16, 23, 7]. The Ehrhart
series of $P$ can be expressed as a rational function:

$\displaystyle\operatorname{Ehr}_{P}(t)=\frac{\delta_{0}+\delta_{1}t+\cdots+\delta_{\rho(d+1)-1}t^{\rho(d+1)-1}}{(1-t^{\rho})^{d+1}},$
$\displaystyle\delta_{i}\in\mathbb{Z}.$

As in the case of lattice polytopes, the $\delta$-vector carries combinatorial
information about $P$ [25, 21, 7, 6, 20]. Given $\rho(d+1)$ terms of the
Ehrhart series, one can recover the $\delta$-vector. This, however, requires
knowing both the dimension $d$ and the quasi-period $\rho$ of $P$.

######  Question 2.

Given a rational polytope $P$, can machine learning recover the dimension $d$
and quasi-period $\rho$ of $P$ from sufficiently many terms of the Ehrhart
series $\operatorname{Ehr}_{P}$?

ML again performs very well here – see §3. We propose a mathematical
explanation for this, in terms of forward differences and affine hyperplanes,
in §3.4.

###  1.2. Some motivating examples

######  Example 1.1.

The $2$-dimensional lattice polytope
$P\coloneq\operatorname{conv}\\{(-1,-1),(-1,2),(2,-1)\\}$ has volume
$\operatorname{Vol}(P)=9$, $|P\cap\mathbb{Z}^{2}|=10$, and
$|P^{\circ}\cap\mathbb{Z}^{2}|=1$. The Ehrhart polynomial of $P$ is:

$L_{P}(m)=\frac{9}{2}m^{2}+\frac{9}{2}m+1,$

and the Ehrhart series of $P$ is generated by:

$\operatorname{Ehr}_{P}(t)=\frac{1+7t+t^{2}}{(1-t)^{3}}=1+10t+28t^{2}+55t^{3}+91t^{4}+\cdots.$

######  Example 1.2.

The smallest dilation of the triangle
$P\coloneq\operatorname{conv}\\{(5,-1),(-1,-1),(-1,1/2)\\}$ giving a lattice
triangle is $2P$. Hence $L_{P}$ can be written as a quasi-polynomial of degree
$2$ and period $2$:

$f_{0}(q)=18q^{2}+9q+1,\qquad f_{1}(q)=18q^{2}+27q+10.$

This is not, however, the minimum possible period. As in Example 1.1,

$L_{P}(m)=\frac{9}{2}m^{2}+\frac{9}{2}m+1$

and thus $P$ has quasi-period $\rho=1$.

This striking example of quasi-period collapse is developed further in Example
4.1.

###  1.3. Code and data availability

The datasets used in this work were generated with V2.25-4 of Magma [9]. We
performed our ML analysis with scikit-learn [28], a standard machine learning
library for Python, using scikit-learn v0.24.1 and Python v3.8.8. All data,
along with the code used to generate it and to perform the subsequent
analysis, is available from Zenodo [11, 10] under a permissive open source
license (MIT for the code, CC0 for the data).

##  2\. Question 1: Dimension

In this section we investigate whether machine learning can predict the
dimension $d$ of a lattice polytope $P$ from sufficiently many terms of the
Ehrhart series $\operatorname{Ehr}_{P}$. We calculate terms $L_{P}(m)$ of the
Ehrhart series, for $0\leq m\leq 1100$, and encode these in a  _logarithmic
Ehrhart vector_ :

$(\log y_{0},\log y_{1},\ldots,\log y_{1100}),\qquad\text{ where
}y_{m}\coloneq L_{P}(m).$

We find that standard ML techniques are extremely effective, recovering the
dimension of $P$ with almost $100\%$ accuracy from the logarithmic Ehrhart
vector. We then ask whether ML can recover $\operatorname{Vol}(P)$ from this
Ehrhart data; again, this is achieved with near $100\%$ success.

###  2.1. Data generation

A dataset [11] containing $2918$ distinct entries, with $2\leq d\leq 8$, was
generated using Algorithm 2.1 below. The distribution of this data is
summarised in Table 1.

######  Algorithm 2.1.

x

Input:

A positive integer $d$.

Output:

A vector

$\big{(}\log y_{0},\log y_{1},\ldots,\log
y_{1100},d,\operatorname{Vol}(P)\big{)}$

for a $d$-dimensional lattice polytope $P$, where $y_{m}\coloneq L_{P}(m)$.

  1. (i)

Choose $d+k$ lattice points $\\{v_{1},\ldots,v_{d+k}\\}$ uniformly at random
in a box $[-5,5]^{d}$, where $k$ is chosen uniformly at random in
$\\{1,\ldots,5\\}$.

  2. (ii)

Set $P\coloneq\operatorname{conv}\\{v_{1},\ldots,v_{d+k}\\}$. If $\dim(P)\neq
d$ return to step (i).

  3. (iii)

Calculate the coefficients $y_{m}\coloneq L_{P}(m)$ of the Ehrhart series of
$P$, for $0\leq m\leq 1100$.

  4. (iv)

Return the vector $\big{(}\log y_{0},\log y_{1},\ldots,\log
y_{1100},d,\operatorname{Vol}(P)\big{)}$.

We deduplicated on the vector $\big{(}\log y_{0},\log y_{1},\ldots,\log
y_{1100},d,\operatorname{Vol}(P)\big{)}$ to get a dataset with distinct
entries. In particular, two polytopes that are equivalent under the group of
affine-linear transformations
$\operatorname{GL}_{d}(\mathbb{Z})\ltimes\mathbb{Z}^{d}$ give rise to the same
point in the dataset.

###  2.2. Machine learning the dimension

We reduced the dimensionality of the dataset by projecting onto the first two
principal components of the logarithmic Ehrhart vector. As one would expect
from Figure 1, a linear support vector machine (SVM) classifier trained on
these features predicted the dimension of $P$ with $100\%$ accuracy. Here we
used a scikit-learn pipeline consisting of a StandardScaler followed by an SVC
classifier with linear kernel and regularization hyperparameter $C=0.1$, using
$50\%$ of the data for training the classifier and tuning the hyperparameter,
and holding out the remaining $50\%$ of the data for model validation.

In §2.4 below we give a mathematical explanation for the structure observed in
Figure 1, and hence for why ML is so effective at predicting the dimension of
$P$. Note that the discussion in §3.4 suggests that one should also be able to
extract the dimension using ML on the  _Ehrhart vector_

$(y_{0},y_{1},\ldots,y_{1100})$

rather than the logarithmic Ehrhart vector

$(\log y_{0},\log y_{1},\ldots,\log y_{1100}).$

This is indeed the case, although here it is important not to reduce the
dimensionality of the data too much111This is consistent with the discussion
in §3.4, which suggests that we should try to detect whether the Ehrhart
vector lies in a union of linear subspaces that have fairly high codimension..
A linear SVM classifier trained on the full Ehrhart vector predicts the
dimension of $P$ with $98.7\%$ accuracy, with a pipeline exactly as above
except that $C=50000$. A linear SVM classifier trained on the first $30$
principal components of the Ehrhart vector (the same pipeline, but with
$C=20$) gives $93.7\%$ accuracy, but projecting to the first two components
reduces accuracy to $43.3\%$.

Dimension | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---
Total | 431 | 787 | 812 | 399 | 181 | 195 | 113
Table 1. The distribution of the dimensions appearing in the dataset for
Question 1. Figure 1. Projection onto the first two principal components of
the logarithmic Ehrhart vector coloured by the dimension of $P$.

###  2.3. Machine learning the volume

To learn the normalised volume of a lattice polytope from its logarithmic
Ehrhart vector, we used a scikit-learn pipeline consisting of a StandardScaler
followed by an SVR regressor with linear kernel and regularization
hyperparameter $C=340$. We restricted attention to the roughly $75\%$ of the
data with volume less than $10\,000$, thereby removing outliers. We used
$50\%$ of that data for training and hyperparameter tuning, selecting the
training set using a shuffle stratified by volume; this corrects for the fact
that the dataset contains a high proportion of polytopes with small volume.
The regression had a coefficient of determination ($R^{2}$) of $0.432$, and
gave a strong hint (see Figure 2) that we should repeat the analysis replacing
the logarithmic Ehrhart vector with the Ehrhart vector.

Using the same pipeline (but with $C=1000$) and the Ehrhart vector gives a
regression with $R^{2}=1.000$; see Figure 2. This regressor performs well over
the full dataset, with volumes ranging up to approximately $4.5$ million: over
the full dataset we still find $R^{2}=1.000$. The fact that Support Vector
Machine methods are so successful in recovering the volume of $P$ from the
Ehrhart vector is consistent with the discussion in §3.4.

Figure 2. Plotting the true versus the predicted normalised volume.

###  2.4. Crude asymptotics for $\log y_{k}$

Since

$L_{P}(m)=\operatorname{vol}(P)m^{d}+\text{lower order terms in~{}$m$}$

we have that

$\log y_{m}\sim d\log m+\log\operatorname{vol}(P).$

For $m\gg 0$, therefore, we see that the different components $\log y_{m}$ of
the logarithmic Ehrhart vector depend approximately affine-linearly on each
other as $P$ varies. It seems intuitively plausible that the first two PCA
components of the logarithmic Ehrhart vector should depend non-trivially on
$\log y_{m}$ for $m\gg 0$, and in fact this is the case – see Figure 3. Thus
the first two PCA components of the logarithmic Ehrhart vector should vary
approximately affine-linearly as $P$ varies, with constant slope and with a
translation that depends only on the dimension of $P$.

Figure 3. Contribution of $\log y_{m}$ to the first two principal components
of the logarithmic Ehrhart vector, as $m$ varies.

##  3\. Question 2: Quasi-period

Here we investigate whether ML can predict the quasi-period $\rho$ of a
$d$-dimensional rational polytope $P$ from sufficiently many terms of the
Ehrhart vector of $P$. Once again, standard ML techniques based on Support
Vector Machines are highly effective, achieving classification accuracies of
up to $95.3\%$. We propose a potential mathematical explanation for this in
§3.4.

###  3.1. Data generation

In this section we consider a dataset [10] containing $84\,000$ distinct
entries. This was obtained by using Algorithm 3.1 below to generate a larger
dataset, followed by random downsampling to a subset with 2000 datapoints for
each pair $(d,\rho)$ with $d\in\\{2,3,4\\}$ and $\rho\in\\{2,3,\ldots,15\\}$.

######  Algorithm 3.1.

x

Input:

A positive integer $d$.

Output:

A vector

$\big{(}\log y_{0},\log y_{1},\ldots,\log y_{1100},d,\rho\big{)}$

for a $d$-dimensional rational polytope $P$ with quasi-period $\rho$, where
$y_{m}\coloneq L_{P}(m)$.

  1. (i)

Choose $r\in\\{2,3,\ldots,15\\}$ uniformly at random.

  2. (ii)

Choose $d+k$ lattice points $\\{v_{1},\ldots,v_{d+k}\\}$ uniformly at random
in a box $[-5r,5r]^{d}$, where $k$ is chosen uniformly at random in
$\\{1,\ldots,5\\}$.

  3. (iii)

Set $P\coloneq\operatorname{conv}\\{v_{1},\ldots,v_{d+k}\\}$. If $\dim(P)\neq
d$ return to step (ii).

  4. (iv)

Choose a lattice point $v\in P\cap\mathbb{Z}^{d}$ uniformly at random and
replace $P$ with the translation $P-v$. (We perform this step to ensure that
the resulting rational polytope always contains a lattice point; this avoids
complications when taking $\log$ in step (viii).)

  5. (v)

Replace $P$ with the dilation $P/r$.

  6. (vi)

Calculate the coefficients $y_{m}\coloneq L_{P}(m)$ of the Ehrhart series of
$P$, for $0\leq m\leq 1100$.

  7. (vii)

Calculate the quasi-period $\rho$.

  8. (viii)

Return the vector $\big{(}\log y_{0},\log y_{1},\ldots,\log
y_{1100},d,\rho\big{)}$.

As before, we deduplicated the dataset on the vector $\big{(}\log y_{0},\log
y_{1},\ldots,\log y_{1100},d,\rho\big{)}$.

###  3.2. Recovering the dimension and volume

Figure 4 shows the first two principal components of the logarithmic Ehrhart
vector. As in §2.2, this falls into widely-separated linear clusters according
to the value of $\dim(P)$, and so the dimension of the rational polytope $P$
can be recovered with high accuracy from its logarithmic Ehrhart vector.
Furthermore, as in §2.3, applying a linear SVR regressor (with $C=1000$) to
the Ehrhart vector predicts the volume of a rational polytope $P$ with high
accuracy ($R^{2}=1.000$).

Figure 4. Projection onto the first two principal components of the
logarithmic Ehrhart vector data coloured by the dimension of $P$.

###  3.3. Machine learning the quasi-period

To learn the quasi-period of a rational polytope from its Ehrhart vector, we
used a scikit-learn pipeline consisting of a StandardScaler followed by a
LinearSVC classifier. We fixed a dimension $d\in\\{2,3,4\\}$, moved to PCA co-
ordinates, and used $50\%$ of the data ($N=14000$) for training the classifier
and hyperparameter tuning, holding out the remaining $50\%$ of the data for
model validation. Results are summarised on the left-hand side of Table 2,
with learning curves in the left-hand column of Figure 5 and confusion
matrices in the left-hand column of Figure 6. The confusion matrices hint at
some structure in the misclassified data.

One could also use the same pipeline but applied to the logarithmic Ehrhart
vector rather than the Ehrhart vector. Results are summarised on the right-
hand side of Table 2, with learning curves in the right-hand column of Figure
5 and confusion matrices in the right-hand column of Figure 6. Using the
logarithmic Ehrhart data resulted in a less accurate classifier, but the
learning curves suggest that this might be improved by adding more training
data. Again there are hints of structure in the misclassified data.

Ehrhart vector |  | log Ehrhart vector
---|---|---
Dimension | $C$ | Accuracy |  | Dimension | $C$ | Accuracy
2 | 0.01 | 80.6% |  | 2 | 1 | 79.7%
3 | 0.001 | 94.0% |  | 3 | 1 | 85.2%
4 | 0.001 | 95.3% |  | 4 | 1 | 83.2%
Table 2. The regularization hyperparameter $C$ and accuracy for a LinearSVC
classifier predicting the quasi-period from the Ehrhart vector and logarithmic
Ehrhart vector of rational polytopes $P$.

Figure 5. Quasi-period learning curves for a linear SVM classifier, for both
the Ehrhart vector (left-hand column) and the logarithmic Ehrhart vector
(right-hand column).

Figure 6. Row-normalised confusion matrices for a linear SVM classifier for
the quasi-period, for both the Ehrhart vector (left-hand column) and the
logarithmic Ehrhart vector (right-hand column).

###  3.4. Forward differences

Let $y$ denote the sequence $(y_{m})_{m=0}^{\infty}$. Recall the forward
difference operator $\Delta$ defined on the space of sequences:

$\Delta y=(y_{m+1}-y_{m})_{m=0}^{\infty}.$

A sequence $y$ depends polynomially on $m$, that is

$\displaystyle y_{m}=a_{0}+a_{1}m+\cdots+a_{d}m^{d}$ for some
$a_{0},\ldots,a_{d}\in\mathbb{R}$,

if and only if $y$ lies in the kernel of $\Delta^{d+1}$. Furthermore, in this
case, $\Delta^{d}y$ is the constant sequence with value $d!a_{d}$.

Thus a sequence $(y_{m})_{m=0}^{\infty}$ is quasi-polynomial of degree $d$ and
period $k$, in the sense of §1.1, if and only if it lies in the kernel of
$\Delta_{k}^{d+1}$ where $\Delta_{k}$ is the $k$-step forward difference
operator,

$\Delta_{k}\left((y_{m})_{m=0}^{\infty}\right)=(y_{m+k}-y_{m})_{m=0}^{\infty},$

and $(y_{m})_{m=0}^{\infty}$ does not lie in the kernel of $\Delta_{k}^{d}$.
Furthermore, in this case, we can determine the leading coefficients of the
polynomials $f_{0},\ldots,f_{k-1}$ by examining the values of the $k$-periodic
sequence $\Delta_{k}^{d}y$. When $y$ arises as the Ehrhart series of a
rational polytope $P$, all of these constant terms equal the volume of $kP$,
and so the value of the constant sequence $\Delta_{k}^{d}y$ determines the
normalised volume $\operatorname{Vol}(P)$.

This discussion suggests that an SVM classifier with linear kernel should be
able to learn the quasi-period and volume with high accuracy from the Ehrhart
vector of a rational polytope, at least if we consider only polytopes of a
fixed dimension $d$. Having quasi-period $k$ amounts to the Ehrhart vector
lying in (a relatively open subset of) a certain subspace
$\ker\Delta_{k}^{d+1}$; these subspaces, being linear objects, should be
easily separable using hyperplanes. Similarly, having fixed normalised volume
amounts to lying in a given affine subspace; such affine subspaces should be
easily separable using affine hyperplanes.

From this point of view it is interesting that an SVM classifier with linear
kernel also learns the quasi-period with reasonably high accuracy from the
logarithmic Ehrhart vector. Passing from the Ehrhart vector to the logarithmic
Ehrhart vector replaces the linear subspaces $\ker\Delta_{k}^{d+1}$ by non-
linear submanifolds. But our experiments above suggest that these non-linear
submanifolds must nonetheless be close to being separable by appropriate
collections of affine hyperplanes.

##  4\. A remark on the Gorenstein index

In this section we discuss a geometric question where machine learning
techniques failed. This involves a more subtle combinatorial invariant called
the Gorenstein index, and was the question that motivated the rest of the work
in this paper. We then suggest why in retrospect we should not have expected
to be able to answer this question using machine learning (or at all).

###  4.1. The Gorenstein index

Fix a rank $d$ lattice $N\cong\mathbb{Z}^{d}$, write $N_{\mathbb{Q}}\coloneq
N\otimes_{\mathbb{Z}}\mathbb{Q}$, and let $P\subset N_{\mathbb{Q}}$ be a
lattice polytope. The  _polar polyhedron_ of $P$ is given by

$P^{*}\coloneq\\{u\in M\otimes_{\mathbb{Z}}\mathbb{Q}\mid u(v)\geq-1\text{ for
all }v\in P\\}$

where $M\coloneq\operatorname{Hom}(N,\mathbb{Z})\cong\mathbb{Z}^{d}$ is the
lattice dual to $N$. The polar polyhedron $P^{*}$ is a convex polytope with
rational vertices if and only if the origin lies in the strict interior of
$P$. In this case, $(P^{*})^{*}=P$. The smallest positive integer $k_{P}$ such
that $k_{P}P^{*}$ is a lattice polytope is called the  _Gorenstein index_ of
$P$.

The Gorenstein index arises naturally in the context of  _Fano toric
varieties_ , and we will restrict our discussion to this setting. Let
$P\subset N_{\mathbb{Q}}$ be a lattice polytope such that the vertices of $P$
are primitive lattice vectors and that the origin lies in the strict interior
of $P$; such polytopes are called  _Fano_ [22]. The spanning fan $\Sigma_{P}$
of $P$ – that is, the complete fan whose cones are generated by the faces of
$P$ – gives rise to a Fano toric variety $X_{P}$ [17, 12]. This construction
gives a one-to-one correspondence between
$\operatorname{GL}_{d}(\mathbb{Z})$-equivalence classes of Fano polytopes and
isomorphism classes of Fano toric varieties. Let $P\subset N_{\mathbb{Q}}$ be
a Fano polytope that corresponds to a Fano toric variety $X\coloneq X_{P}$.
The polar polytope $P^{*}\subset M_{\mathbb{Q}}$ then corresponds to a divisor
on $X$ called the  _anticanonical divisor_ , which is denoted by $-K_{X}$. In
general $-K_{X}$ is an ample $\mathbb{Q}$-Cartier divisor, and $X$ is
$\mathbb{Q}$-Gorenstein. The Gorenstein index $k_{P}$ of $P$ is equal to the
smallest positive multiple $k$ of the anticanonical divisor such that
$-kK_{X}$ is Cartier.

Under the correspondence just discussed, the Ehrhart series
$\operatorname{Ehr}_{P^{*}}$ of the polar polytope $P^{*}$ coincides with the
Hilbert series $\operatorname{Hilb}_{X}(-K_{X})$. The Hilbert series is an
important numerical invariant of $X$, and it makes sense to ask whether the
Gorenstein index of $X$ is determined by the Hilbert series. Put differently:

######  Question 3.

Given a Fano polytope $P$, can ML recover the Gorenstein index $k_{P}$ of $P$
from sufficiently many terms of the Ehrhart series
$\operatorname{Ehr}_{P^{*}}$ of the polar polytope $P^{*}$?

There are good reasons, as we discuss below, to expect the answer to Question
3 to be ‘no’. But part of the power of ML in mathematics is that it can detect
or suggest structure that was not known or expected previously (see e.g.
[13]). That did not happen on this occasion: applying the techniques discussed
in §2 and §3 did not allow us to predict the Gorenstein index of $P$ from the
Ehrhart series of the polar polytope $\operatorname{Ehr}_{P^{*}}$.

###  4.2. Should we have expected this?

The Hilbert series is preserved under an important class of deformations
called  _qG-deformations_ [24]. But the process of  _mutation_ [2] can
transform a Fano polytope $P$ to a Fano polytope $Q$ with
$\operatorname{Ehr}_{P^{*}}=\operatorname{Ehr}_{Q^{*}}$. Mutation gives rise
to a qG-deformation from $X_{P}$ to $X_{Q}$ [1], but need not preserve the
Gorenstein index: $k_{P}$ need not be equal to $k_{Q}$. Thus the Gorenstein
index is not invariant under qG-deformation. It might have been unrealistic to
expect that a qG-deformation invariant quantity (the Hilbert series) could
determine an invariant (the Gorenstein index) which can vary under qG-
deformation.

### Quasiperiod collapse

Although the phenomenon of quasi-period collapse remains largely mysterious
from a combinatorial view-point, in the context of toric geometry one possible
explanation arises from mutation and qG-deformation [23]. The following
example revisits Examples 1.1 and 1.2 from this point of view, and illustrates
why Question 3 cannot have a meaningful positive answer.

$\longrightarrow$

Figure 7. The mutation from $P_{(1,1,1)}$ to $P_{(1,1,2)}$ gives a scissors
congruence between the polar polytopes.

######  Example 4.1.

Let $P_{(a,b,c)}\subset N_{\mathbb{Q}}$ denote the $2$-dimensional Fano
polytope associated with weighted projective space
$\mathbb{P}(a^{2},b^{2},c^{2})$, where $a$, $b$, $c$ are pairwise coprime
positive integers. Then
$P_{(1,1,1)}\coloneq\operatorname{conv}\\{(1,0),(0,1),(-1,-1)\\}$ is the Fano
polygon associated with $\mathbb{P}^{2}$, with polar polygon
$P^{*}_{(1,1,1)}\subset M_{\mathbb{Q}}$ the lattice triangle appearing in
Example 1.1.

The graph of mutations of $P_{(1,1,1)}$ has been completely described [19, 3].
Up to $\operatorname{GL}_{2}(\mathbb{Z})$-equivalence, there is exactly one
mutation from $P_{(1,1,1)}$: this gives $P_{(1,1,2)}$ corresponding to
$\mathbb{P}(1,1,4)$. The polar polygon $P^{*}_{(1,1,2)}$ is the rational
triangle in Example 1.2. As Figure 7 illustrates, a mutation between polytopes
gives rise to a scissors congruence [18] between polar polytopes. Mutation
therefore preserves the Ehrhart series of the polar polytope, and this
explains why we have quasi-period collapse in this example.

We can mutate $P_{(1,1,2)}$ in two ways that are distinct up to the action of
$\operatorname{GL}_{2}(\mathbb{Z})$: one returns us to $P_{(1,1,1)}$, whilst
the other gives $P_{(1,2,5)}$. Continuing to mutate, we obtain an infinite
graph of triangles $P_{(a,b,c)}$, where the $(a,b,c)$ are the  _Markov
triples_ , that is, the positive integral solutions to the Markov equation:

$3xyz=x^{2}+y^{2}+z^{2}.$

The Gorenstein index of $P_{(a,b,c)}$ is $abc$. In particular, the Gorenstein
index can be made arbitrarily large whilst the Ehrhart series, and hence
quasi-period, of $P^{*}_{(a,b,c)}$ is fixed. See [23] for details.

##  5\. Conclusion

We have seen that Support Vector Machine methods are very effective at
extracting the dimension and volume of a lattice or rational polytope $P$, and
the quasi-period of a rational polytope $P$, from the initial terms of its
Ehrhart series. We have also seen that ML methods are unable to reliably
determine the Gorenstein index of a Fano polytope $P$ from the Ehrhart series
of its polar polytope $P^{*}$. The discussions in §2.4, §3.4, and §4.2 suggest
that these results are as expected: that ML is detecting known and understood
structure in the dimension, volume, and quasiperiod cases, and that there is
probably no structure to detect in the Gorenstein index case. But there is a
more useful higher-level conclusion to draw here too: when applying ML methods
to questions in pure mathematics, one needs to think carefully about methods
and results. Questions 2 and 3 are superficially similar, yet one is amenable
to ML and the other is not. Furthermore, applying standard ML recipes in a
naive way would have led to false negative results. For example, since the
Ehrhart series grows so fast, it would have been typical to suppress the
growth rate by taking logarithms, and also to pass to principal components.
Taking logarithms is a good idea for some of our questions but not for others;
this reflects the mathematical realities underneath the data, and not just
whether the vector-components $y_{m}$ involved grow rapidly with $m$ or not.
Passing to principal components is certainly a useful tool, but naive feature
extraction would have retained only the first principal component, which is
responsible for more than $99.999\%$ of the variation (in both the logarithmic
Ehrhart vector and the Ehrhart vector). This would have left us unable to
detect the positive answers to Questions 1 and 2, in the former case because
projection to one dimension amalgamates clusters (see Figure 1) and in the
latter case because we need to detect whether the Ehrhart vector lies in a
certain high-codimension linear subspace and that structure is destroyed by
projection to a low-dimensional space.

### Acknowledgments

TC is supported by ERC Consolidator Grant 682603 and EPSRC Programme Grant
EP/N03189X/1. JH is supported by a Nottingham Research Fellowship. AK is
supported by EPSRC Fellowship EP/N022513/1.

## References

  * [1] Mohammad Akhtar, Tom Coates, Alessio Corti, Liana Heuberger, Alexander M. Kasprzyk, Alessandro Oneto, Andrea Petracci, Thomas Prince, and Ketil Tveiten, _Mirror symmetry and the classification of orbifold del Pezzo surfaces_ , Proc. Amer. Math. Soc. 144 (2016), no. 2, 513–527.
  * [2] Mohammad Akhtar, Tom Coates, Sergey Galkin, and Alexander M. Kasprzyk, _Minkowski polynomials and mutations_ , SIGMA Symmetry Integrability Geom. Methods Appl. 8 (2012), Paper 094, 17.
  * [3] Mohammad E. Akhtar and Alexander M. Kasprzyk, _Mutations of fake weighted projective planes_ , Proc. Edinb. Math. Soc. (2) 59 (2016), no. 2, 271–285.
  * [4] Jiakang Bao, Yang-Hui He, Edward Hirst, Johannes Hofscheier, Alexander M. Kasprzyk, and Suvajit Majumder, _Polytopes and Machine Learning_ , arXiv:2109.09602 [math.CO], 2021.
  * [5] by same author, _Hilbert series, machine learning, and applications to physics_ , Phys. Lett. B 827 (2022), Paper No. 136966, 8.
  * [6] Matthias Beck, Benjamin Braun, and Andrés R. Vindas-Meléndez, _Decompositions of Ehrhart $h^{*}$-Polynomials for Rational Polytopes_, Discrete Comput. Geom. 68 (2022), no. 1, 50–71.
  * [7] Matthias Beck, Sophia Elia, and Sophie Rehberg, _Rational Ehrhart theory_ , arXiv:2110.10204 [math.CO], 2022\.
  * [8] Matthias Beck, Steven V. Sam, and Kevin M. Woods, _Maximal periods of (Ehrhart) quasi-polynomials_ , J. Combin. Theory Ser. A 115 (2008), no. 3, 517–525.
  * [9] Wieb Bosma, John Cannon, and Catherine Playoust, _The Magma algebra system. I. The user language_ , J. Symbolic Comput. 24 (1997), no. 3-4, 235–265, Computational algebra and number theory (London, 1993).
  * [10] Tom Coates, Johannes Hofscheier, and Alexander M. Kasprzyk, _Ehrhart series coefficients and quasi-period for random rational polytopes_ , Zenodo, doi:10.5281/zenodo.6614829, 2022\.
  * [11] by same author, _Ehrhart series coefficients for random lattice polytopes_ , Zenodo, doi:10.5281/zenodo.6614821, 2022\.
  * [12] David A. Cox, John B. Little, and Henry K. Schenck, _Toric varieties_ , Graduate Studies in Mathematics, vol. 124, American Mathematical Society, Providence, RI, 2011.
  * [13] Alex Davies, Petar Veličković, Lars Buesing, Sam Blackwell, Daniel Zheng, Nenad Tomašev, Richard Tanburn, Peter Battaglia, Charles Blundell, András Juhász, Marc Lackenby, Geordie Williamson, Demis Hassabis, and Pushmeet Kohli, _Advancing mathematics by guiding human intuition with AI_ , Nature 600 (2021), 70–74.
  * [14] E. Ehrhart, _Sur un problème de géométrie diophantienne linéaire. II. Systèmes diophantiens linéaires_ , J. Reine Angew. Math. 227 (1967), 25–49.
  * [15] Eugène Ehrhart, _Sur les polyèdres homothétiques bordés à $n$ dimensions_, C. R. Acad. Sci. Paris 254 (1962), 988–990.
  * [16] Matthew H. J. Fiset and Alexander M. Kasprzyk, _A note on palindromic $\delta$-vectors for certain rational polytopes_, Electron. J. Combin. 15 (2008), no. 1, Note 18, 4.
  * [17] William Fulton, _Introduction to toric varieties_ , Annals of Mathematics Studies, vol. 131, Princeton University Press, Princeton, NJ, 1993, The William H. Roever Lectures in Geometry.
  * [18] Christian Haase and Tyrrell B. McAllister, _Quasi-period collapse and ${\rm GL}_{n}(\mathbb{Z})$-scissors congruence in rational polytopes_, Integer points in polyhedra—geometry, number theory, representation theory, algebra, optimization, statistics, Contemp. Math., vol. 452, Amer. Math. Soc., Providence, RI, 2008, pp. 115–122.
  * [19] Paul Hacking and Yuri Prokhorov, _Smoothable del Pezzo surfaces with quotient singularities_ , Compos. Math. 146 (2010), no. 1, 169–192.
  * [20] Girtrude Hamm, Johannes Hofscheier, and Alexander M. Kasprzyk, _Half-integral polygons with a fixed number of lattice points_ , in preparation.
  * [21] Andrew John Herrmann, _Classification of Ehrhart quasi-polynomials of half-integral polygons_ , Master’s thesis, San Francisco State University, 2010\.
  * [22] Alexander M. Kasprzyk and Benjamin Nill, _Fano polytopes_ , Strings, gauge fields, and the geometry behind, World Sci. Publ., Hackensack, NJ, 2013, pp. 349–364.
  * [23] Alexander M. Kasprzyk and Ben Wormleighton, _Quasi-period collapse for duals to Fano polygons: an explanation arising from algebraic geometry_ , arXiv:1810.12472 [math.CO], 2018.
  * [24] J. Kollár and N. I. Shepherd-Barron, _Threefolds and deformations of surface singularities_ , Invent. Math. 91 (1988), no. 2, 299–338.
  * [25] Tyrrell B. McAllister, _Coefficient functions of the Ehrhart quasi-polynomials of rational polygons_ , Proceedings of the 2008 International Conference on Information Theory and Statistical Learning, ITSL 2008, Las Vegas, Nevada, USA, July 14-17, 2008 (Matthias Dehmer, Michael Drmota, and Frank Emmert-Streib, eds.), CSREA Press, 2008.
  * [26] Tyrrell B. McAllister and Kevin M. Woods, _The minimum period of the Ehrhart quasi-polynomial of a rational polytope_ , J. Combin. Theory Ser. A 109 (2005), no. 2, 345–352.
  * [27] P. McMullen, _Lattice invariant valuations on rational polytopes_ , Arch. Math. (Basel) 31 (1978/79), no. 5, 509–516.
  * [28] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, _Scikit-learn: Machine learning in Python_ , Journal of Machine Learning Research 12 (2011), 2825–2830.
  * [29] Richard P. Stanley, _Decompositions of rational convex polytopes_ , Ann. Discrete Math. 6 (1980), 333–342.