marginparsep has been altered.
topmargin has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Stopping Bayesian Optimization with Probabilistic Regret Bounds

James T. Wilson

††footnotetext: Correspondence to<EMAIL_ADDRESS>
Preliminary work. Under review by the International Conference on Machine
Learning (ICML).

###### Abstract

Bayesian optimization is a popular framework for efficiently finding high-
quality solutions to difficult problems based on limited prior information. As
a rule, these algorithms operate by iteratively choosing what to try next
until some predefined budget has been exhausted. We investigate replacing this
de facto stopping rule with an $(\epsilon,\delta)$-criterion: _stop when a
solution has been found whose value is within $\epsilon>0$ of the optimum with
probability at least $1-\delta$ under the model._ Given access to the prior
distribution of problems, we show how to verify this condition in practice
using a limited number of draws from the posterior. For Gaussian process
priors, we prove that Bayesian optimization with the proposed criterion stops
in finite time and returns a point that satisfies the
$(\epsilon,\delta)$-criterion under mild assumptions. These findings are
accompanied by extensive empirical results which demonstrate the strengths and
weaknesses of this approach.

##  1 Introduction

In the real world, we are often interested in finding a nearly optimal
solution for a given problem. More specifically, we would like to find a
solution that is sufficient for the task at hand Simon (1956). Many of these
problems are not only expensive to solve but difficult to reason about without
extensive background knowledge—such as discovering new chemicals Griffiths &
Hernández-Lobato (2020), designing better experiments von Kügelgen et al.
(2019), or configuring machine learning algorithms Snoek et al. (2012).

A common approach is therefore to construct models for these problem and use
them to predict real-world outcomes. In recent years, Bayesian optimization
(BO) has emerged as a leading approach for accomplishing these tasks. Precise
definitions vary, but BO methods are frequently characterized by their use of
probabilistic models to guide the search for good solutions. The idea is for
these models to provide predictive distributions over unknown results, which
can then be used to simulate the usefulness of trying different things. These
methods enjoy continued success on an increasingly diverse set of problems,
many of which are beyond the scope of our presentation here. For a recent
review of these topics, see Garnett (2023).

Despite these successes, an ongoing issue for practitioners has been the
persistent lack of user-friendly stopping rules. The vast majority of BO runs
proceed until a predetermined budget (usually, a number of trials or amount of
resources) has been exhausted. While there are many reasons for this trend,
the primary cause is likely mistrust in the underlying model. After all, if
the model is bad then the algorithm is liable to stop much too soon or far too
late. To avoid potential disappointment, let us say upfront that this
criticism is valid and will not be resolved in this work. We will return to
this point when discussing experiments in Section 6.

At the same time, there are many applications where models have repeatedly
proved to be reliable. Especially in these settings, there is much to be
gained by using models to help us decide when to stop. One benefit of model-
based stopping is its ability to adapt to the data. Sometimes we will get
lucky and stumble upon a good solution early on. Other times, we will struggle
to make progress. If the model accounts for these events, then stopping
decisions can be tailored to reflect the peculiarities of individual runs.

Another benefit of model-based stopping is that it can simplify the user
experience by allowing us to better communicate our goals. Choosing a budget
is fairly unintuitive. Unless one has dealt with many similar problems before,
it can be difficult to gauge the amount of resources to allocate for search.
One of our aims is therefore to express stopping conditions in a language that
is easier to understand. To this end, we focus on the case where the user
deems a solution sufficient if its value is within $\epsilon>0$ of the best
possible outcome with probability at least $1-\delta$. We show how this goal
can be accomplished in a principled and practical way.

The remaining text is organized as follows. Section 2 reviews background
material. Section 3 presents the idea behind our approach. Section 4 uses this
premise to construct a practical algorithm. Section 5 analyzes this
algorithm’s properties when the true model is known. Finally, Section 6
examines performance under ideal and realistic conditions.

Figure 1:  _Left:_ Posterior mean and two standard (blue) given 8 observations
(black dots) of $y(\cdot)\sim\mathcal{N}\left(f(\cdot),10^{-2}\right)$. The
goal is to find an $\boldsymbol{x}\in\mathcal{X}$ whose true function (dashed
black line) value is within $\epsilon>0$ of the optimum (orange star).
_Middle:_ Draws of $f_{t}\sim\mathcal{GP}(\mu_{t},k_{t})$ are shown in blue
alongside corresponding samples of $f^{*}_{t}$ (orange stars). _Right:_
Various approximations for (4) under the given posterior. Ground truth (dashed
black) was established by assessing (6) with location-scale-based samples of
$f_{t}$ on a dense grid. All other methods simulated $f_{t}$ using the
approximate sampling method from Wilson et al. (2020).  The resulting
approximation to (6) is portrayed in blue. The remaining methods analytically
integrated out $f_{t}(\boldsymbol{x})$, while potentially conditioning on
$f^{*}_{t}$. $\operatorname{\mathbb{P}}(r_{t}(\boldsymbol{x})\geq\epsilon)$ is
shown in orange,
$\operatorname{\mathbb{P}}(r_{t}(\boldsymbol{x})\geq\epsilon\mid
f_{t}(\boldsymbol{x})\leq f^{*}_{t})$ in green, and
$\operatorname{\mathbb{P}}(r_{t}(\boldsymbol{x})\geq\epsilon\mid
f_{t}(\boldsymbol{x})\leq f^{*}_{t},f_{t}(\boldsymbol{x}^{*})=f^{*}_{t})$ in
red where
$\boldsymbol{x}^{*}\in\operatorname*{\arg\min}_{\boldsymbol{x}\in\mathcal{X}}f_{t}(\boldsymbol{x})$
is sampled.

##  2 Background

We will use lowercase and uppercase boldface characters to denote vectors and
matrices, respectively. By minor abuse of notation, we will write, e.g.,
$\mathbf{X}_{n}=\left\\{\boldsymbol{x}_{i}:i=1,\ldots,n\right\\}$ for a set
consisting of $n\in\mathbb{N}$ items, typically the first $n$ elements of a
sequence $(\boldsymbol{x}_{n})$. For a function $h:\mathcal{X}\to\mathbb{R}$,
we will similarly write
$h(\mathbf{X}_{n})=\left\\{h(\boldsymbol{x}_{i}):i=1,\ldots,n\right\\}$ and
sometimes treat this set as an $n$-dimensional vector.

The left panel of Figure 1 portrays a prototypical problem, the details of
which are expanded upon below. We are given a function
$f:\mathcal{X}\to\mathbb{R}$ and tasked with finding a point
$\boldsymbol{x}\in\mathcal{X}$ whose value is within $\epsilon>0$ of the
optimum $f^{*}=\sup_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})$. Such a
point will be said to be _$\epsilon$ -optimal_.

At time $t\in\mathbb{N}$, our understanding of $f$ will be based on prior
knowledge and any data that we have collected so far. We will therefore use a
model to help us make informed decisions. Different types of models are
eligible and methods introduced below will only require the ability to
simulate whether points are $\epsilon$-optimal. We focus on the most popular
model family in this setting: Gaussian processes (GPs).

A Gaussian process is a random function
$f\\!:\\!\mathcal{X}\times\Omega\to\mathbb{R}$ with the property that, for any
finite subset $\mathbf{X}\subseteq\mathcal{X}$, the random variable
$f(\mathbf{X}):\Omega\to\mathbb{R}^{\left|\mathbf{X}\right|}$ is Gaussian in
distribution. We will write $f\sim\mathcal{GP}(0,k)$ for a centered GP with
covariance function $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$. We model
the observed outcomes of trials $\mathbf{X}_{t}\in\mathcal{X}^{t}$ as function
values corrupted by independent Gaussian noise and denote them by
$y(\mathbf{X}_{t})\sim\mathcal{N}\left(f(\mathbf{X}_{t}),\gamma^{2}\mathbf{I}\right)$.
Given $y(\mathbf{X}_{t})$, our understanding of $f$ is represented by the
posterior process $f_{t}\sim\mathcal{GP}(\mu_{t},k_{t})$, where
$\mathbf{\Lambda}=k(\mathbf{X}_{t},\mathbf{X}_{t})+\gamma^{2}\mathbf{I}$
defines

$\displaystyle\mu_{t}(\boldsymbol{x})$
$\displaystyle=k(\boldsymbol{x},\mathbf{X}_{t})\mathbf{\Lambda}^{-1}y(\mathbf{X}_{t})$
(1) $\displaystyle k_{t}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=k(\boldsymbol{x},\boldsymbol{x}^{\prime})-k(\boldsymbol{x},\mathbf{X}_{t})\mathbf{\Lambda}^{-1}k(\mathbf{X}_{t},\boldsymbol{x}^{\prime}).$
(2)

An example of a GP posterior can be seen in Figure 1.

Throughout, we assume that $\mathcal{X}$ is compact and that $\mu_{t}$ and
$k_{t}$ are both continuous so that they attain their limits on $\mathcal{X}$.
Among other things, this assumption allows us to write
$\boldsymbol{s}_{t}\in\operatorname*{\arg\max}_{\boldsymbol{x}\in\mathbf{S}_{t}}\mu_{t}(\boldsymbol{x})$
for a preferred solution at time $t$, where $\mathbf{S}_{t}$ is either the set
of queried points $\mathbf{X}_{t}$ or space $\mathcal{X}$.

Figure 2:  Simulated stopping behavior using $\operatorname{PRB}$ when
$f:[0,1]^{2}\to\mathbb{R}$ is drawn from a known prior with noise variance
$\gamma^{2}=10^{-4}$. Regret bounds $\epsilon>0$ control how close
$f(\boldsymbol{x})$ must be to the optimal value $f^{*}$ for the point
$\boldsymbol{x}\in\mathcal{X}$ to be deemed satisfactory. Risk tolerances
$\delta>0$ upper bound the probability that the returned point fails to be
$\epsilon$-optimal. _Left:_ Percent of runs that stopped before time $T=128$.
_Middle:_ Percent of stopped runs that returned $\epsilon$-optimal points.
_Right:_ Median number of trials performed by stopped runs.

##  3 Premise

In practice, we will rarely ever know if a point
$\boldsymbol{x}\in\mathcal{X}$ is truly $\epsilon$-optimal. Instead, we must
rely upon a model to provide an educated guess. If this model is
probabilistic, then we would like to know its view on the probability that
$\boldsymbol{x}$ is $\epsilon$-optimal. More formally, define the _regret_
incurred by a solution $\boldsymbol{x}\in\mathcal{X}$ under the model at time
$t$ as

$\displaystyle r_{t}(\boldsymbol{x})=f^{*}_{t}-f_{t}(\boldsymbol{x}),$ (3)

where $f^{*}_{t}=\sup_{\boldsymbol{x}\in\mathcal{X}}f_{t}(\boldsymbol{x})$ is
the best predicted outcome. The model-based probability of $\boldsymbol{x}$
being $\epsilon$-optima is then111We will often omit $\epsilon$ as an argument
of $\Psi_{t}$ to ease notation.

$\displaystyle\Psi_{t}(\boldsymbol{x};\epsilon)=\operatorname{\mathbb{P}}\left(r_{t}(\boldsymbol{x})\leq\epsilon\right).$
(4)

We propose to search for better solutions until we have found a point that
satisfies a _probabilistic regret bound_

$\displaystyle\operatorname{PRB}_{t}(\boldsymbol{x};\epsilon,\delta)=\begin{cases}1&\Psi_{t}(\boldsymbol{x};\epsilon)\geq
1-\delta\\\ 0&\mathrel{\text{missing}}{otherwise,}\end{cases}$ (5)

where _risk tolerance_ $\delta>0$ upper bounds the probability that we fail to
return an $\epsilon$-optimal solution. Points that meet this requirement will
be said to be _$(\epsilon,\delta)$ -optimal_.

As a rule, the $\epsilon$-optimality probability (4) is intractable. We
therefore assume that samples of $f^{*}_{t}$ and $f_{t}(\boldsymbol{x})$ will
be used to generate a Monte Carlo estimate

$\displaystyle\Psi_{t}^{n}(\boldsymbol{x};\epsilon)=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\left(r_{t}^{i}(\boldsymbol{x})\leq\epsilon\right),$
(6)

where $r_{t}^{i}(\boldsymbol{x})$ is the $i$-th draw of the model-based regret
(3). Figure 1 illustrates this procedure by visualizing draws of $f_{t}$ and
comparing different approximations to (4).

We now would like to use $\Psi_{t}^{n}$ to test (5). While
$\Psi_{t}^{n}(\boldsymbol{x})\geq 1-\delta$ does not guarantee that
$\Psi_{t}(\boldsymbol{x})\geq 1-\delta$, we may bound the estimation error
inherent to $\Psi_{t}^{n}(\boldsymbol{x})$ with high probability and use this
bound to infer whether the latter event is likely to hold. We accomplish this
task with the help of the following inequality from Audibert et al. (2007), a
variant of which was also studied by Maurer & Pontil (2009),

######  Theorem 1 (Empirical Bernstein bound).

Let $X_{1},\ldots,X_{n}$ be independent and identically distributed copies of
a random variable $X\in[a,b]$ with empirical mean and variance

$\displaystyle\overline{X}_{n}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$
$\displaystyle\overline{\sigma}_{n}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left(X_{i}-\overline{X}_{n}\right)^{2}.$
(7)

It follows with probability at least $1-\delta$ that

$\displaystyle\left|\overline{X}_{n}-\operatorname{\mathbb{E}}[X]\right|$
$\displaystyle\leq\Delta_{n},$ (8)

where
$\Delta_{n}=\overline{\sigma}_{n}\sqrt{\frac{2\log(3/\delta)}{n}}+\frac{3(b-a)\log(3/\delta)}{n}$.

This bound was chosen because it decays faster than, e.g., Hoeffding’s
inequality when the sample variance is significantly smaller than the sample
range. We are interested in inferring whether an event occurs with high
probability, so this variance will vanish as the true probability goes to one.

In order to test the criterion (5), choose any nonzero probabilities
$\delta_{\mathrm{mod}}$ and $\delta_{\mathrm{est}}$ so that
$(1-\delta_{\mathrm{mod}})(1-\delta_{\mathrm{est}})\geq 1-\delta$. As we will
see shortly, $\lambda=1-\delta_{\mathrm{mod}}$ lower bounds the chance that
$\boldsymbol{x}$ is $\epsilon$-optimal under the model
$\Psi_{t}(\boldsymbol{x})=\operatorname{\mathbb{E}}[X]$. Meanwhile,
$1-\delta_{\mathrm{est}}$ lower bounds the probability that we correctly infer
whether $\Psi_{t}(\boldsymbol{x})\geq\lambda$ based on a draw of
$\Psi_{t}^{n}(\boldsymbol{x})=\overline{X}_{n}$. The rationale for this
inference is explained below and previously appeared in Bardenet et al.
(2014).

Suppose we draw $\Psi_{t}^{n}(\boldsymbol{x})$. If this draw is closer to
$\Psi_{t}(\boldsymbol{x})$ than it is to the decision boundary $\lambda$, then
$\Psi_{t}^{n}(\boldsymbol{x})$ and $\Psi_{t}(\boldsymbol{x})$ must be
positioned on the same side of $\lambda$. Consequently, we would be able to
tell if $\Psi_{t}(\boldsymbol{x})\geq\lambda$ by looking at
$\Psi_{t}^{n}(\boldsymbol{x})$. According to Theorem 1, the distance
$\left|\Psi_{t}^{n}(\boldsymbol{x})-\Psi_{t}(\boldsymbol{x})\right|$ vanishes
as $n\to\infty$ with high probability. For now, assume we have chosen
$n\in\mathbb{N}$ in such a way that

$\left|\Psi_{t}^{n}(\boldsymbol{x})-\Psi_{t}(\boldsymbol{x})\right|<\left|\Psi_{t}^{n}(\boldsymbol{x})-\lambda\right|$
(9)

holds in at least $1-\delta_{\mathrm{est}}$ percent of cases.222We may assume
that $\Psi_{t}(\boldsymbol{x})\neq\lambda$ with probability one. Note that
$\Psi_{t}(\boldsymbol{x})$ is only random prior to making the first $t$
observations. Then, $\Psi_{t}^{n}(\boldsymbol{x})\geq\lambda$ implies the
probability that $\boldsymbol{x}$ is $\epsilon$-optimal is greater than or
equal to $(1-\delta_{\mathrm{mod}})(1-\delta_{\mathrm{est}})\geq 1-\delta$.

This line of reasoning is the basic premise of our method. We will generate
enough draws of $f_{t}$ so that our estimate is likely closer to
$\Psi_{t}(\boldsymbol{x})$ than it is to $\lambda$. We will then use this
estimate to decide whether or not $\Psi_{t}(\boldsymbol{x})\geq\lambda$.

Figure 2 illustrates how $\operatorname{PRB}$ behaves for different regret
bounds $\epsilon$ and risk tolerances $\delta$. Data was generated by running
BO a hundred times and jointly sampling $f^{*}_{t}$ and
$f_{t}(\boldsymbol{s}_{t})$ a thousand times at steps $t=5,\ldots,127$. These
samples were used to evaluate $\Psi_{t}^{n}(\boldsymbol{s}_{t},\epsilon)$ and
stopping was simulated by comparing these estimates with
$\lambda=1-\delta_{\mathrm{mod}}$, where
$\delta_{\mathrm{mod}}=\frac{\delta}{2}$. As $\epsilon$ and $\delta$ go to
zero, $\operatorname{PRB}$ becomes increasingly strict and stopping times are
seen to increase. Encouragingly, the success rates of stopped runs are
consistently close to or in excess of $1-\delta$. The same is not true when
estimates are directly compared to $1-\delta$, since this does not account for
estimation errors. Note that $n=1000$ samples is often not enough to ensure
correct decisions are made with high probability when
$\left|\Psi_{t}^{n}(\boldsymbol{s}_{t},\epsilon)-\lambda\right|$ is small.

A remaining question is how to choose sample sizes $n\in\mathbb{N}$. In
practice, drawing $f^{*}_{t}$ can often be fairly expensive. Our immediate
goal is therefore to design an algorithm that allows us to avoid paying for a
gratuitous number of samples.

##  4 Method

Figure 3:  _Left:_ Median number of draws used by Algorithm 1 to decide if the
expectation of a Bernoulli random variable $X\sim\operatorname{Bern}(p)$
exceeds $\lambda=10^{-5}$ (chosen arbitrarily). _Middle:_ Empirical CDFs of
$\Psi_{t}^{n}$ when running BO with known priors $\mathcal{GP}(0,k)$ in two,
four, and six dimensions given noise variance $\gamma^{2}=10^{-6}$ (solid,
$\circ$) or $\gamma^{2}=10^{-2}$ (dashed, $\times$). Here, $\epsilon=0.1$ and
$\delta_{\mathrm{mod}}=\delta_{\mathrm{est}}=0.025$. _Right:_ Runtimes of
Algorithm 1 using the sampling procedure described in Section 4 and an
approximate time limit of one thousand seconds.

This section shows how to convert the preceding ideas into practical
algorithms. At each time $t\in\mathbb{N}$, we will use a Monte Carlo estimate
(6) to test the $\operatorname{PRB}$ criterion. The number of samples used to
construct these estimates will be chosen online by Algorithm 1 (discussed
below) such that each estimate may be used to confidently infer whether a
point is $\epsilon$-optimal. We review this algorithm first and then give
details on how, when, and where to test $\operatorname{PRB}$.

Algorithm 1 operates by drawing i.i.d. copies of a random variable $X$ until a
bound on the resulting estimation error is sufficiently tight for it to
confidently decide whether or not $\operatorname{\mathbb{E}}[X]\geq\lambda$.
To better understand this, start by defining a pair of sequences: sample sizes
$(n_{j})$ and risk tolerances $(\delta_{j})$. The sizes should be increasing,
while the tolerances should be positive and satisfy
$\sum_{j=1}^{\infty}\delta_{j}\leq\delta$.

For every $j\in N$, use the first $n_{j}$ copies of $X$ to define an event
that holds with probability at least $1-\delta_{j}$

$\displaystyle\mathcal{E}_{j}=\left\\{\left|\overline{X}_{n_{j}}-\operatorname{\mathbb{E}}[X]\right|\leq\Delta_{n_{j}}\right\\},$
(10)

where $\Delta_{n_{j}}$ is the confidence bound given in Theorem 1. Each event
corresponds to a simple test: if $\mathcal{E}_{j}$ obtains and
$\left|\overline{X}_{n_{j}}-\lambda\right|\geq\Delta_{n_{j}}$,
$\overline{X}_{n_{j}}$ is on the same side of $\lambda$ as
$\operatorname{\mathbb{E}}[X]$. Suppose we generate draws of $X$ until we
encounter the first $j\in\mathbb{N}$ such that
$\left|\overline{X}_{n_{j}}-\lambda\right|>\Delta_{n_{j}}$ and use this
estimate to determine if we think that
$\operatorname{\mathbb{E}}[X]\geq\lambda$. By definition of $(\delta_{j})$ and
the union bound, we have

$\displaystyle\operatorname{\mathbb{P}}\left(\mathcal{E}_{1},\ldots,\mathcal{E}_{j}\right)\geq
1-\sum_{i=1}^{\infty}\operatorname{\mathbb{P}}(\mathcal{E}_{i}^{c})\geq
1-\delta.$ (11)

Hence, our decision is correct with probability at least $1-\delta$.

Algorithm 1 Adaptive Empirical Bernstein Classifier

1: Given a random variable $X:\Omega\to[a,b]$, a boundary $\lambda\in[a,b]$,
an increasing sequence of natural numbers $\left(n_{j}\right)$, and a positive
sequence
$\left(\delta_{j}\right)\mathrel{\text{missing}}{s.t.}\sum_{j=1}^{\infty}\delta_{j}\leq\delta$.

2: Initialize $j=0$, $n=0$, $\overline{\mu}=0$, $\overline{\nu}=0$, and
$\Delta=\infty$.

3: while $\left|\overline{\mu}-\lambda\right|<\Delta$ do

4: $j\leftarrow j+1$

5:
$\boldsymbol{x}\leftarrow\left[X(\omega_{n+1}),\ldots,X(\omega_{n_{j}})\right]^{\top}$
// Get $n_{j}-n$ draws

6: $\boldsymbol{y}\leftarrow\boldsymbol{x}-\overline{\mu}$

7: $n\leftarrow n_{j}$

8:
$\overline{\mu}\leftarrow\overline{\mu}+n^{-1}\boldsymbol{1}^{\top}\boldsymbol{y}$
// Update mean

9:
$\overline{\nu}\leftarrow\overline{\nu}+(\boldsymbol{x}-\overline{\mu})^{\top}\boldsymbol{y}$
// Update squared error

10: $\Delta\leftarrow
n^{-1}\left(\sqrt{2\log(3/\delta_{j})\overline{\nu}}+3(b-a)\log(3/\delta_{j})\right)$

11: end while

12: Return $\mathbbm{1}\left(\overline{\mu}\geq\lambda\right)$ //
$\operatorname{\mathbb{P}}\left[\mathbbm{1}\left(\overline{\mu}\geq\lambda\right)\neq\mathbbm{1}\left(\operatorname{\mathbb{E}}\left[X\right]\geq\lambda\right)\right]\leq\delta$

This algorithm is not new. Algorithm 1 was inspired by Mnih et al. (2008) and
references contained therein, who studied how concentration inequalities like
Theorem 1 can be used to iteratively test for convergence in the sense that
$\operatorname{\mathbb{P}}\left(\left|\overline{X}_{n}-\operatorname{\mathbb{E}}[X]\right|\leq\epsilon\operatorname{\mathbb{E}}[X]\right)\geq
1-\delta.$ As mentioned earlier however, our approach is closer to Bardenet et
al. (2014). There, a variant of Algorithm 1 was used to decide whether or not
to accept Metropolis-Hastings proposals based on subsampled log likelihood
estimates. This application is quite different from ours, but the underlying
task of inferring whether an expected value exceeds a level is the same.

We follow Mnih et al. (2008) by setting
$\delta_{j}=j^{-\alpha}\frac{(\alpha-1)}{\alpha}\delta$ where $\alpha=1.1$.
Given an initial sample size $n_{0}\in\mathbb{N}$, we similarly define
$n_{j}=\lceil\beta^{j-1}n_{0}\rceil$ where $\beta=1.5$. We set $n_{0}=64$ when
estimating $\Psi_{t}$. Using a geometric schedule for $(n_{j})$ helps to avoid
cases where $\delta_{j}$ shrinks unnecessarily due to a large number of tests
being performed with a small number of samples. In exchange, this schedule can
result in nearly $\beta$ times too many samples being drawn.

Figure 3 shows how Algorithm 1 behaves when deciding whether
$\operatorname{\mathbb{E}}[X]\geq\lambda$, where
$X\sim\operatorname{Bern}(p)$. As $p\to 1$, the distance between
$\overline{X}_{n}$ and $\lambda$ tends to increase and a decision can be made
with larger confidence bounds constructed using fewer draws of $X$. As
$\delta\to 1$, these bounds decrease and a decision can similarly be made
using fewer samples.

###  4.1 How to sample $X=\mathbbm{1}(r_{t}(\boldsymbol{x})\leq\epsilon)$

We propose to sample $f_{t}$ via pathwise conditioning Wilson et al. (2020;
2021). More specifically, we suggest to approximate the prior with a Bayesian
linear model

$\displaystyle\hat{f}(\cdot)=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$
$\displaystyle\boldsymbol{w}\sim\mathcal{N}(\boldsymbol{0},\mathbf{I}),$ (12)

where the feature map $\boldsymbol{\phi}:\mathcal{X}\to\mathbb{R}^{\ell}$ is
randomly generated so that
$\boldsymbol{\phi}(\boldsymbol{x})^{\top}\boldsymbol{\phi}(\boldsymbol{x}^{\prime})\approx
k(\boldsymbol{x},\boldsymbol{x}^{\prime})$ Rahimi & Recht (2007).

Denoting
$\mathbf{\Lambda}=k(\mathbf{X}_{t},\mathbf{X}_{t})+\gamma^{2}\mathbf{I}$ and
$\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},\gamma^{2}\mathbf{I})$,
we may then use Matheron’s rule to approximate draws from the posteriors as

$f_{t}(\cdot)\stackrel{{\scriptstyle
d}}{{\approx}}\hat{f}(\cdot)+k(\cdot,\mathbf{X}_{t})\mathbf{\Lambda}^{-1}\big{[}\boldsymbol{y}_{t}-\hat{f}(\mathbf{X}_{t})-\boldsymbol{\varepsilon}\big{]}.$
(13)

Given a draw of $f_{t}$, indicators
$\mathbbm{1}(r_{t}(\boldsymbol{x})\leq\epsilon)$ can then be obtained using
multi-start gradient ascent. A helpful insight is that generating indicators
with function draws eliminates the need to find $f^{*}_{t}$. For each draw of
$f_{t}$, it suffices to determine if there exists a point
$\boldsymbol{x}^{\prime}\in\mathcal{X}$ such that
$f_{t}(\boldsymbol{x}^{\prime})-f_{t}(\boldsymbol{x})>\epsilon$. This property
can be exploited to speed up the process of sampling
$\mathbbm{1}(r_{t}(\boldsymbol{x})\leq\epsilon)$, however these benefits
typically fall off as $\Psi_{t}(\boldsymbol{x})$ increases.

It should be said that this sampling procedure introduces a yet-to-be-
determined amount of error in practice, since draws of $f_{t}$ are not only
approximate but non-convex. Initial results suggest these errors are small
(see Figure 1), however we leave this as a topic for future investigation.

Figure 3 reports empirical distributions for
$\Psi_{t}^{n}(\boldsymbol{s}_{t})$ and wall times for Algorithm 1 under the
given sampling procedure. Runtimes were found to primarily depended on the
distance from $\Psi_{t}^{n}(\boldsymbol{s}_{t})$ to $\lambda$ and
dimensionality of $\mathcal{X}=[0,1]^{D}$. The latter trend largely occurred
because the number of starting positions we considered was exponential in $D$.
Overall, it is clear that letting Algorithm 1 pick sample sizes
$n\in\mathbb{N}$ leads to substantial savings over using a large static value.

###  4.2 When and where to test for $\epsilon$-optimality

At each step $t\in\mathbb{N}$ (or a subsequence thereof), we will use
Algorithm 1 to infer $(\epsilon,\delta_{\mathrm{mod}})$-optimality with
accepted risk less than or equal to $\delta_{\mathrm{est}}^{t}>0$. Defining
$\sum_{t=1}^{\infty}\delta_{\mathrm{est}}^{t}\leq\delta_{\mathrm{est}}$ allows
us to (union) bound the probability that any test fails and, hence, that BO
stops prematurely.

When testing multiple points per step, we must therefore divide
$\delta_{\mathrm{est}}^{t}$ in order to retain this bound. Draws of $f_{t}$
can be shared, but decreasing $\delta_{\mathrm{est}}^{t}$ may cause the
algorithm to request additional samples (see Figure 3). Results from the next
section suggest it suffices to test $\boldsymbol{s}_{t}$, but this practice
may lead to suboptimal stopping times. If solutions are restricted to
$\mathbf{X}_{t}$, then we may safely ignore any
$\boldsymbol{x}\in\mathbf{X}_{t}$ for which
$\operatorname{\mathbb{P}}\left(f_{t}(\boldsymbol{s}_{t})-f_{t}(\boldsymbol{x})\geq\epsilon\right)\geq\delta_{\mathrm{mod}}$.
In practice, we found that this heuristic usually eliminates most candidates.

More generally, we may use gradient information to maximize an estimator for
$\Psi_{t}$ based on one of the (approximate) conditioning strategies from
Figure 1. The resulting points can then be tested with an independent run of
Algorithm 1 (independent in the sense that draws of $f_{t}$ are not recycled).

##  5 Analysis

We will show that Bayesian optimization with the $\operatorname{PRB}$ stopping
criterion not only terminates, but does so correctly with high probability. We
begin by discussing the assumptions made throughout this section:

  1. A1.

The search space $\mathcal{X}=[0,1]^{D}$ is a unit hypercube.

  2. A2.

$f\sim\mathcal{GP}(0,k)$, where $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$
is a known covariance function with Lipschitz constant $L_{k}$, i.e.
$\left|k(\boldsymbol{x},\\!\boldsymbol{x})\\!-\\!k(\boldsymbol{x},\\!\boldsymbol{x}^{\prime})\right|\leq
L_{k}\left\lVert\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\rVert_{\infty}$
for all $\boldsymbol{x},\boldsymbol{x}^{\prime}\\!\in\\!\mathcal{X}$.

  3. A3.

The sequence $(\boldsymbol{x}_{t})$ is almost surely dense in $\mathcal{X}$.

A1 and A2 guarantee it is possible for the maximum posterior variance to
become arbitrarily small based on a finite number of observations. A3 then
implies that this event will almost surely occur in finite time. More
generally, A3 is necessary to ensure convergence when the only things we know
are that $\mathcal{X}$ is compact and $f$ is continuous Törn & Žilinskas
(1989).

When A1 and A2 hold, certain strategies are known to produce almost surely
dense iterates. For instance, Vazquez & Bect (2010) showed that Expected
Improvement Šaltenis (1971), Probability of Improvement Kushner (1962), and
Knowledge Gradient Močkus (1975); Frazier et al. (2009) all exhibit this
behavior when $f$ is observed. This condition can also be enforced by
introducing a small chance that $\boldsymbol{x}_{t}$ will be chosen at random
Sutton & Barto (2018).

We first prove that Bayesian optimization satisfies the $\operatorname{PRB}$
criterion with probability one when solutions may be freely chosen on
$\mathcal{X}$. An extension of this result for cases where solutions are
restricted to $\mathbf{X}_{t}$ is given in Appendix A.

######  Theorem 2.

Under assumptions A1–A3 and for all regret bounds $\epsilon>0$ and risk
tolerances $\delta>0$, there almost surely exists a $T\in\mathbb{N}$ so that,
at every time $t\geq T$, each point
$\boldsymbol{s}_{t}\in\operatorname*{\arg\max}_{\boldsymbol{x}\in\mathcal{X}}\mu_{t}(\boldsymbol{x})$
satisfies

$\displaystyle\operatorname{\mathbb{P}}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f_{t}\left(\boldsymbol{x}\right)-f_{t}\left(\boldsymbol{s}_{t}\right)\leq\epsilon\right]\geq
1-\delta.$ (14)

###### Proof.

See Appendix A. ∎

To understand the proof of Theorem 2, consider the process

$g_{t}(\cdot)=\underbracket{f_{t}(\cdot)-f_{t}(\boldsymbol{s}_{t})}_{\mathrel{\text{missing}}{performancegap}}+\underbracket{\mu_{t}(\boldsymbol{s}_{t})-\mu_{t}(\cdot)}_{\geq
0}.$ (15)

By construction, $g_{t}$ is centered and its supremum satisfies

$\displaystyle g_{t}^{*}$
$\displaystyle=\sup_{\boldsymbol{x}\in\mathcal{X}}g_{t}(\boldsymbol{x})\geq\sup_{\boldsymbol{x}\in\mathcal{X}}f_{t}(\boldsymbol{x})-f_{t}(\boldsymbol{s}_{t})=r_{t}(\boldsymbol{s}_{t}).$
(16)

Hence, it suffices to upper bound the probability that $g_{t}^{*}$ is greater-
equal to $\epsilon$. Such a bound may be obtain using the Borell-TIS
inequality Borell (1975); Tsirelson et al. (1976)

$\displaystyle\operatorname{\mathbb{P}}\left(g^{*}_{t}\geq\epsilon\right)\leq\exp\left(-\frac{1}{2}\left[\frac{\epsilon-\operatorname{\mathbb{E}}\left(g^{*}_{t}\right)}{2\sigma_{t}}\right]^{2}\right)$
(17)

which holds $\forall\epsilon>\operatorname{\mathbb{E}}\left(g^{*}_{t}\right)$,
where
$\sigma^{2}_{t}=\max_{\boldsymbol{x}\in\mathcal{X}}k_{t}(\boldsymbol{x},\boldsymbol{x})$.
In (17), we divide by $2\sigma_{t}$ since
$\max_{\boldsymbol{x}\in\mathcal{X}}\\!\operatorname{Var}\left[g_{t}(\boldsymbol{x})\right]\leq
4\sigma^{2}_{t}$. From here, one can prove that
$\operatorname{\mathbb{E}}(g_{t}^{*})$ decreases as $\sigma_{t}\to 0$. The
remaining task is therefore to show that $\sigma_{t}$ becomes arbitrarily
small as $\mathbf{X}_{t}$ becomes increasingly dense in $\mathcal{X}$.

Similar ideas can be found in Grünewälder et al. (2010), who proved that the
expected regret goes to zero under like assumptions. These authors showed that
BO will converge; we show that we will know it has done so. Equipped with this
knowledge, we now demonstrate that the algorithm stops successfully with high
probability.

######  Theorem 3.

Given a regret bound $\epsilon>0$ and risk tolerance $\delta>0$, let
$\Psi_{t}(\boldsymbol{s}_{t})$ be the probability under the model at time
$t\in\mathbb{N}$ that
$\boldsymbol{s}_{t}\in\operatorname*{\arg\max}_{\boldsymbol{x}\in\mathcal{X}}\mu_{t}(\boldsymbol{x})$
is $\epsilon$-optimal. Suppose that at each step $t$ (or a subsequence
thereof) an independent run of Algorithm 1 is used to produce a random,
unbiased estimate $\Psi_{t}^{n}(\boldsymbol{s}_{t})$ that is closer to
$\Psi_{t}(\boldsymbol{s}_{t})$ than it is to $1-\delta_{\mathrm{mod}}$ with
probability at least $1-\delta_{\mathrm{est}}^{t}$, where
$\delta_{\mathrm{mod}}$ and $\delta_{\mathrm{est}}$ are nonzero probabilities
so that $\delta_{\mathrm{mod}}+\delta_{\mathrm{est}}\leq\delta$ and
$(\delta_{\mathrm{est}}^{t})$ is a positive sequence satisfying
$\sum_{i=1}^{\infty}\delta_{\mathrm{est}}^{t}\leq\delta_{\mathrm{est}}$.

If assumptions A1–A3 hold and the first $\boldsymbol{s}_{t}$ for which
$\Psi_{t}^{n}(\boldsymbol{s}_{t})\geq 1-\delta_{\mathrm{mod}}$ is returned,
then the algorithm almost surely finishes in finite time with an
$(\epsilon,\delta)$-optimal solution.

###### Proof.

See Appendix A. ∎

The proof of Theorem 3 is straightforward. Since A1–A3 hold and estimates to
$\Psi_{t}(\boldsymbol{s}_{t})$ at each step are independent and unbiased,
Theorem 2 implies there almost surely exists a first time $T\in\mathbb{N}$
where $\Psi_{T}(\boldsymbol{s}_{T})\geq\lambda$ and
$\Psi_{T}^{n}(\boldsymbol{s}_{T})\geq\lambda$. The probability of terminating
prior to this point is at most
$\sum_{t=1}^{T-1}\delta_{\mathrm{est}}^{t}\leq\delta_{\mathrm{est}}$ by the
union bound, and the probability that $\boldsymbol{s}_{T}$ fails to be
$\epsilon$-optimal is at most $\delta_{\mathrm{mod}}$. As a result, the
algorithm stops and returns a solution that is $\epsilon$-optimal with
probability at least $1-\delta_{\mathrm{mod}}-\delta_{\mathrm{est}}\geq
1-\delta$.

Let us now see how well this theory holds up in practice.

Figure 4:  Results on functions drawn from known priors $\mathcal{GP}(0,k)$ in
two, four, and six dimensions with noise variance $\gamma^{2}=10^{-6}$ (solid)
or $\gamma^{2}=10^{-2}$ (dashed). _Left_ : Medians and interquartile ranges of
log regrets
$\log_{10}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})-f(\cdot)\right]$.
_Middle_ : Medians and interquartile ranges of estimated $\epsilon$-optimality
probabilities
$\max_{\boldsymbol{x}\in\mathbf{X}_{t}}\Psi_{t}^{n}(\boldsymbol{x})$. _Right:_
Empirical probabilities of evaluating $t$ or more trials. In the first two
panels, sequences were padded with their final values to accommodate different
lengths.

##  6 Experiments

Experiments discussed here were performed by first running BO with
conservatively chosen budgets $T\in\mathbb{N}$. We then stepped through each
saved run with different stopping rules to establish stopping times and
terminal performance. This paradigm ensured fair comparisons and reduced
overheads. For each problem, a hundred independent runs of BO were used to
generate the reported results. Note that despite the general notation of the
paper, all problems were defined as minimization tasks. Full details can be
found in Appendix B and extended results are presented in Appendix C. Code is
available online at https://github.com/j-wilson/trieste_stopping.

For $\operatorname{PRB}$, we divided $\delta$ evenly between
$\delta_{\mathrm{est}}$ and $\delta_{\mathrm{mod}}$. Since experiments were
carried out using preexisting BO runs that each began with five random trials
and ended at times $T$, we employed a constant schedule
$\delta_{\mathrm{est}}^{t}=\frac{1}{T-5}\delta_{\mathrm{est}}$ for the
accepted amount of risk at steps $t\in\mathbb{N}$. Parameter schedules for
Algorithm 1 are described in Section 4.

As a practical concession, we limited each run of Algorithm 1 to 1000 draws of
$f_{t}$ and used the resulting estimate to decide whether to stop—even if the
corresponding bound was not sufficiently tight to provide guarantees.
Preliminary results under this setup were consistent with those discussed in
Figure 3, so we opted to limit the number of samples in order to run more
experiments on a wider range of problems. Additionally, when optimizing
$f\sim\mathcal{GP}(0,k)$ in six dimensions with noise $\gamma^{2}=10^{-2}$, we
evaluated $\operatorname{PRB}$ once every five steps to further expedite these
experiments.

|  GP-2∗, 10-6 |  GP-2∗, 10-2 |  GP-4∗, 10-6 |  GP-4∗, 10-2 |  GP-6∗, 10-6 |  GP-6∗, 10-2
---|---|---|---|---|---|---
Oracle∗ |  $10\,(100\%)$ |  $11\,(100\%)$ |  $27\,(100\%)$ |  $29.5\,(100\%)$ |  $40.5\,(99\%)$ |  $65\,(100\%)$
Budget∗ |  $18\,(96\%)$ |  ${\color[rgb]{0,0,1}\mathbf{23\,(96\%)}}$ |  $65\,(95\%)$ |  $95\,(95\%)$ |  $125\,(95\%)$ |  ${\color[rgb]{0,0,1}\mathbf{228\,(96\%)}}$
Budget |  $33\,(98\%)$ |  $33\,(100\%)$ |  $65\,(95\%)$ |  $65\,(85\%)$ |  $97\,(86\%)$ |  $97\,(70\%)$
Acq. Cutoff∗ |  $26\,(99\%)$ |  $60.5\,(100\%)$ |  $76.5\,(100\%)$ |  $89\,(96\%)$ |  ${\color[rgb]{0,0,1}\mathbf{115.5\,(97\%)}}$ |  $145\,(92\%)$
UCB – LCB |  $12\,(89\%)$ |  $82\,(100\%)$ |  $22.5\,(66\%)$ |  $256\,(100\%)$ |  $32.5\,(50\%)$ |  $512\,(100\%)$
PRB (ours) |  ${\color[rgb]{0,0,1}\mathbf{17\,(97\%)}}$ |  ${\color[rgb]{0,0,1}\mathbf{23\,(99\%)}}$ |  ${\color[rgb]{0,0,1}\mathbf{64\,(99\%)}}$ |  ${\color[rgb]{0,0,1}\mathbf{86.5\,(96\%)}}$ |  $133.5\,(98\%)$ |  $235\,(100\%)$

|  GP-4, 10-6 |  GP-4, 10-2 |  Branin-2 |  Hartmann-3 |  Hartmann-6 | Rosenbrock-4
---|---|---|---|---|---|---
Oracle∗ |  $35\,(100\%)$ |  $51\,(100\%)$ |  $19\,(100\%)$ |  $14\,(100\%)$ |  $36\,(67\%)$ |  $34\,(100\%)$
Budget∗ |  $80\,(95\%)$ |  $158\,(95\%)$ |  ${\color[rgb]{0,0,1}\mathbf{26\,(95\%)}}$ |  $23\,(96\%)$ |  $256\,(67\%)$ |  ${\color[rgb]{0,0,1}\mathbf{47\,(95\%)}}$
Budget |  $65\,(85\%)$ |  $65\,(61\%)$ |  $33\,(100\%)$ |  $49\,(100\%)$ |  $97\,(67\%)$ |  $65\,(100\%)$
Acq. Cutoff∗ |  ${\color[rgb]{0,0,1}\mathbf{77.5\,(97\%)}}$ |  ${\color[rgb]{0,0,1}\mathbf{112\,(96\%)}}$ |  $64\,(100\%)$ |  $25\,(100\%)$ |  $38\,(67\%)$ |  $62\,(94\%)$
UCB – LCB |  $18\,(30\%)$ |  $224\,(80\%)$ |  $31\,(99\%)$ |  $15\,(83\%)$ |  $26\,(46\%)$ |  $73.5\,(100\%)$
PRB (ours) |  $61\,(88\%)$ |  $100.5\,(92\%)$ |  $33\,(99\%)$ |  ${\color[rgb]{0,0,1}\mathbf{19\,(100\%)}}$ |  $40\,(64\%)$ |  $84\,(100\%)$

Table 1:  Median stopping times and success rates when seeking
$(\epsilon,\delta)$-optimal solutions, where $\delta=0.05$ and $\epsilon=0.1$.
For Rosenbrock, $\epsilon=10^{-4}$. Numbers atop each column denote the
dimensionality of $\mathcal{X}=[0,1]^{D}$ followed by the noise variance
$\gamma^{2}$ (where present). Asterisks indicate that information which would
normally be unavailable was used to set parameters to their optimal values.
Non-oracle methods that returned $\epsilon$-optimal solutions at least
$1-\delta$ percent of the time using the fewest function evaluations are shown
in blue.

#### Baselines

We tested several baselines, many of which we granted access to information
that is usually unavailable (indicated by an asterisk). We summarize these as
follows:

  1. M1.

Oracle∗: stops once an $\epsilon$-optimal point is observed.

  2. M2.

Budget: stops after a predetermined number of trials.

  3. M3.

Budget∗: same as above, but with limits retroactively set to the smallest
values such that at least $1-\delta$ percent of runs returned
$\epsilon$-optimal results (where possible).

  4. M4.

Acquisition cutoff: stops once the acquisition value of the next trial falls
below a given threshold.

  5. M5.

UCB – LCB: stops when the gap between upper and lower confidence bounds is
less than or equal to $\epsilon$, i.e.

$\displaystyle\max_{\boldsymbol{x}\in\mathcal{X}}\operatorname{UCB}_{t}(\boldsymbol{x})-\max_{\boldsymbol{x}^{\prime}\in\mathbf{X}_{t}}\operatorname{LCB}_{t}(\boldsymbol{x}^{\prime})\leq\epsilon,$

where
$\operatorname{[U/L]CB}_{t}(\boldsymbol{x})=\mu_{t}(\boldsymbol{x})\pm\sqrt{\beta_{t}k_{t}(\boldsymbol{x},\boldsymbol{x})}$.

M1 is the optimal stopping rule, but requires perfect information for $f$.
Similarly, M3 is the optimal fixed budget in hindsight. For M2, we allocated
$16D$ trials per experiment in order to balance simplicity and performance.

M4 is discussed by Jones (2001), Nguyen et al. (2017), and references
contained therein for Expected Improvement, which is closely related to the
$1$-step optimal strategy used in our experiments. We stopped when acquisition
values fell below a threshold of $10^{-5}$, chosen post-hoc by walking powers
of ten and selecting the best performing value.

Lastly, M5 was proposed by Makarova et al. (2022). This method could seemingly
provide guarantees akin to those discussed in Section 5. However, the
suggested schedule
$\beta_{t}=\frac{2}{5}\log\left(Dt^{2}\pi^{2}/6\delta\right)$ does not afford
such a promise.333 Srinivas et al. (2009, Lemma 5.1) uses a union bound over
$\mathcal{X}$, so $\beta_{t}$ usually depends on cardinality
$\left|\mathcal{X}\right|$ not dimensionality $D$.

Below, we discuss results when optimizing draws from known priors
$f\sim\mathcal{GP}(0,k)$ and then proceed to the case where models were fit
online to the observed data.

###  6.1 Results with known models

Figure 4 shows how $\operatorname{PRB}$ runs evolved over time. Median
stopping times increased with dimensionality $D$ and noise variance
$\gamma^{2}$. While not evident in this plot, $\epsilon$-optimality
probabilities
$\max_{\boldsymbol{x}\in\mathbf{X}_{t}}\Psi_{t}^{n}(\boldsymbol{x})$ tended to
increase over time and occasionally dip back down when encountering better-
than-expected outcomes that caused the model to reconsider how $f$ behaves on
different regions of $\mathcal{X}$.

The top half of Table 1 summarizes the stopping times and success rates of
each method on synthetic problems. Here, $\operatorname{PRB}$ compared
favorably and either outperformed or kept pace with all (non-oracle)
baselines. As was true for the simulations presented in Section 3,
$\operatorname{PRB}$ success rates were consistently close to or above the
$95\%$ target rate.

M3 and M4 performed well in our experiments—save for on Branin, where M4 often
reached the limit of $T=64$ trials. This performance can largely be explained
for by the fact that we used their optimal parameter values. It is unclear how
these value should be set in practice. We found that M5 struggled to terminate
in the presence of moderate noise. This tendency can be understood by
comparing the various estimators in Figure 1. Decoupling $f^{*}_{t}$ and
$f_{t}(\boldsymbol{s}_{t})$ causes excess competition between
$f_{t}(\boldsymbol{s}_{t})$ and positively correlated process values
(including itself).

We also observe that the gap between $\operatorname{PRB}$ and oracle stopping
times often grew much faster than oracle stopping times themselves did. This
pattern emerges because we only need to learn about $f$ at an appropriately
chosen point to find an $\epsilon$-optimal solution. In contrast, we must
learn about $f$ globally to verify $\epsilon$-optimality. This interpretation
agrees with empirical analysis in which we found that $\operatorname{PRB}$
stopping times correlated more with how function values concentrated around
$f^{*}$ than with the first time at which an $\epsilon$-optimal point was
observed. We speculate that strategies which better account for the goal of
finding an $(\epsilon,\delta)$-optimal point may help to bridge this gap.
Similar ideas were previously explored by Cai et al. (2021).

###  6.2 Results with maximum likelihood models

Results for this setting occupy the bottom half of Table 1. As seen in
Appendix C, BO performance was consistently strong. However, this performance
did not always translate to the satisfactory outcomes at stopping time.

On some tasks (Branin, Hartmann-3 and Rosenbrock-4), models were sufficiently
well-calibrated and model-based stopping rules generally performed well. On
others (GP-4 and Hartmann-6), overconfident models caused optimazation runs or
model-based stopping to miss the mark. For 4D functions drawn from GP priors,
the variance of $f$ was often underestimated. For Hartmann-6, the smoothness
of $f$ was often overestimated. It is worth mentioning that $33\%$ of the
underlying optimization runs on Hartmann-6 converged to a local minimum valued
at $-3.20$ rather than the global minimum of $-3.32$ (see Appendix C). Hence,
each of these runs and the vast majority of stopping rule experiments would
have succeeded if $\epsilon=0.1$ was slightly larger.

These results highlight the importance of uncertainty calibration. Our
hyperpriors were chosen to avoid strong inductive biases and this flexibility
increased the chances of errant model fits. In practice, we suggest using
conservative priors that favor smaller lengthscales and larger variances.
Calibration may also be improved by marginalizing over hyperparameters Simpson
et al. (2021) or utilizing more expressive models, such as deep GPs Salimbeni
& Deisenroth (2017) or Bayesian neural networks Li et al. (2023). Note that
while our analysis is specific to GPs, $\operatorname{PRB}$ only assumes that
we can sample $\epsilon$-optimality indicators.

##  7 Discussion

To the best of our knowledge, results presented here are among the first of
their kind for Bayesian optimization. Supposing the prior distribution of
problems is known, we have shown that the proposed stopping criterion
($\operatorname{PRB}$) guarantees BO will return a solution that the user
finds satisfactory with high probability. Further, we have demonstrated
empirically that these findings carry over in practice for well-calibrated
models.

Nevertheless, there is still more work to be done. These same empirical
results echo initial skepticism about placing too much faith in our models.
The success of the proposed stopping rule (and others like it) will depend in
large part on our ability to overcome this challenge. To the extent that
models prove trustworthy however, we suggest that model-based stopping may one
day become as common place as model-based optimization.

## Impact Statement

The goal of this work is to introduce a model-based stopping rule for Bayesian
optimization that saves resources and alleviates the burden on the user. As is
usually the case, automating real-world decisions comes with a certain degree
of risk. Still, we remain optimistic that our work will have a positive impact
on theoretical and applied fronts.

## References

  * Abadi et al. (2015) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G. S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., and Zheng, X.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Adler & Taylor (2009) Adler, R. J. and Taylor, J. E.  _Random fields and geometry_.  Springer Science & Business Media, 2009.
  * Audibert et al. (2007) Audibert, J.-Y., Munos, R., and Szepesvári, C.  Tuning bandit algorithms in stochastic environments.  In _International conference on algorithmic learning theory_ , pp. 150–165. Springer, 2007.
  * Balandat et al. (2020) Balandat, M., Karrer, B., Jiang, D. R., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  In _Advances in Neural Information Processing Systems 33_ , 2020.  URL http://arxiv.org/abs/1910.06403.
  * Bardenet et al. (2014) Bardenet, R., Doucet, A., and Holmes, C.  Towards scaling up Markov chain Monte Carlo: an adaptive subsampling approach.  In _International conference on machine learning_ , pp. 405–413. PMLR, 2014.
  * Bishop (2006) Bishop, C.  Pattern recognition and machine learning.  _Springer google schola_ , 2:531–537, 2006.
  * Borell (1975) Borell, C.  The Brunn-Minkowski inequality in Gauss space.  _Inventiones mathematicae_ , 30(2):207–216, 1975\.
  * Byrd et al. (1995) Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C.  A limited memory algorithm for bound constrained optimization.  _SIAM Journal on scientific computing_ , 16(5):1190–1208, 1995.
  * Cai et al. (2021) Cai, X., Gomes, S., and Scarlett, J.  Lenient regret and good-action identification in Gaussian process bandits.  In _International Conference on Machine Learning_ , pp. 1183–1192. PMLR, 2021.
  * Frazier et al. (2009) Frazier, P., Powell, W., and Dayanik, S.  The knowledge-gradient policy for correlated normal beliefs.  _INFORMS journal on Computing_ , 21(4):599–613, 2009.
  * Garnett (2023) Garnett, R.  _Bayesian optimization_.  Cambridge University Press, 2023.
  * Griffiths & Hernández-Lobato (2020) Griffiths, R.-R. and Hernández-Lobato, J. M.  Constrained Bayesian optimization for automatic chemical design using variational autoencoders.  _Chemical science_ , 11(2):577–586, 2020.
  * Grünewälder et al. (2010) Grünewälder, S., Audibert, J.-Y., Opper, M., and Shawe-Taylor, J.  Regret bounds for Gaussian process bandit problems.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pp. 273–280. JMLR Workshop and Conference Proceedings, 2010.
  * Hansen (2016) Hansen, N.  The CMA evolution strategy: A tutorial.  _arXiv preprint arXiv:1604.00772_ , 2016.
  * Jones (2001) Jones, D. R.  A taxonomy of global optimization methods based on response surfaces.  _Journal of global optimization_ , 21:345–383, 2001.
  * Kushner (1962) Kushner, H. J.  A versatile stochastic model of a function of unknown and time varying form.  _Journal of Mathematical Analysis and Applications_ , 5(1):150–167, 1962.
  * Lederer et al. (2019a) Lederer, A., Umlauft, J., and Hirche, S.  Posterior variance analysis of Gaussian processes with application to average learning curves.  _arXiv preprint arXiv:1906.01404_ , 2019a.
  * Lederer et al. (2019b) Lederer, A., Umlauft, J., and Hirche, S.  Uniform error bounds for Gaussian process regression with application to safe control.  _Advances in Neural Information Processing Systems_ , 32, 2019b.
  * Li et al. (2023) Li, Y. L., Rudner, T. G., and Wilson, A. G.  A study of Bayesian neural network surrogates for Bayesian optimization.  _arXiv preprint arXiv:2305.20028_ , 2023.
  * Makarova et al. (2022) Makarova, A., Shen, H., Perrone, V., Klein, A., Faddoul, J. B., Krause, A., Seeger, M., and Archambeau, C.  Automatic termination for hyperparameter optimization.  In _International Conference on Automated Machine Learning_ , pp. 7–1. PMLR, 2022.
  * Massart (2007) Massart, P.  _Concentration inequalities and model selection: Ecole d’Eté de Probabilités de Saint-Flour XXXIII-2003_.  Springer, 2007.
  * Matthews et al. (2017) Matthews, A. G. d. G., van der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León-Villagrá, P., Ghahramani, Z., and Hensman, J.  GPflow: A Gaussian process library using TensorFlow.  _Journal of Machine Learning Research_ , 18(40):1–6, apr 2017.  URL http://jmlr.org/papers/v18/16-537.html.
  * Maurer & Pontil (2009) Maurer, A. and Pontil, M.  Empirical Bernstein bounds and sample variance penalization.  _arXiv preprint arXiv:0907.3740_ , 2009.
  * Mnih et al. (2008) Mnih, V., Szepesvári, C., and Audibert, J.-Y.  Empirical Bernstein stopping.  In _Proceedings of the 25th international conference on Machine learning_ , pp. 672–679, 2008.
  * Močkus (1975) Močkus, J.  On Bayesian methods for seeking the extremum.  In _Optimization Techniques IFIP Technical Conference: Novosibirsk, July 1–7, 1974_ , pp. 400–404. Springer, 1975.
  * Nguyen et al. (2017) Nguyen, V., Gupta, S., Rana, S., Li, C., and Venkatesh, S.  Regret for expected improvement over the best-observed value and stopping condition.  In _Asian conference on machine learning_ , pp. 279–294. PMLR, 2017\.
  * Picheny et al. (2023) Picheny, V., Berkeley, J., Moss, H. B., Stojic, H., Granta, U., Ober, S. W., Artemev, A., Ghani, K., Goodall, A., Paleyes, A., Vakili, S., Pascual-Diaz, S., Markou, S., Qing, J., Loka, N. R. B. S., and Couckuyt, I.  Trieste: Efficiently exploring the depths of black-box functions with tensorflow, 2023.  URL https://arxiv.org/abs/2302.08436.
  * Rahimi & Recht (2007) Rahimi, A. and Recht, B.  Random features for large-scale kernel machines.  _Advances in neural information processing systems_ , 20, 2007.
  * Salimbeni & Deisenroth (2017) Salimbeni, H. and Deisenroth, M.  Doubly stochastic variational inference for deep Gaussian processes.  _Advances in neural information processing systems_ , 30, 2017.
  * Šaltenis (1971) Šaltenis, V. R.  On a method of multi-extremal optimization.  _Automatics and Computers (Avtomatika i Vychislitelnayya Tekchnika)_ , 3:33–38, 1971.
  * Scott et al. (2011) Scott, W., Frazier, P., and Powell, W.  The correlated knowledge gradient for simulation optimization of continuous parameters using Gaussian process regression.  _SIAM Journal on Optimization_ , 21(3):996–1026, 2011.
  * Simon (1956) Simon, H. A.  Rational choice and the structure of the environment.  _Psychological review_ , 63(2):129, 1956.
  * Simpson et al. (2021) Simpson, F., Lalchand, V., and Rasmussen, C. E.  Marginalised Gaussian processes with nested sampling.  _Advances in neural information processing systems_ , 34:13613–13625, 2021.
  * Snoek et al. (2012) Snoek, J., Larochelle, H., and Adams, R. P.  Practical Bayesian optimization of machine learning algorithms.  _Advances in neural information processing systems_ , 25, 2012.
  * Srinivas et al. (2009) Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M.  Gaussian process optimization in the bandit setting: No regret and experimental design.  _arXiv preprint arXiv:0912.3995_ , 2009.
  * Sutton & Barto (2018) Sutton, R. S. and Barto, A. G.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Törn & Žilinskas (1989) Törn, A. and Žilinskas, A.  _Global optimization_ , volume 350.  Springer, 1989.
  * Tsirelson et al. (1976) Tsirelson, B., Ibragimov, I., and Sudakov, V.  Norms of Gaussian sample functions.  In _Proceedings of the Third Japan—USSR Symposium on Probability Theory_ , volume 550, pp. 20–41. Tashkent, 1976.
  * Vazquez & Bect (2010) Vazquez, E. and Bect, J.  Convergence properties of the expected improvement algorithm with fixed mean and covariance functions.  _Journal of Statistical Planning and inference_ , 140(11):3088–3095, 2010.
  * von Kügelgen et al. (2019) von Kügelgen, J., Rubenstein, P. K., Schölkopf, B., and Weller, A.  Optimal experimental design via Bayesian optimization: active causal structure learning for Gaussian process networks.  _arXiv preprint arXiv:1910.03962_ , 2019.
  * Wilson et al. (2018) Wilson, J., Hutter, F., and Deisenroth, M.  Maximizing acquisition functions for Bayesian optimization.  _Advances in neural information processing systems_ , 31, 2018.
  * Wilson et al. (2020) Wilson, J., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M.  Efficiently sampling functions from Gaussian process posteriors.  In _International Conference on Machine Learning_ , pp. 10292–10302. PMLR, 2020.
  * Wilson et al. (2021) Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P.  Pathwise conditioning of Gaussian processes.  _The Journal of Machine Learning Research_ , 22(1):4741–4787, 2021.

##  Appendix A Technical proofs

Findings presented in the body of this paper rely upon a number of supporting
results, some of which borrow heavily from earlier works. The main results of
this section are as follows. Lemma 5 shows that the posterior variance becomes
arbitrarily small as $\mathbf{X}_{t}$ becomes increasingly dense in
$\mathcal{X}$. Lemma 7 upper bounds the expected supremum of
$f\sim\mathcal{GP}(0,k)$ based on its maximum variance. Theorem 2 and
Corollary 8 prove that the $\operatorname{PRB}$ stopping criterion (5) almost
surely converges. Finally, Theorem 3 shows that Bayesian optimization using
$\operatorname{PRB}$ terminates and returns an $(\epsilon,\delta)$-optimal
solution.

Where relevant, we will attribute credit at the beginning of the corresponding
proof.

######  Proposition 4.

Let $\mathcal{X}\subseteq\mathbb{R}^{D}$ be convex and suppose that
$\mathbf{X}\subseteq\mathcal{X}$ generates an $\varepsilon$-cover of
$\mathcal{X}$. For every $\boldsymbol{x}\in\mathcal{X}$ and
$\rho\geq\varepsilon$, the intersection of the set $\mathbf{X}$ and the ball
$B(\boldsymbol{x},\rho)=\left\\{\boldsymbol{x}^{\prime}\in\mathcal{X}:\left\lVert\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\rVert_{\infty}\leq\rho\right\\}$
generates a $2\varepsilon$-cover of $B(\boldsymbol{x},\rho)$.

###### Proof.

Consider the ball $B(\boldsymbol{x},r)$ with radius $r=\rho-\varepsilon$.
Since $\mathcal{X}$ is convex, for every point $\boldsymbol{a}\in
B(\boldsymbol{x},\rho)$ there exists a $\boldsymbol{b}\in B(\boldsymbol{x},r)$
such that
$\left\lVert\boldsymbol{a}-\boldsymbol{b}\right\rVert_{\infty}\leq\varepsilon$.
Moreover, because $\mathbf{X}$ generates an $\varepsilon$-cover of
$\mathcal{X}$, for every point $\boldsymbol{b}\in B(\boldsymbol{x},r)$ there
exists a $\boldsymbol{c}\in\mathbf{X}$ so that
$\left\lVert\boldsymbol{b}-\boldsymbol{c}\right\rVert_{\infty}\leq\varepsilon$,
which implies that $\boldsymbol{c}\in B(\boldsymbol{x},\rho)$. It follows by
the triangle inequality that for every point $\boldsymbol{a}\in
B(\boldsymbol{x},\rho)$ there exists a pair of points
$\boldsymbol{b},\boldsymbol{c}\in
B(\boldsymbol{x},r)\times\left[B(\boldsymbol{x},\rho)\cap\mathbf{X}\right]$
such that

$\left\lVert\boldsymbol{a}-\boldsymbol{c}\right\rVert_{\infty}\leq\left\lVert\boldsymbol{a}-\boldsymbol{b}\right\rVert_{\infty}+\left\lVert\boldsymbol{b}-\boldsymbol{c}\right\rVert_{\infty}\leq\varepsilon+\varepsilon=2\varepsilon,$
(18)

which completes the proof. ∎

######  Lemma 5.

Under assumptions A1 and A2, if
$y(\cdot)\sim\mathcal{N}\left(f(\cdot),\gamma^{2}\right)$ is observed on a set
of points $\mathbf{X}\subseteq\mathcal{X}$ that generates an
$\varepsilon$-cover of $\mathcal{X}$,
$0\leq\varepsilon\leq\min\left\\{1,\nicefrac{{k(\boldsymbol{x},\boldsymbol{x})}}{{L_{k}}}\right\\},$
then

$\displaystyle\operatorname{Var}\left[f(\boldsymbol{x})\mid
y(\mathbf{X})\right]\leq\kappa_{\varepsilon}(\boldsymbol{x}),$ (19)

where

$\kappa_{\varepsilon}(\boldsymbol{x})=\frac{\left[4L_{k}\rho(\varepsilon)k(\boldsymbol{x},\boldsymbol{x})-L_{k}^{2}\rho(\varepsilon)^{2}\right]\eta(\varepsilon)+\gamma^{2}k(\boldsymbol{x},\boldsymbol{x})}{\left[k(\boldsymbol{x},\boldsymbol{x})+2L_{k}\rho(\varepsilon)\right]\eta(\varepsilon)+\gamma^{2}}$

is given in terms of
$\eta(\varepsilon)=\max\left\\{1,\nicefrac{{\rho(\varepsilon)}}{{4\varepsilon}}\right\\}^{D}$
and $\rho(\varepsilon)=\varepsilon^{\varepsilon}$ for any $0<\varepsilon<1$.

###### Proof.

This result extends Lederer et al. (2019a, Theorem 3.1), who showed that, for
all $0\leq\rho\leq\nicefrac{{k(\boldsymbol{x},\boldsymbol{x})}}{{L_{k}}}$,

$\displaystyle\operatorname{Var}\left[f(\boldsymbol{x})\mid
y\left(\mathbf{B}_{\rho}(\boldsymbol{x})\right)\right]\leq\frac{\left(4L_{k}\rho
k(\boldsymbol{x},\boldsymbol{x})-L_{k}^{2}\rho^{2}\right)\left|\mathbf{B}_{\rho}(\boldsymbol{x})\right|+\gamma^{2}k(\boldsymbol{x},\boldsymbol{x})}{\left(k(\boldsymbol{x},\boldsymbol{x})+2L_{k}\rho\right)\left|\mathbf{B}_{\rho}(\boldsymbol{x})\right|+\gamma^{2}}$
$\displaystyle\mathbf{B}_{\rho}(\boldsymbol{x})=B(\boldsymbol{x},\rho)\cap\mathbf{X}.$
(20)

We would like to translate this upper bound as a function of
$0\leq\varepsilon\leq 1$. To this end, begin by noticing that the upper bound
in (20) increases monotonically on
$0\leq\rho\leq\nicefrac{{k(\boldsymbol{x},\boldsymbol{x})}}{{L_{k}}}$ and
decreases monotonically on
$n=\left|\mathbf{B}_{\rho}(\boldsymbol{x})\right|\in\mathbb{W}$. Substituting
$\rho(\varepsilon)$ for $\rho$ and $\eta(\varepsilon)$ for $n$ therefore
yields a valid bound so long as
$\rho\leq\rho(\varepsilon)\leq\nicefrac{{k(\boldsymbol{x},\boldsymbol{x})}}{{L_{k}}}$
and $0\leq\eta(\varepsilon)\leq n$. As in (20), $\rho(\varepsilon)$ defines
the radius of a ball around $\boldsymbol{x}$ and $\eta(\varepsilon)$ conveys
the number of elements from $\mathbf{X}$ that lie within this ball.

Starting with the latter, lower bounds on the cardinality of
$\mathbf{B}_{\rho}(\boldsymbol{x})$ may be obtained from the fact that
$\mathbf{X}$ is assumed to generate an $\varepsilon$-cover of $\mathcal{X}$.
By Proposition 4, it follows that $\mathbf{B}_{\rho}(\boldsymbol{x})$
generates a $2\varepsilon$-cover of $B(\boldsymbol{x},\rho)$. Accordingly,
$\left|\mathbf{B}_{\rho}(\boldsymbol{x})\right|$ must be greater-equal to the
minimum number of points required to construct such a cover. Under the
$\left\lVert\cdot\right\rVert_{\infty}$ norm, the $\varepsilon$-covering
number of a ball

$B(\boldsymbol{x},\rho)=\prod_{d=1}^{D}\left[\max(x_{d}-\rho,0),\min(x_{d}+\rho,1)\right]$
(21)

is given by

$M\big{(}B(\boldsymbol{x},\rho),\left\lVert\cdot\right\rVert_{\infty},\varepsilon\big{)}=\prod_{d=1}^{D}\bigg{\lceil}\frac{\min(x_{d}+\rho,1)-\max(0,x_{d}-\rho)}{2\varepsilon}\bigg{\rceil}.$
(22)

This number is minimized when $B(\cdot,\rho)$ is placed in a corner, such as
$B(\boldsymbol{0},\rho)=[0,\rho]^{D}$. Choosing

$\eta(\varepsilon)=1\vee\left[\frac{\rho(\varepsilon)}{4\varepsilon}\right]^{D}\leq\bigg{\lceil}\frac{\rho(\varepsilon)}{4\varepsilon}\bigg{\rceil}^{D}$
(23)

therefore ensures that $\eta(\varepsilon)$ lower bounds the cardinality of
every $\mathbf{B}_{\rho}(\cdot)$. Note that there are two factors of two at
play here: one accounts for the fact that $\mathbf{B}_{\rho}(\cdot)$ is only
guaranteed to provide a $2\varepsilon$-cover of $B(\cdot,\rho)$, and the other
accounts for the fact that the corner balls are up to $2^{D}$ times smaller
than other balls with the same radius.

Turning our attention to the choice of function $\rho(\varepsilon)$, some
desiderata come into focus. First, we require
$\rho(\varepsilon)\geq\varepsilon$ so that every $\mathbf{B}_{\rho}(\cdot)$ is
nonempty. Second, we desire $\lim_{\varepsilon\to 0^{+}}\rho(\varepsilon)=0$
because the resulting posterior variance bound will increase monotonically in
$\rho(\varepsilon)$. Lastly, we want the ratio of $\rho(\varepsilon)$ to
$\varepsilon$ to diverge to infinity as $\varepsilon$ approaches zero from
above so that $\lim_{\varepsilon\to 0^{+}}\eta(\varepsilon)=\infty$. Based on
these criteria, a convenient choice when $\mathcal{X}=[0,1]^{D}$ is

$\displaystyle\rho(\varepsilon)=\varepsilon^{\alpha}$ $\displaystyle
0<\alpha<1.$ (24)

In summary, the claim follows by expressing $\rho$ as a function of
$\varepsilon$ and using it to lower bound
$\left|\mathbf{B}_{\rho}(\cdot)\right|$ with $\eta(\varepsilon)$:

$\operatorname{Var}\left[f(\boldsymbol{x})\mid
y(\mathbf{X})\right]\leq\frac{\left(4L_{k}\rho(\varepsilon)k(\boldsymbol{x},\boldsymbol{x})-L_{k}^{2}\rho(\varepsilon)^{2}\right)\eta(\varepsilon)+\gamma^{2}k(\boldsymbol{x},\boldsymbol{x})}{\left(k(\boldsymbol{x},\boldsymbol{x})+2L_{k}\rho(\varepsilon)\right)\eta(\varepsilon)+\gamma^{2}}.$
(25)

∎

######  Proposition 6.

For any choice of constants $a>0$, $b\geq 0$, $c\geq 0$,

$\int_{0}^{c}\sqrt{\log(1+b\varepsilon^{-\nicefrac{{1}}{{a}}})}d\varepsilon\leq
c\sqrt{a^{-1}+\log\left(1+bc^{-\nicefrac{{1}}{{a}}}\right)}.$ (26)

###### Proof.

The following proof originally appeared in Grünewälder et al. (2010, Appendix
A). Let $\xi=(1+\sqrt[a]{c}b^{-1})^{a}$ so that

$\displaystyle\int_{0}^{c}\sqrt{\log\left(1+b\varepsilon^{-\nicefrac{{1}}{{a}}}\right)}d\varepsilon\leq\int_{0}^{c}\sqrt{\log\left(\xi^{\nicefrac{{1}}{{a}}}b\varepsilon^{-\nicefrac{{1}}{{a}}}\right)}d\varepsilon.$
(27)

Next, define auxiliary functions

$\displaystyle f(u)=\sqrt{\log\left(u^{-\nicefrac{{1}}{{a}}}\right)}$
$\displaystyle g(\varepsilon)=\frac{\varepsilon}{\xi b^{a}}$ (28)

such that
$f(g(\varepsilon))=\sqrt{\log\left(\xi^{\nicefrac{{1}}{{a}}}b\varepsilon^{-\nicefrac{{1}}{{a}}}\right)}$,
and use them to integrate by substitution:

$\displaystyle\int_{0}^{c}\sqrt{\log\left(\xi^{\nicefrac{{1}}{{a}}}b\varepsilon^{-\nicefrac{{1}}{{a}}}\right)}d\varepsilon=\xi
b^{a}\int_{0}^{g(c)}\sqrt{\log\left(u^{-\nicefrac{{1}}{{a}}}\right)}du=\frac{\xi
b^{a}}{\sqrt{a}}\int_{0}^{g(c)}\sqrt{-\log(u)}du.$ (29)

The Cauchy-Schwarz inequality now gives

$\displaystyle\int_{0}^{g(c)}\sqrt{-\log(u)}du\leq\left(\int_{0}^{g(c)}du\right)^{\nicefrac{{1}}{{2}}}\left(-\int_{0}^{g(c)}\log(u)du\right)^{\nicefrac{{1}}{{2}}}=\frac{c}{\xi
b^{a}}\sqrt{1-\log\left(\frac{c}{\xi b^{a}}\right)}.$ (30)

Hence, the claim follows

$\displaystyle\int_{0}^{c}\sqrt{\log\left(1+b\varepsilon^{-\nicefrac{{1}}{{a}}}\right)}d\varepsilon\leq
c\sqrt{\frac{1+\log\left(\xi
b^{a}c^{-1}\right)}{a}}=c\sqrt{a^{-1}+\log\left(1+bc^{-\nicefrac{{1}}{{a}}}\right)}.$
(31)

For comparison with Grünewälder et al. (2010), we have

$\displaystyle
b^{a}=2c\implies\xi=\left(1+2^{-\nicefrac{{1}}{{a}}}\right)^{a}\leq
2^{a}\implies c\sqrt{\frac{1+\log\left(\xi
b^{a}c^{-1}\right)}{a}}=c\sqrt{\frac{1+\log\left(2\xi\right)}{a}}\leq
c\sqrt{\frac{\log\left(e2^{a+1}\right)}{a}}.$ (32)

∎

######  Lemma 7.

Let $f\sim\mathcal{GP}(0,k)$ be a Gaussian process with an $L_{k}$-Lipschitz
continuous covariance function $k:\mathcal{X}^{2}\to\mathbb{R}$ on
$\mathcal{X}=[0,r]^{D}$ having maximum variance
$\sigma^{2}=\max_{\boldsymbol{x}\in\mathcal{X}}k(\boldsymbol{x},\boldsymbol{x})$.
Then,

$\operatorname{\mathbb{E}}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})\right]\leq
12\sigma\sqrt{2D+D\log\left(1+4L_{k}r\sigma^{-2}\right)}.$ (33)

###### Proof.

This proof paraphrases parts of Grünewälder et al. (2010, Section 4.3).

Massart (2007, Theorem 3.18) proved that the expected supremum of $f$ is upper
bounded by

$\operatorname{\mathbb{E}}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})\right]\leq
12\int_{0}^{\sigma}\sqrt{\log N(\mathcal{X},d_{k},\varepsilon)}d\varepsilon,$
(34)

where $N(\mathcal{X},d_{k},\varepsilon)$ is defined as the
$\varepsilon$-packing number—i.e. the largest number of points that can be
“packed” inside of $\mathcal{X}$ without any two points being within
$\varepsilon$ of one another—under the canonical pseudo-metric444While $d_{k}$
has most of the properties of a proper metric,
$d_{k}(\boldsymbol{x},\boldsymbol{x}^{\prime})=0$ does not always imply that
$\boldsymbol{x}=\boldsymbol{x}^{\prime}$ Adler & Taylor (2009).

$d_{k}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\operatorname{\mathbb{E}}\left[(f(\boldsymbol{x})-f(\boldsymbol{x}^{\prime}))^{2}\right]^{\nicefrac{{1}}{{2}}}=\sqrt{k(\boldsymbol{x},\boldsymbol{x})-2k(\boldsymbol{x},\boldsymbol{x}^{\prime})+k(\boldsymbol{x}^{\prime},\boldsymbol{x}^{\prime})}.$
(35)

We may use (34) by upper bounding the right-hand side with a known quantity.
We will bound the $\epsilon$-packing number
$N(\mathcal{X},d_{k},\varepsilon)$, translate this bound from the $d_{k}$
pseudo-metric to the infinity norm, and then integrate the result.

The first step follows immediately from the fact that the
$\varepsilon$-packing number is smaller than the
$\tfrac{\varepsilon}{2}$-covering number—defined as the minimum number of
balls $B(\cdot,\tfrac{\varepsilon}{2})$ required to cover $\mathcal{X}$. The
second is accomplished by using Lipschitz continuity of $k$ to show that the
squared pseudo-metric $d_{k}(\cdot,\cdot)^{2}$ is $2L_{k}$-Lipschitz:

$\displaystyle
d_{k}(\boldsymbol{x},\boldsymbol{x}^{\prime})^{2}=\big{[}k(\boldsymbol{x},\boldsymbol{x})-k(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{]}+\big{[}k(\boldsymbol{x}^{\prime},\boldsymbol{x}^{\prime})-k(\boldsymbol{x}^{\prime},\boldsymbol{x})\big{]}\leq
2L_{k}\left\lVert\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\rVert_{\infty}$
$\displaystyle\forall\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}.$
(36)

Given a set $\mathbf{X}\subseteq\mathcal{X}$ such that
$\max_{\boldsymbol{x}\in\mathcal{X}}\min_{\boldsymbol{x}^{\prime}\in\mathbf{X}}\left\lVert\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\rVert_{\infty}\leq
C$, it follows that
$\max_{\boldsymbol{x}\in\mathcal{X}}\min_{\boldsymbol{x}^{\prime}\in\mathbf{X}}d_{k}(\boldsymbol{x},\boldsymbol{x}^{\prime})\leq\sqrt{2L_{k}C}$.
An $\nicefrac{{\varepsilon^{2}}}{{8L_{k}}}$-cover under the infinity norm
therefore guarantees an $\nicefrac{{\varepsilon}}{{2}}$-cover under $d_{k}$.
The former may be constructed from a grid of uniformly spaced points with
elements at intervals of $\nicefrac{{\varepsilon^{2}}}{{4L_{k}}}$. Such a grid
will consist of $\big{\lceil}4L_{k}r\varepsilon^{-2}\big{\rceil}^{D}$ points
assuming $\mathcal{X}=[0,r]^{D}$, meaning that

$N(\mathcal{X},d_{k},\varepsilon)<\left(1+4L_{k}r\varepsilon^{-2}\right)^{D}.$
(37)

To complete the proof, use Proposition 6 with $a=\tfrac{1}{2}$, $b=4L_{k}r$,
and $c=\sigma$ to show that

$\operatorname{\mathbb{E}}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})\right]\leq
12\sqrt{D}\int_{0}^{\sigma}\sqrt{\log\left(1+4L_{k}r\varepsilon^{-2}\right)}d\varepsilon\leq
12\sigma\sqrt{2D+D\log\left(1+4L_{k}r\sigma^{-2}\right)}.$ (38)

∎

See 2

###### Proof.

Consider the centered process

$g_{t}(\cdot)=\big{[}f_{t}(\cdot)-\mu_{t}(\cdot)\big{]}-\big{[}f_{t}(\boldsymbol{s}_{t})-\mu_{t}(\boldsymbol{s}_{t})\big{]},$
(39)

with covariance

$\displaystyle c_{t}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=k_{t}(\boldsymbol{x},\boldsymbol{x}^{\prime})-k_{t}(\boldsymbol{x},\boldsymbol{s}_{t})-k_{t}(\boldsymbol{s}_{t},\boldsymbol{x}^{\prime})+k_{t}(\boldsymbol{s}_{t},\boldsymbol{s}_{t}).$
(40)

The term $\mu_{t}(\boldsymbol{s}_{t})-\mu_{t}(\cdot)$ is nonnegative by
construction such that

$\displaystyle
g_{t}^{*}=\sup_{\boldsymbol{x}\in\mathcal{X}}g_{t}(\boldsymbol{x})\geq
f^{*}_{t}-f_{t}(\boldsymbol{s}_{t})$ $\displaystyle
f^{*}_{t}=\sup_{\boldsymbol{x}\in\mathcal{X}}f_{t}(\boldsymbol{x})$ (41)

and, therefore,

$\operatorname{\mathbb{P}}\left(g_{t}^{*}\geq\epsilon\right)\geq\operatorname{\mathbb{P}}\left(f^{*}_{t}-f_{t}(\boldsymbol{s}_{t})\geq\epsilon\right).$
(42)

We would now like to use the Borell-TIS inequality Borell (1975); Tsirelson et
al. (1976) to show that: if
$\epsilon>\operatorname{\mathbb{E}}\left(g_{t}^{*}\right)$, then

$\displaystyle\operatorname{\mathbb{P}}\left(g^{*}_{t}\geq\epsilon\right)\leq\exp\left(-\frac{1}{2}\left[\frac{\epsilon-\operatorname{\mathbb{E}}\left(g^{*}_{t}\right)}{2\sigma_{t}}\right]^{2}\right),$
(43)

where
$\sigma_{t}=\max_{\boldsymbol{x}\in\mathcal{X}}\sqrt{k_{t}(\boldsymbol{x},\boldsymbol{x})}$
and $2\sigma_{t}$ appears in the denominator (rather than $\sigma_{t}$)
because
$\max_{\boldsymbol{x}\in\mathcal{X}}c_{t}(\boldsymbol{x},\boldsymbol{x})\leq
4\sigma^{2}_{t}$. Since (43) is an increasing, continuous function of both
$\sigma_{t}\geq 0$ and
$0\leq\operatorname{\mathbb{E}}\left(g^{*}_{t}\right)<\epsilon$, the claim
will hold if these quantities vanish as the (global) fill distance
$h_{t}=\max_{\boldsymbol{x}\in\mathcal{X}}\min_{1\leq i\leq
t}\left\lVert\boldsymbol{x}-\boldsymbol{x}_{i}\right\rVert_{\infty}$ goes to
zero.

The former result is an immediate consequence of Lemma 5. Regarding the
latter, $f_{t}(\boldsymbol{s}_{t})-\mu_{t}(\boldsymbol{s}_{t})$ is a centered
random variable. It follows by linearity of expectation that

$\operatorname{\mathbb{E}}\left(g_{t}^{*}\right)=\operatorname{\mathbb{E}}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f_{t}(\boldsymbol{x})-\mu_{t}(\boldsymbol{x})\right].$
(44)

Next, denote the canonical pseudo-metric at time $t$ by

$d_{k_{t}}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\operatorname{\mathbb{E}}\left[(f_{t}(\boldsymbol{x})-f_{t}(\boldsymbol{x}^{\prime}))^{2}\right]^{\nicefrac{{1}}{{2}}}=\sqrt{k_{t}(\boldsymbol{x},\boldsymbol{x})-2k_{t}(\boldsymbol{x},\boldsymbol{x}^{\prime})+k_{t}(\boldsymbol{x}^{\prime},\boldsymbol{x}^{\prime})}.$
(45)

This pseudo-metric is non-increasing in $t$. To see this, let
$\beta=k(\boldsymbol{x}_{t+1},\boldsymbol{x})-k(\boldsymbol{x}_{t+1},\boldsymbol{x}^{\prime})$
and write

$\displaystyle\begin{split}d_{k_{t+1}}(\boldsymbol{x},\boldsymbol{x}^{\prime})^{2}&=k_{t+1}(\boldsymbol{x},\boldsymbol{x})-2k_{t+1}(\boldsymbol{x},\boldsymbol{x}^{\prime})+k_{t+1}(\boldsymbol{x}^{\prime},\boldsymbol{x}^{\prime})\\\
&=\underbracket{k_{t}(\boldsymbol{x},\boldsymbol{x})-2k_{t}(\boldsymbol{x},\boldsymbol{x}^{\prime})+k_{t}(\boldsymbol{x}^{\prime},\boldsymbol{x}^{\prime})}_{d_{k_{t}}(\boldsymbol{x},\boldsymbol{x}^{\prime})^{2}}-\underbracket{\beta^{2}\left[k_{t}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+1})+\gamma^{2}\right]^{-1}}_{\geq
0}.\end{split}$ (46)

As $t$ increases, points therefore become closer together under the
$d_{k_{t}}$ pseudo-metric. For this reason, the posterior
$\varepsilon$-packing number $N(\mathcal{X},d_{k_{t}},\varepsilon)$ is less-
equal to the prior $\varepsilon$-packing number
$N(\mathcal{X},d_{k},\varepsilon)$. Applying Lemma 7, we now have

$\operatorname{\mathbb{E}}\left(g_{t}^{*}\right)\leq\int_{0}^{\sigma_{t}}\sqrt{\log
N(\mathcal{X},d_{k_{t}},\varepsilon)}d\varepsilon\leq\int_{0}^{\sigma_{t}}\sqrt{\log
N(\mathcal{X},d_{k},\varepsilon)}d\varepsilon\leq
12\sigma_{t}\sqrt{2D+D\log\left(1+4L_{k}\sigma_{t}^{-2}\right)}.$ (47)

From here, note that (47) is an increasing, continuous function of
$\sigma_{t}$ that vanishes as $\sigma_{t}\to\infty$. By Lemma 5, the same is
true of $\sigma_{t}$ as a function of $h_{t}$. As a result, (43) becomes
arbitrarily small as $h_{t}\to 0$ and there exists a constant $h_{*}>0$ such
that this upper bound is less-equal to $\delta$ whenever $h_{t}\leq h_{*}$.
Finally, since $(\boldsymbol{x}_{t})$ is almost surely dense in $\mathcal{X}$,
there almost surely exists a time $T\in\mathbb{N}$ such that

$t\geq T\implies h_{t}\leq
h_{*}\implies\operatorname{\mathbb{P}}\left(g_{t}^{*}>\epsilon\right)\leq\delta\implies\operatorname{\mathbb{P}}\left(f_{t}^{*}-f_{t}(\boldsymbol{s}_{t})>\epsilon\right)\leq\delta.$
(48)

∎

######  Corollary 8.

Suppose assumptions A1–A3 hold and that there exists a constant
$\epsilon^{\prime}>0$ so that, with probability one,

$\displaystyle\limsup_{t\to\infty}\left[\max_{\boldsymbol{x}\in\mathcal{X}}\mu_{t}\left(\boldsymbol{x}\right)-\mu_{t}(\boldsymbol{s}_{t})\right]\leq\epsilon^{\prime},$
(49)

where
$\boldsymbol{s}_{t}\in\operatorname*{\arg\max}_{\boldsymbol{x}\in\mathbf{X}_{t}}\mu_{t}(\boldsymbol{x})$.
Then, for every $\epsilon>\epsilon^{\prime}$ and $\delta\in(0,1]$, there
almost surely exists a time $T\in\mathbb{N}$ such that, for all $t\geq T$,

$\displaystyle\operatorname{\mathbb{P}}\left[\sup_{\boldsymbol{x}\in\mathcal{X}}f_{t}(\boldsymbol{x})-f_{t}(\boldsymbol{s}_{t})\leq\epsilon\right]\geq
1-\delta.$ (50)

###### Proof.

Per (49), there almost surely exists an $S\in\mathbb{N}$ such that $t\geq
S\implies\max_{\boldsymbol{x}\in\mathcal{X}}\mu_{t}\left(\boldsymbol{x}\right)-\mu_{t}(\boldsymbol{s}_{t})\leq\epsilon^{\prime}$.
Theorem 2 therefore implies there almost surely exists a $T\geq S$ so that

$t\geq
T\implies\operatorname{\mathbb{P}}\big{[}f^{*}_{t}-f_{t}(\boldsymbol{s}_{t})\geq\epsilon-\epsilon^{\prime}\big{]}\leq\delta,$
(51)

which completes the proof. ∎

The assumption that the posterior mean approaches its maximum on
$(\boldsymbol{x}_{t})$ protects against adversarial cases where—no matter how
densely we observe $f$—there is always an
$\boldsymbol{x}\in\mathcal{X}\setminus\mathbf{X}_{t}$ whose expected value
$\mu_{t}\left(\boldsymbol{x}\right)$ exceeds
$\mu_{t}\left(\boldsymbol{s}_{t}\right)$ by at least $\epsilon$. Note that
(49) becomes a necessary condition when $\delta<\frac{1}{2}$. Nevertheless, it
is unclear how to ensure this condition without making stronger assumptions
for $f$ and $(\boldsymbol{x}_{t})$. One can use A2 and the Cauchy-Schwarz
inequality to show that the posterior mean is Lipschitz continuous Lederer et
al. (2019b); but, its Lipschitz constant may continue to grow as $t\to\infty$,
so (49) may not hold.

See 3

###### Proof.

Under assumptions A1–A3, Theorem 2 implies that there almost surely exists an
$S\in\mathbb{N}$ so that, at all times $t\geq S$, $\boldsymbol{s}_{t}$ is
$(\epsilon,\delta_{\mathrm{mod}})$-optimal. Since this event occurs infinitely
often with probability one and each estimate
$\Psi_{t}^{n}(\boldsymbol{s}_{t})$ is independent and unbiased, there almost
surely exists a least $T\in\mathbb{N}$ such that
$\Psi_{T}(\boldsymbol{s}_{T})\geq 1-\delta_{\mathrm{mod}}$ and
$\Psi_{T}^{n}(\boldsymbol{s}_{T})\geq 1-\delta_{\mathrm{mod}}$.

By assumption, each event

$\displaystyle\mathcal{A}_{t}=\left\\{\left|\Psi_{t}^{n}(\boldsymbol{s}_{t})-\Psi_{t}(\boldsymbol{s}_{t})\right|\leq\left|\Psi_{t}^{n}(\boldsymbol{s}_{t})-\lambda\right|\right\\}$
$\displaystyle t=1,2,\ldots$ (52)

occurs with probability at least $1-\delta_{\mathrm{est}}^{t}$. If
$\mathcal{A}_{t}$ holds, then $\Psi_{t}^{n}(\boldsymbol{s}_{t})\geq
1-\delta_{\mathrm{mod}}$ implies $\Psi_{t}(\boldsymbol{s}_{t})\geq
1-\delta_{\mathrm{mod}}$. Per the union bound, the probability of terminating
prior to $T$ is therefore no larger than
$\sum_{t=1}^{T-1}\delta_{\mathrm{est}}^{t}\leq\delta_{\mathrm{est}}$.

If the algorithm reaches $T$, then the probability that $\boldsymbol{s}_{T}$
fails to be $\epsilon$-optimal is less-equal to $\delta_{\mathrm{mod}}$ by
construction. Hence, the algorithm almost surely terminates and returns an
$\epsilon$-optimal solution with probability at least
$1-\delta_{\mathrm{mod}}-\delta_{\mathrm{est}}\geq 1-\delta$. ∎

##  Appendix B Experiment Details

#### Implementation

Code is available online at https://github.com/j-wilson/trieste_stopping.
Experiments were run using a combination of GPFlow Matthews et al. (2017) and
Trieste Picheny et al. (2023). Runtimes reported in Figure 3 were measured on
an Apple M1 Pro Chip using an off-the-shelf version of TensorFlow Abadi et al.
(2015).

#### Model specification

We employed Gaussian process priors $f\sim\mathcal{GP}(\mu,k)$ with constant
mean functions $\mu(\cdot)=c$ and Matérn-$\nicefrac{{5}}{{2}}$ covariance
functions equipped with ARD lengthscales.

When optimizing functions drawn from GP priors, we set the prior mean to zero
and used unit variance kernels with lengthscales
$\ell_{i}=\frac{1}{4}\sqrt{D}$. Noise variances are indicated alongside
reported results.

When optimizing black-box functions, we employed broad and uninformative
hyperpriors. Let $[\mathcal{X}]_{i}=[a_{i},b_{i}]$ be the range of the $i$-th
design variable, $q_{t}:[0,1]\to\mathbb{R}$ be the empirical quantile function
of $y$ at time $t$, and
$\nu_{t}=\overline{\operatorname{Var}}\left[\boldsymbol{y}_{t-1}\right]$ be
the empirical variance of observations
$\boldsymbol{y}_{t-1}=\left\\{y(\boldsymbol{x}_{1}),\ldots,y(\boldsymbol{x}_{t-1})\right\\}$.
Our hyperpriors are then as follows:

Name | Distribution |  Parameters
---|---|---
Constant Mean | $\operatorname{Uniform}(a,b)$ | $a=q_{t}(0.05)$ |  $b=q_{t}(0.95)$
Log Kernel Variance | $\operatorname{Uniform}\left(a,b\right)$ | $a=\log(10^{-1}\nu_{t})$ |  $b=\log(10\nu_{t})$
Log Noise Variance | $\operatorname{Uniform}\left(a,b\right)$ | $a=\log(10^{-9}\nu_{t})$ |  $b=\log(10\nu_{t})$
$i$-th Lengthscale | $\operatorname{LogNormal}\left(\mu,\sigma\right)$ | $\mu=\frac{1}{2}(b_{i}-a_{i})$ |  $\sigma=1$

Note that we directly parameterize certain hyperparamters in log-space and
that, e.g., $\log(\theta)\sim\operatorname{Uniform}(a,b)$ is not the same as
$\theta\sim\operatorname{LogUniform}(e^{a},e^{b})$ Bishop (2006).

#### Acquisition function

In our experiments, we defined the set of feasible solutions at time $t$ as
the set of evaluated points $\mathbf{X}_{t}$. Under these circumstances, one
can show that the optimal one-step policy is given by an “in-sample” version
of the Knowledge Gradient strategy Močkus (1975); Frazier et al. (2009). Let

$\mu_{t+1}\left(\cdot;\boldsymbol{x},z\right)=\mu_{t}(\cdot)+\frac{k_{t}\left(\cdot,\boldsymbol{x}\right)z}{\sqrt{k_{t}(\boldsymbol{x},\boldsymbol{x})+\gamma^{2}}}$
(53)

be the posterior mean of $f$ at time $t+1$ if we observe
$y_{t+1}=\mu_{t}(\boldsymbol{x})+\sqrt{k_{t}(\boldsymbol{x},\boldsymbol{x})+\gamma^{2}}z$,
where $z\sim\mathcal{N}(0,1)$. Then, the aforementioned acquisition function
is given by

$\displaystyle\operatorname{ISKG}_{t}(\boldsymbol{x})$
$\displaystyle=\operatorname{\mathbb{E}}_{z}\Big{[}\max\mu_{t+1}\left(\mathbf{X}_{t}\cup\left\\{\boldsymbol{x}\right\\};\boldsymbol{x},z\right)\Big{]}-\max\mu_{t}(\mathbf{X}_{t})$
$\displaystyle z\sim\mathcal{N}(0,1).$ (54)

$\operatorname{ISKG}$ is identical to the Expected Improvement function when
$\gamma^{2}=0$ Scott et al. (2011), but avoids pathologies (such as re-
evaluated previously observed points) when $\gamma^{2}>0$.

In practice, we estimated (54) with Gauss-Hermite quadrature and maximized it
using multi-start gradient ascent Wilson et al. (2018); Balandat et al.
(2020). Starting positions we obtained by running CMA-ES Hansen (2016) several
times to partial convergence. The best point from each run was then combined
with a large number of random points and the top $16$ points were fine-tuned
using L-BFGS-B Byrd et al. (1995).

##  Appendix C Extended Results

This section provides extended results for experiments presented in the body.
Results are given for each problem in the order they appeared in Tables 1 and
1. In each of the following table, we report success rates along with medians
and interquartile ranges of stopping times, log regrets, and log excess
regrets. Here, excess regret is defined as the amount by which regrets
exceeded $\epsilon$ for runs that failed to return $\epsilon$-optimal
solutions, i.e.
$\sup_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})-f(\boldsymbol{s}_{T})+\epsilon$.
For compactness, table numbers and problems names are reported in the upper
left corner.

2\. GP-2∗, $10^{-6}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  8 |  9 |  12 |  -5.88 |  -5.34 |  -4.81 |  – |  – | –
Budget∗ | 0.96 |  18 |  18 |  18 |  -4.25 |  -3.16 |  -2.21 |  -1.32 |  -0.98 | -0.92
Budget | 0.98 |  33 |  33 |  33 |  -4.99 |  -4.30 |  -3.78 |  -1.01 |  -0.99 | -0.97
Acq. Cutoff∗ | 0.99 |  23 |  26 |  33 |  -4.78 |  -4.25 |  -3.81 |  -1.03 |  -1.03 | -1.03
UCB – LCB | 0.89 |  10 |  12 |  15 |  -2.54 |  -2.00 |  -1.43 |  -1.00 |  -0.86 | -0.21
PRB (ours) | 0.97 |  14 |  17 |  21 |  -3.72 |  -2.96 |  -2.36 |  -1.30 |  -1.03 | -0.94
3\. GP-2∗, $10^{-2}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  8 |  10 |  12 |  -4.35 |  -3.89 |  -3.20 |  – |  – | –
Budget∗ | 0.96 |  23 |  23 |  23 |  -2.49 |  -1.83 |  -1.42 |  -1.84 |  -1.66 | -1.41
Budget | 1.00 |  33 |  33 |  33 |  -2.91 |  -2.08 |  -1.60 |  – |  – | –
Acq. Cutoff∗ | 1.00 |  34.75 |  60.5 |  90.5 |  -3.08 |  -2.50 |  -2.20 |  – |  – | –
UCB – LCB | 1.00 |  70 |  82 |  99 |  -3.33 |  -2.77 |  -2.25 |  – |  – | –
PRB (ours) | 0.99 |  17 |  23 |  32 |  -2.51 |  -1.87 |  -1.58 |  -0.11 |  -0.11 | -0.11
4\. GP-4∗, $10^{-6}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  17.75 |  26 |  34 |  -4.24 |  -3.73 |  -3.33 |  – |  – | –
Budget∗ | 0.95 |  65 |  65 |  65 |  -3.56 |  -3.11 |  -2.67 |  -1.70 |  -1.18 | -0.62
Budget | 0.95 |  65 |  65 |  65 |  -3.56 |  -3.11 |  -2.67 |  -1.70 |  -1.18 | -0.62
Acq. Cutoff∗ | 1.00 |  59.5 |  76.5 |  90 |  -3.66 |  -3.22 |  -2.93 |  – |  – | –
UCB – LCB | 0.66 |  17.75 |  22.5 |  30 |  -1.98 |  -1.49 |  -0.81 |  -1.26 |  -0.70 | -0.45
PRB (ours) | 0.99 |  46 |  64 |  78 |  -3.49 |  -3.10 |  -2.70 |  -2.67 |  -2.67 | -2.67
5\. GP-4∗, $10^{-2}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  19 |  28.5 |  42 |  -2.66 |  -2.36 |  -1.92 |  – |  – | –
Budget∗ | 0.95 |  95 |  95 |  95 |  -2.35 |  -1.82 |  -1.52 |  -1.63 |  -1.45 | -1.30
Budget | 0.85 |  65 |  65 |  65 |  -1.94 |  -1.57 |  -1.22 |  -1.73 |  -1.35 | -1.02
Acq. Cutoff∗ | 0.96 |  69.75 |  89 |  111 |  -2.24 |  -1.84 |  -1.53 |  -1.75 |  -1.47 | -1.23
UCB – LCB | 1.00 |  223.75 |  256 |  256 |  -2.44 |  -1.97 |  -1.73 |  – |  – | –
PRB (ours) | 0.96 |  67 |  86.5 |  113.25 |  -2.21 |  -1.81 |  -1.50 |  -1.92 |  -1.49 | -1.14
6\. GP-4, $10^{-6}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  21.5 |  34 |  50.5 |  -4.17 |  -3.66 |  -3.15 |  – |  – | –
Budget∗ | 0.95 |  80 |  80 |  80 |  -3.59 |  -3.13 |  -2.58 |  -1.70 |  -1.47 | -0.88
Budget | 0.85 |  65 |  65 |  65 |  -3.35 |  -2.88 |  -1.88 |  -1.11 |  -0.79 | -0.55
Acq. Cutoff∗ | 0.97 |  61.75 |  77.5 |  91 |  -3.49 |  -3.15 |  -2.77 |  -1.58 |  -1.47 | -1.43
UCB – LCB | 0.30 |  15 |  18 |  21.25 |  -1.25 |  -0.34 |  -0.08 |  -0.53 |  -0.25 | -0.04
PRB (ours) | 0.88 |  48.25 |  61 |  75.25 |  -3.27 |  -2.85 |  -2.24 |  -1.41 |  -0.66 | -0.49
7\. GP-4, $10^{-2}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  29.75 |  50 |  74.75 |  -2.67 |  -2.32 |  -1.93 |  – |  – | –
Budget∗ | 0.95 |  158 |  158 |  158 |  -2.24 |  -1.79 |  -1.51 |  -2.83 |  -1.39 | -0.63
Budget | 0.61 |  65 |  65 |  65 |  -1.72 |  -1.30 |  -0.69 |  -1.26 |  -0.61 | -0.31
Acq. Cutoff∗ | 0.96 |  88.75 |  112 |  142 |  -2.25 |  -1.75 |  -1.46 |  -2.08 |  -1.75 | -1.68
UCB – LCB | 0.80 |  121.5 |  224 |  256 |  -2.30 |  -1.76 |  -1.19 |  -0.84 |  -0.31 | -0.01
PRB (ours) | 0.92 |  73.75 |  100.5 |  129.5 |  -2.07 |  -1.69 |  -1.39 |  -1.78 |  -1.55 | -1.31
8\. GP-6∗, $10^{-6}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 0.99 |  25 |  39.5 |  77.25 |  -3.71 |  -3.15 |  -2.69 |  -2.74 |  -2.74 | -2.74
Budget∗ | 0.95 |  125 |  125 |  125 |  -3.04 |  -2.69 |  -2.26 |  -1.54 |  -1.16 | -0.84
Budget | 0.86 |  97 |  97 |  97 |  -2.93 |  -2.42 |  -1.70 |  -1.51 |  -1.15 | -0.74
Acq. Cutoff∗ | 0.97 |  86 |  115.5 |  153.5 |  -2.98 |  -2.60 |  -2.31 |  -2.05 |  -1.67 | -1.26
UCB – LCB | 0.50 |  24 |  32.5 |  41.25 |  -1.63 |  -1.03 |  -0.47 |  -1.04 |  -0.62 | -0.20
PRB (ours) | 0.98 |  97 |  133.5 |  178 |  -3.15 |  -2.77 |  -2.47 |  -2.47 |  -2.21 | -1.94
9\. GP-6∗, $10^{-2}$ |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  33 |  64 |  96.5 |  -2.41 |  -2.00 |  -1.72 |  – |  – | –
Budget∗ | 0.96 |  228 |  228 |  228 |  -1.97 |  -1.70 |  -1.37 |  -1.73 |  -1.60 | -1.39
Budget | 0.70 |  97 |  97 |  97 |  -1.53 |  -1.25 |  -0.87 |  -1.34 |  -0.98 | -0.71
Acq. Cutoff∗ | 0.92 |  110.5 |  145 |  187.5 |  -1.76 |  -1.53 |  -1.30 |  -1.87 |  -1.56 | -1.33
UCB – LCB | 1.00 |  499 |  512 |  512 |  -2.25 |  -1.91 |  -1.60 |  – |  – | –
PRB (ours) | 1.00 |  170 |  235 |  295 |  -2.02 |  -1.76 |  -1.53 |  – |  – | –
10\. Branin-2 |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  13 |  18 |  21.25 |  -6.41 |  -6.34 |  -6.17 |  – |  – | –
Budget∗ | 0.95 |  26 |  26 |  26 |  -2.79 |  -1.97 |  -1.45 |  -1.62 |  -1.54 | -1.15
Budget | 1.00 |  33 |  33 |  33 |  -3.44 |  -3.02 |  -2.61 |  – |  – | –
Acq. Cutoff∗ | 1.00 |  64 |  64 |  64 |  -6.23 |  -5.99 |  -5.69 |  – |  – | –
UCB – LCB | 0.99 |  26.75 |  31 |  32 |  -3.08 |  -2.75 |  -2.38 |  0.16 |  0.16 | 0.16
PRB (ours) | 0.99 |  31 |  33 |  35 |  -3.44 |  -3.07 |  -2.72 |  0.16 |  0.16 | 0.16
11\. Hartmann-3 |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  11 |  13 |  15.25 |  -6.66 |  -6.64 |  -6.58 |  – |  – | –
Budget∗ | 0.96 |  23 |  23 |  23 |  -4.39 |  -3.81 |  -3.32 |  -0.80 |  -0.37 | -0.18
Budget | 1.00 |  49 |  49 |  49 |  -6.58 |  -6.39 |  -6.05 |  – |  – | –
Acq. Cutoff∗ | 1.00 |  23 |  25 |  27.25 |  -4.89 |  -4.50 |  -4.02 |  – |  – | –
UCB – LCB | 0.83 |  14 |  15 |  16 |  -2.20 |  -1.98 |  -1.61 |  -0.17 |  -0.17 | -0.17
PRB (ours) | 1.00 |  17 |  19 |  21 |  -3.81 |  -3.41 |  -2.90 |  – |  – | –

12\. Hartmann-6 |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 0.67 |  25.75 |  35 |  256 |  -5.70 |  -5.69 |  -0.92 |  -1.72 |  -1.72 | -1.72
Budget∗ | 0.67 |  256 |  256 |  256 |  -5.70 |  -5.69 |  -0.92 |  -1.72 |  -1.72 | -1.72
Budget | 0.67 |  97 |  97 |  97 |  -5.66 |  -5.59 |  -0.92 |  -1.72 |  -1.72 | -1.72
Acq. Cutoff∗ | 0.67 |  35 |  38 |  42 |  -3.28 |  -2.81 |  -0.92 |  -1.71 |  -1.70 | -1.68
UCB – LCB | 0.46 |  18 |  26 |  31 |  -1.75 |  -0.85 |  0.03 |  -1.09 |  -0.59 | 0.47
PRB (ours) | 0.64 |  37 |  40 |  43 |  -3.29 |  -2.82 |  -0.92 |  -1.71 |  -1.70 | -1.67
13\. Rosenbrock-4 |  $\boldsymbol{\epsilon}$-optimal | Stopping Time | Log10 Regret |  Log10 Excess Regret
---|---|---|---|---
Method | % |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 |  Q3 |  Q1 |  Q2 | Q3
Oracle∗ | 1.00 |  26 |  33 |  41 |  – |  – |  – |  – |  – | –
Budget∗ | 0.95 |  47 |  47 |  47 |  -4.89 |  -4.53 |  -4.27 |  -5.35 |  -5.23 | -4.60
Budget | 1.00 |  65 |  65 |  65 |  -6.06 |  -5.68 |  -5.44 |  – |  – | –
Acq. Cutoff∗ | 0.94 |  60 |  62 |  64 |  -6.06 |  -5.59 |  -4.88 |  -4.40 |  -4.33 | -4.26
UCB – LCB | 1.00 |  71 |  73.5 |  76 |  – |  – |  – |  – |  – | –
PRB (ours) | 1.00 |  81 |  84 |  88 |  – |  – |  – |  – |  – | –


This article has been removed by arXiv administrators because the submitter
did not have the right to agree to the license at the time of submission.
[2020-11-19]