Bayesian Optimization Meets Bayesian Optimal Stopping

Bayesian Optimization Meets Bayesian Optimal Stopping
Zhongxiang Dai 1 Haibin Yu 1 Bryan Kian Hsiang Low 1 Patrick Jaillet 2
Abstract Silver et al., 2016). However, a major difficulty faced by ML

Bayesian optimization (BO) is a popular paradigm practitioners is the choice of model hyperparameters which
for optimizing the hyperparameters of machine significantly impacts the predictive performance. This calls
learning (ML) models due to its sample efficiency. for the need to develop hyperparameter optimization algo-
Many ML models require running an iterative rithms that have to be sample-efficient since the training
training procedure (e.g., stochastic gradient de- of many modern ML models consumes massive computa-
scent). This motivates the question whether in- tional resources. To this end, Bayesian optimization (BO)
formation available during the training process is a popular paradigm due to its high sample efficiency and
(e.g., validation accuracy after each epoch) can be strong theoretical performance guarantee (Shahriari et al.,
exploited for improving the epoch efficiency of 2016). In particular, the BO algorithm based on the Gaus-
BO algorithms by early-stopping model training sian process-upper confidence bound (GP-UCB) acquisition
under hyperparameter settings that will end up function has been shown to achieve no regret asymptotically
under-performing and hence eliminating unneces- and perform competitively in practice (Srinivas et al., 2010).
sary training epochs. This paper proposes to unify Many ML models require running an iterative training proce-
BO (specifically, Gaussian process-upper confi- dure for some number of epochs such as stochastic gradient
dence bound (GP-UCB)) with Bayesian optimal descent for neural networks (LeCun et al., 2015) and boost-
stopping (BO-BOS) to boost the epoch efficiency ing procedure for gradient boosting machines (Friedman,
of BO. To achieve this, while GP-UCB is sample- 2001). During BO, any query of a hyperparameter setting
efficient in the number of function evaluations, usually involves training the ML model for a fixed number of
BOS complements it with epoch efficiency for epochs. Information typically available during the training
each function evaluation by providing a principled process (e.g., validation accuracy after each epoch) is rarely
optimal stopping mechanism for early stopping. exploited for improving the epoch efficiency of BO algo-
BO-BOS preserves the (asymptotic) no-regret per- rithms, specifically, by early-stopping model training under
formance of GP-UCB using our specified choice hyperparameter settings that will end up under-performing,
of BOS parameters that is amenable to an ele- hence eliminating unnecessary training epochs. Note that
gant interpretation in terms of the exploration- this objective is different from that of standard early stop-
exploitation trade-off. We empirically evaluate ping during the training of neural networks, which is used
the performance of BO-BOS and demonstrate its to prevent overfitting.
generality in hyperparameter optimization of ML
models and two other interesting applications. To address this challenging issue, a number of works have
been proposed to make BO more epoch-efficient: Freeze-
thaw BO (Swersky et al., 2014) explores a diverse collection
1. Introduction
of hyperparameter settings in the initial stage by training
The state-of-the-art machine learning (ML) models have their ML models with a small number of epochs, and then
recently reached an unprecedented level of predictive per- gradually focuses on (exploits) a small number of promising
formance in several applications such as image recognition, settings. Despite its promising epoch efficiency, its perfor-
complex board games, among others (LeCun et al., 2015; mance is not theoretically guaranteed and its computational
1
cost can be excessive. Multi-fidelity BO1 (Kandasamy et al.,
Department of Computer Science, National University of 2016; 2017), as well as its extensions using predictive en-
Singapore, Republic of Singapore. 2 Department of Electrical
Engineering and Computer Science, Massachusetts Institute of tropy search (Zhang et al., 2017) and value of information
Technology, USA. Correspondence to: Bryan Kian Hsiang Low for fidelity selection (Wu & Frazier, 2017), reduces the
<lowkh@comp.nus.edu.sg>. resource consumption of BO by utilizing low-fidelity func-
1
Proceedings of the 36 th International Conference on Machine A closely related counterpart is multi-fidelity active learn-
Learning, Long Beach, California, PMLR 97, 2019. Copyright ing (Zhang et al., 2016).
2019 by the author(s).
tions which can be obtained by using a subset of the training objectives of GP-UCB vs. BOS (respectively, objective func-
data or by training the ML model for a small number of tion v.s. expected loss), BO-BOS can preserve the trademark
epochs. However, in each BO iteration, since the fidelity (asymptotic) no-regret performance of GP-UCB with our
(e.g., number of epochs) is determined before function eval- specified choice of BOS parameters that is amenable to an el-
uation, it is not influenced by information that is typically egant interpretation in terms of the exploration-exploitation
available during the training process (e.g., validation accu- trade-off (Section 4). Though the focus of our work here is
racy after each epoch). In addition to BO, attempts have on epoch-efficient BO for hyperparameter tuning, we addi-
also been made to improve the epoch efficiency of other tionally evaluate the performance of BO-BOS empirically
hyperparameter optimization algorithms: Some heuristic in two other interesting applications to demonstrate its gen-
methods (Baker et al., 2017; Domhan et al., 2015; Klein erality: policy search for reinforcement learning, and joint
et al., 2017) predict the final training outcome based on hyperparameter tuning and feature selection (Section 5).
partially trained learning curves in order to identify hyper-
parameter settings that are predicted to under-perform and 2. Background and Notations
early-stop their model training. Hyperband (Li et al., 2017),
which dynamically allocates the computational resource 2.1. Bayesian Optimization (BO) and GP-UCB
(e.g., training epochs) through random sampling and elimi- Consider the problem of sequentially maximizing an un-
nates under-performing hyperparameter settings by succes- known objective function f : D → R representing the val-
sive halving, has been proposed and shown to perform well idation accuracy over a compact input domain D ⊆ Rd of
in practice. Both the learning curve prediction methods and different hyperparameter settings for training an ML model:
Hyperband can be combined with BO to further improve In each iteration t = 1, . . . , T , an input query zt , [xt , nt ]
the epoch efficiency (Domhan et al., 2015; Falkner et al., of a hyperparameter setting (comprising N0 < nt ≤ N
2018; Klein et al., 2017), but their resulting performances training epochs and a vector xt of the other hyperparam-
are not theoretically guaranteed. Despite these recent ad- eters) is selected for evaluating the validation accuracy f
vances, we still lack an epoch-efficient algorithm that can of the ML model to yield a noisy observed output (valida-
incorporate early stopping into BO (i.e., by exploiting infor- tion accuracy) yt , f (zt ) + with i.i.d. Gaussian noise
mation available during the training process) and yet offer a ∼ N (0, σ 2 ) and noise variance σ 2 . Since every eval-
theoretical performance guarantee, the design of which is uation of f is costly (Section 1), our goal is to strategi-
likely to require a principled decision-making mechanism cally select input queries to approach the global maximizer
for determining the optimal stopping time. z∗ , arg maxz∈D f (z) as rapidly as possible. This can be
achieved by minimizing a standard BO objective such as
Optimal stopping is a classic research topic in statistics and
the simple regret ST , f (z∗ ) − maxt∈{1,...,T } f (zt ). A
operations research regarding sequential decision-making
BO algorithm is said to guarantee no regret asymptotically
problems whose objective is to make the optimal stopping
if it satisfies limT →∞ ST = 0, thus implying that it will
decision with a small number of observations (Ferguson,
eventually converge to the global maximum.
2006). In Bayesian optimal stopping (BOS) or Bayesian
sequential design, the decision between stopping vs. continu- To guarantee no regret, the belief of f is modeled by a
ing is made to maximize the expected utility or, equivalently, Gaussian process (GP). Let {f (z)}z∈D denote a GP, that
minimize the expected loss (Powell & Ryzhov, 2012). BOS is, every finite subset of {f (z)}z∈D follows a multivariate
has found success in application domains such as finance Gaussian distribution (Rasmussen & Williams, 2006). Then,
(Longstaff & Schwartz, 2001), clinical design (Brockwell & a GP is fully specified by its prior mean µ(z) and covariance
Kadane, 2003; Müller et al., 2007; Wathen & Thall, 2008), k(z, z0 ) for all z, z0 ∈ D, which are, respectively, assumed
and economics (Davis & Cairns, 2012). The capability of w.l.o.g. to be µ(z) = 0 and k(z, z0 ) ≤ 1 for notational
BOS in providing a principled optimal stopping mechanism simplicity. Given a column vector yT , [yt ]> t=1,...,T of
makes it a prime candidate for introducing early stopping noisy outputs observed from evaluating f at the selected
into BO in a theoretically sound and rigorous way. input queries z1 , . . . , zT ∈ D after T iterations, the GP
posterior belief of f at some input z ∈ D is a Gaussian with
This paper proposes to unify Bayesian optimization (specifi-
the following posterior mean µT (z) and variance σT2 (z):
cally, GP-UCB) with Bayesian optimal stopping (BO-BOS)
to boost the epoch efficiency of BO (Section 3). Intu- µT (z) , kT (z)> (KT + σ 2 I)−1 yT ,
itively, GP-UCB is acclaimed for being sample-efficient (1)
σT2 (z) , k(z, z) − kT (z)> (KT + σ 2 I)−1 kT (z)
in the number of function evaluations while BOS can reduce
the required number of epochs for each function evalua- where KT , [k(zt , zt0 )]t,t0 =1,...,T and kT (z) ,
tion. BO-BOS unifies the best of both worlds to yield an [k(zt , z)]>
t=1,...,T . In each iteration t of BO, an input query
epoch-efficient hyperparameter optimization algorithm. In- zt ∈ D is selected to maximize the GP-UCB acquisition
terestingly, in spite of the seemingly disparate optimization function (Srinivas et al., 2010) that trades off between ob-
serving an expected maximum (i.e., with a large GP poste- Solving the BOS problem (2) yields a Bayes-optimal de-
rior mean µt−1 (z)) given the current GP belief of f (i.e., cision rule in each epoch n: Take the Bayes-optimal stop-
exploitation) vs. that of high predictive uncertainty (i.e., ping decision if the expected loss of either terminal deci-
2
with a large GP posterior variance σt−1 (z)) to improve sion d1 or d2 is at most that of the continuation decision
the GP belief of f over D (i.e.,√ exploration). That is, d0 , that is, min{ Eθt |yt,n [l(d1 , θt )], Eθt |yt,n [l(d2 , θt )] } ≤
zt , arg maxz∈D µt−1 (z) + βt σt−1 (z) where the pa- cd0 + Eyt,n+1 |yt,n [ρt,n+1 (yt,n+1 )]. Otherwise, continue
rameter βt > 0 is set to trade off between exploitation vs. model training to yield the noisy observed output (validation
exploration for guaranteeing no regret asymptotically with accuracy) yt,n+1 and repeat this rule in the next epoch n + 1.
high probability.
3. BO-BOS Algorithm
2.2. Bayesian Optimal Stopping (BOS)
In this section, we will describe our proposed BO-BOS algo-
BOS provides a principled mechanism for making the
rithm (Section 3.1) and define the loss function l in BOS (2)
Bayes-optimal stopping decision with a small number of
such that it can serve as an effective early-stopping mech-
observations. As shall be seen in Algorithm 1, in each itera-
anism in BO (Section 3.2). We focus on problem settings
tion t of BO-BOS, BOS is used to early-stop model training
where the objective function f is bounded and monotoni-
under the selected input hyperparameters xt that will end
cally increasing in n:
up under-performing, hence reducing the required number
of training epochs. In a BOS problem, the goal is to de- Assumption 1. (a) f (z) ∈ [0, 1] for all z ∈ D and (b)
cide whether to stop and conclude either hypothesis/event f ([x, n]) ≤ f ([x, n + 1]) for all x and n = 1, . . . , N − 1.
θt = θt,1 or θt = θt,2 corresponding to terminal decision d1 Assumption 1a is not restrictive since it applies to any
or d2 , or to gather one more observation via the continuation bounded f with a proper transformation. Assumption 1b
decision d0 . Let yt,n0 be the noisy output (validation accu- holds reasonably well in a number of important ML prob-
racy) observed in epoch n0 and yt,n , [yt,n0 ]> n0 =1,...,n be a lems: (a) f represents the validation accuracy of an ML
vector of noisy outputs observed up till epoch n in iteration model and n denotes the number of training epochs or the
t. Recall from Section 2.1 that in iteration t, the ML model number of selected features during feature selection, and (b)
is trained using the selected input hyperparameter setting f represents the (discounted) cumulative rewards (assuming
[xt , nt ] to yield the noisy observed output (validation accu- non-negative rewards) in reinforcement learning (RL) and n
racy) yt . So, yt = yt,nt for t = 1, . . . , T . After each epoch denotes the number of steps taken by the agent in the envi-
n, the posterior belief of event θt is updated to Pr(θt |yt,n ) ronment. Our experiments in Section 5 will demonstrate that
which will be used to compute the expected losses of termi- BO-BOS outperforms the state-of-the-art hyperparameter
nal decisions d1 and d2 . Such a loss function l has to encode optimization algorithms in these ML problems.
the cost of making a wrong decision. Define ρt,n (yt,n ) as
the minimum expected loss among all decisions in epoch n: 3.1. Algorithm Description
ρt,n (yt,n ) , min{ Eθt |yt,n [l(d1 , θt )], Eθt |yt,n [l(d2 , θt )],
In each iteration t of BO-BOS (Algorithm 1), the input
cd0 + Eyt,n+1 |yt,n [ρt,n+1 (yt,n+1 )] } hyperparameters xt are selected to maximize the GP-UCB
(2)
acquisition function with the input dimension of training
for n = N0 + 1, . . . , N − 1 where the first two terms are
epochs fixed at N (line 2). The ML model is trained using
the expected losses of terminal decisions d1 and d2 , the last
xt for N0 initial training epochs to yield the noisy observed
term sums the immediate cost cd0 and expected future loss
outputs (validation accuracies) yt,N0 (line 3). After that,
of making the continuation decision d0 to continue model
the BOS problem is solved (line 5) to obtain Bayes-optimal
training in the next epoch n + 1 to yield the noisy observed
decision rules (see Sections 2.2 and 3.2). Then, in each
output (validation accuracy) yt,n+1 , and ρt,N (yt,N ) ,
epoch n > N0 , model training continues under xt to yield
min{ Eθt |yt,N [l(d1 , θt )], Eθt |yt,N [l(d2 , θt )] }. Since ρt,n
the noisy observed output (validation accuracy) yt,n (line
depends on ρt,n+1 , it naturally prompts the use of backward
8). If both of the following conditions are satisfied (line 9):
induction to solve the BOS problem (2) exactly, which is
unfortunately intractable due to an uncountable set of possi- C1. BOS decision rule in epoch n outputs stopping decision;
ble observed outputs yt,n+1 . This computational difficulty
C2. σt−1 ([xt , N ]) ≤ κ σt−1 ([xt , n]) ,
can be overcome using approximate backward induction
techniques (Brockwell & Kadane, 2003; Müller et al., 2007) then model training is early-stopped in epoch nt = n. Other-
whose main ideas include using summary statistics to repre- wise, the above procedure is repeated in epoch n+1. If none
sent the posterior beliefs, discretizing the space of summary of the training epochs n = N0 + 1, . . . , N − 1 satisfy both
statistics, and approximating the expectation terms via sam- C1 and C2, then nt = N (i.e., no early stopping). Finally,
pling. Appendix A describes a commonly-used approximate the GP posterior belief (1) is updated with the selected input
backward induction algorithm (Müller et al., 2007). hyperparameter setting zt = [xt , nt ] and the corresponding
noisy observed output (validation accuracy) yt = yt,nt (line We define l(d1 , θt ) and l(d2 , θt ) of the respective terminal
11). BO-BOS then proceeds to the next iteration t + 1. decisions d1 and d2 as 0-K loss functions which are com-
monly used in clinical designs (Jiang et al., 2013; Lewis &
To understand the rationale of our choices of C1 and C2, the
Berry, 1994) due to their simplicity and interpretablility:
BOS decision rule in C1 recommends the stopping decision
to early-stop model training in epoch n if it concludes that l(d1 , θt ) , K1 1θt =θt,2 and l(d2 , θt ) , K2 1θt =θt,1 (3)
model training under xt for N epochs will produce a valida-
tion accuracy not exceeding the currently found maximum where the parameters K1 > 0 and K2 > 0 represent the
in iterations 1, . . . , t − 1; this will be formally described in costs of making the wrong terminal decisions d1 and d2 ,
Section 3.2. On the other hand, C2 prefers to evaluate the respectively. Since f ([xt , N ]) is not known in epoch n <
validation accuracy f of the ML model with the input query N , the expected losses of terminal decisions d1 and d2 have
[xt , n] of fewer training epochs n < N than [xt , N ] if the to be evaluated instead:
uncertainty of the validation accuracy f ([xt , N ]) achieved
by model training under xt for N epochs is not more than a Eθt |yt,n [l(d1 , θt )] = K1 Pr(θt = θt,2 |yt,n ) ,
(4)
factor of κ ≥ 1 of that of f ([xt , n]) for n epochs; the degree Eθt |yt,n [l(d2 , θt )] = K2 Pr(θt = θt,1 |yt,n ) .
of preference is controlled by parameter κ. Thus, by satis-
fying both C1 and C2, C2 lends confidence to the resulting According to (4), if K1 (K2 ) is set to +∞, then
performance of model training under xt for N epochs that is Eθt |yt,n [l(d1 , θt )] = +∞ (Eθt |yt,n [l(d2 , θt )] = +∞). Con-
concluded by C1 to be underwhelming. So, model training sequently, terminal decision d1 (d2 ) is never recommended.
can be early-stopped in epoch nt = n. More importantly, The above definitions are plugged into (2) to derive the mini-
both C1 and C2 are necessary for theoretically guaranteeing mum expected loss ρt,n (yt,n ) in epoch n = N0 + 1, . . . , N .
the no-regret performance of BO-BOS (Section 4). Our formulation of the BOS problem (2) for early stopping
3.2. BOS for Early Stopping in BO in BO can be solved using an adapted approximate back-
ward induction algorithm: To account for Assumption 1b,
Let the currently found maximum in iterations 1, . . . , t − 1
∗ a kernel with a prior bias towards exponentially decaying
be denoted as yt−1 , maxt0 ∈{1,...,t−1} yt0 . In the context
learning curves (Swersky et al., 2014) is used to fit a GP
of early stopping in BO, BOS has to decide in each epoch
model to the validation errors 1 − yt,N0 of the ML model
n of iteration t whether model training under xt for N
trained for N0 initial epochs. Samples are then drawn from
epochs will produce a validation accuracy not more than
the resulting GP posterior belief for forward simulation of
the currently found maximum (offset by a noise correction
∗ sample paths from epochs N0 + 1 to N , which are used
term ξt ), i.e., f ([xt , N ]) ≤ yt−1 − ξt where y0∗ , 0 and
to estimate the Pr(θt |yt,n ) and Pr(yt,n+1 |yt,n ) terms neces-
ξt is defined later in Theorem 1. To achieve this, the ter-
sary for approximate backward induction. Following some
minal decisions d1 and d2 and the continuation decision
applications of BOS (Jiang et al., 2013; Müller et al., 2007),
d0 in BOS are defined as follows: d1 stops and concludes
∗ the average validation error is used as the summary statistic.
that f ([xt , N ]) ≤ yt−1 − ξt , d2 stops and concludes that
∗ Our adapted approximate backward induction algorithm is
f ([xt , N ]) > yt−1 −ξt , and d0 continues model training for
explained in detail in Appendix B. Note that the use of the
one more epoch. Then, the event θ (Section 2.2) becomes
( kernel favoring exponentially decaying learning curves in
∗
θt,1 if f ([xt , N ]) ≤ yt−1 − ξt , generating the forward simulation samples is critical for
θt =
θt,2 otherwise . incorporating our prior knowledge about the behavior of
learning curves, which gives BO-BOS an advantage over
multi-fidelity BO algorithms which do not exploit this prior
Algorithm 1 BO-BOS knowledge, thus contributing to the favorable performance
1: for t = 1, 2, . . . , T do √ of BO-BOS.
2: xt ← arg maxx µt−1 ([x, N ]) + βt σt−1 ([x, N ])
3: Train model using xt for N0 epochs to yield yt,N0 After our BOS problem is solved, Bayes-optimal decision
4: n ← N0 rules are obtained and used by C1 in BO-BOS (Algorithm 1):
5: Solve BOS problem (2) to obtain decision rules Specifically, after model training to yield validation accu-
6: repeat racyPynt,n (line 8), the summary statistic is first updated
7: n←n+1 to n0 =1 yt,n0 /n and the BOS decision rule in epoch n
8: Continue model training using xt to yield yt,n recommends a corresponding optimal decision. If the rec-
9: until (n = N ) ∨ (C1 ∧ C2) ommended decision is d1 (i.e., stopping and concluding
∗
10: nt = n f ([xt , N ]) ≤ yt−1 − ξt ), then model training under xt
11: Update GP posterior belief (1) with zt = [xt , nt ] and is early-stopped in epoch n (assuming that C2 is satis-
yt = yt,nt fied). Otherwise, model training continues under xt for
one more epoch and the above procedure is repeated in
epoch n + 1 until the last epoch n = N is reached. Note Theorem 2 below states how the BOS parameters should be
that terminal decision d2 (i.e., stopping and concluding chosen to make BO-BOS asymptotically no-regret:
∗
f ([xt , N ]) > yt−1 − ξt ) does not align with the BO ob-
Theorem 2. In Theorem 1, if K1,t is an increasing sequence
jective of sequentially maximizing f . So, when the recom-
such that K1,1 ≥ K2 + cd0 and that K1,t goes to +∞ in
mended decision is d2 , there is no early stopping and model
finite number of BO iterations, then, with probability of at
training continues under xt for one more epoch.
least 1 − δ − δ 0 − δ 00 , E[ST ] goes to 0 asymptotically.
4. Theoretical Analysis The proof of both theorems is presented in Appendix C.

The first term in the upper bound of E[ST ] in Theorem 1
The goal of the theoretical analysis is to characterize the matches that of the simple regret of the GP-UCB algorithm
growth of the simple regret ST of the proposed BO-BOS (up to the constant κ). Note that the theoretical results rely
algorithm and thus show how the algorithm should be de- on the exact solution of the BOS problems; however, in
signed in order for ST to asymptotically go to 0, i.e., for practice, a trade-off exists between the quality of the approx-
the algorithm to be no regret. To account for the addi- imate backward induction and the computational efficiency.
tional uncertainty introduced by BOS, we analyze the ex- In particular, increasing the number of forward simulation
pected regret, in which the expectation is taken with re- samples and making the grid of summary statistics more
spect to the posterior probabilities from the BOS algorithms: fine-grained both lead to better approximation quality, while
∗
Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt ). increasing the computational cost. Recommended approxi-
We make the following assumption on the smoothness of mation parameters that work well in all our experiments and
the objective function f : thus strike a reasonable balance between these two aspects
are given in Section 5.
Assumption 2. Assume for the kernel k, for some a and b,
Interestingly, the choice of an increasing K1,t sequence
P r(supz∈D |∂f /∂zj | > L) ≤ a exp(−(L/b)2 ) as required by Theorem 2 is well justified in terms of the
exploration-exploitation trade-off. As introduced in section
for j = 1, . . . , d where zj is the j-th component of input z. 3.2, K1 represents how much we would like to penalize
the BOS algorithm for falsely early-stopping (taking deci-
Assumption 2 is satisfied by some commonly used kernels sion d1 ). Therefore, increasing values of K1 implies that,
such as the Square Exponential kernel and Matérn kernel as the BO-BOS algorithm progresses, we become more
with ν > 2 (Srinivas et al., 2010). For simplicity, we as- and more cautious at early-stopping. In other words, the
sume that the underlying domain D is discrete, i.e. |D| < ∞. preference of BOS for early stopping diminishes over BO
However, it is straightforward to extend the analysis to gen- iterations. Interestingly, this corresponds to sequentially
eral compact domain by following similar analysis strategies shifting our preference from exploration (using small num-
as those in Appendix A.2. of (Srinivas et al., 2010). Theo- ber of epochs) to exploitation (using large number of epochs)
rem 1 below shows an upper bound on the expected simple throughout all runs of the BOS algorithms, which is an im-
regret of the BO-BOS algorithm. portant decision-making principle followed by many sequen-
Theorem 1. Suppose that Assumptions 1 and 2 hold. Let tial decision-making algorithms such as BO, multi-armed
δ, δ 0 , δ 00 ∈ (0, 1), βt , 2 log(|D|t2 π 2 /(6δ)), and τT , bandit, reinforcement learning, among others.
t=1 1nt <N be the number of BO iterations in which early
PT
Another intriguing interpretation of the theoretical results
stopping happens from iterations 1 to T . Let κ be the pa- is that the growth rate of the K1,t sequence implicitly de-
rameter used in C2. At iteration t, the BOS algorithm is run termines the trade-off between faster convergence of the
with the corresponding fixed cost parameters K2 and cd0 , BO algorithm (smaller number of BO iterations) and more
as well as iteration-dependent
p cost parameter K1,t , ξ1 , 0, computational saving in each BO iteration (smaller number
and ξt , 2σ 2 log(π 2 t2 (t − 1)/(6δ 00 )) for t > 1. Then of training epochs on average). In particular, if K1,t grows
∀T ≥ 1, with probability of at least 1 − δ − δ 0 − δ 00 , quickly, the second and third terms in the upper bound in
√ PT r Theorem 1 both decay fast, since ηt is inversely related to
κ T C1 βT γT ηt 1 da K1,t and large penalty for early stopping results in small τT ;
E[ST ] ≤ + t=1 + N b log 0 τT
T T T δ as a result, a large number of hyperparameters are run with
K +c
N epochs and the resulted BO-BOS algorithm behaves simi-
in which ηt , K 2
1,t
d0
, C1 = 8/log(1 + σ −2 ), γT is the larly to GP-UCB, which is corroborated by the upper bound
maximum information gain about the function f from any on E[ST ] in Theorem 1 since the first term dominates. On
Q of observations of size T , and the expectation
set
∗
is w.r.t. the other hand, if the K1,t sequence grows slowly, then the
t∈{t0 |t0 =1,...,T,nt0 <N } Pr(f ([x t , N ]) > yt−1 − ξ t |yt,nt ) second and third terms in Theorem 1 decay slowly; conse-
used in the BOS algorithm. quently, these two terms dominate the regret and the resulted
algorithm early-stops very often, thus leading to smaller

number of epochs on average, at the potential expense of
requiring more BO iterations. Furthermore, the constant κ
used in C2 also implicitly encodes our relative preference
for early-stopping. Specifically, large values of κ favor early
stopping by relaxing C2: σt−1 ([xt , n]) ≥ σt−1 ([xt , N ])/κ;
however, more early-stopped function evaluations might
incur larger number of required BO iterations as can be ver-
ified by the fact that larger κ increases the first regret term Figure 1. Best-found validation error of logistic regression v.s. the
total number of epochs (averaged over 10 random initializations).
in Theorem 1 (which matches the regret of GP-UCB). In
K2 = 99 and cd0 = 1 are fixed; K1,1 = 100 for all four
practice, as a result of the above-mentioned trade-offs, the BO-BOS algorithms; for t > 1, the different K1,t sequences
best choices of the BOS parameters and κ are application- K1,t−1
are: (K1,t )a = 0.89
K1,t−1
; (K1,t )b = 0.95
K1,t−1
; (K1,t )c = 0.99 ;
dependent. In addition, to ensure desirable behaviors of K1,t−1
(K1,t )d = 1.0 = K1,1 .
BOS, the BOS parameters should be chosen with additional
care. In particular, cd0 should be small, whereas K1,t should
the BOS algorithms, and as the incumbent validation error
be of similar order with K2 initially. In Section 5, we recom-
decreases, the preference of the algorithm for early stopping
mend some parameters that work well in all our experiments
will increase, thus preventing the resulting algorithm from
and thus are believed to perform robustly in practice.
beginning to exploit (running promising hyperparameter
settings with N epochs). This observation demonstrates the
5. Experiments and Discussion necessity of having an increasing sequence of K1,t values,
thus substantiating the practical relevance of our theoretical
The performance of BO-BOS is empirically compared with
analysis (Theorem 2). The sequence (K1,t )b , as well as
four other hyperparameter optimization algorithms: GP-
the values of K2 = 99 and cd0 = 1, will be used in the
UCB (Srinivas et al., 2010), Hyperband (Li et al., 2017),
following experiments if not further specified.
multi-fidelity BO algorithm called BO with continuous ap-
proximations (BOCA) (Kandasamy et al., 2017), and GP- 5.2. Hyperparameter Optimization of CNN
UCB with learning curve prediction using an ensemble
of Bayesian parametric regression models (LC Prediction) In this section, we tune six hyperparameters of CNN using
(Domhan et al., 2015). Freeze-thaw BO is not included in two image datasets: CIFAR-10 (Krizhevsky, 2009) and
the comparison since its implementation details are compli- Street View House Numbers (SVHN) (Netzer et al., 2011).
cated and not fully available. We empirically evaluate the Note that the goal of the experiments is not to compete with
performance of BO-BOS in hyperparameter optimization of the state-of-the-art models, but to compare the efficiency
logistic regression (LR) and convolutional neural networks of different hyperparameter tuning algorithms, so no data
(CNN), respectively, in Sections 5.1 and 5.2, and demon- augmentation or advanced network architectures are used.
strate its generality in two other interesting applications in As shown in Figures 2a and 2b, BO-BOS outperforms all
Section 5.3. Due to lack of space, additional experimental other algorithms under comparison in terms of the run-time
details are deferred to Appendix D. efficiency. Note that the horizontal axis, which represents
the wall-clock time, includes all the time spent during the
5.1. Hyperparameter Optimization of LR algorithms, including the running time of machine learning
We first tune three hyperparameters of LR trained on the models, the approximate backward induction used to solve
MNIST image dataset. Although both the K1,t sequence BOS, etc. Although Hyperband is able to quickly reduce the
and κ determine the trade-off between the number of BO validation error in the initial stage, it eventually converges
iterations and the number of epochs on average as mentioned to sub-optimal hyperparameter settings compared with both
in section 4, for simplicity, we fix κ = 2 and investigate the GP-UCB and BO-BOS. Similar findings have been reported
impact of different sequences of K1,t values. in previous works (Klein et al., 2017; Falkner et al., 2018)
and they might be attributed to the pure-exploration nature
As shown in Fig. 1, the sequences (K1,t )a , (K1,t )b and of the algorithm. Our BO-BOS algorithm, which trades
(K1,t )c lead to similar performances, all of which outper- off exploration and exploitation, is able to quickly surpass
form GP-UCB. On the other hand, the algorithm with fixed the performance of Hyperband and eventually converges to
K1 values ((K1,t )d ), despite having fast performance im- significantly better hyperparameter settings. In addition, we
provement initially, eventually finds a worse hyperparameter also ran the BOHB algorithm (Falkner et al., 2018) which
setting than all other algorithms. This undesirable perfor- combines Hyperband and BO; however, BOHB did not
mance results from the fact that fixed K1 values give con- manage to reach comparable performance with the other
stant penalty to falsely early-stopping throughout all runs of algorithms in these two tasks. We observed that the sub-
is usually estimated by running the agent in the environment

sequentially for a fixed number of steps and calculating the
cumulative rewards (Wilson et al., 2014). Thus, the sequen-
tial nature of policy evaluation makes BO-BOS an excellent
fit to improve the efficiency of BO for policy search in RL.
We apply our algorithm to the Swimmer-v2 task from Ope-
nAI Gym, MuJoCo (Brockman et al., 2016; Todorov et al.,
2012), and use a linear policy consisting of 16 parame-
(a) CIFAR-10.
ters. Each episode consists of 1000 steps, and we treat
every m consecutive steps as one single epoch such that
N = 1000/m. Direct application of BO-BOS in this task
is inappropriate since the growth pattern of cumulative re-
wards differs significantly from the evolution of the learning
curves of ML models (Appendix D.3). Therefore, the re-
wards are discounted (by γ) when calculating the objective
function, because the pattern of discounted return (cumu-
lative rewards) bears close resemblance to that of learning
(b) SVHN. curves. Note that although the value of the objective func-
Figure 2. Best-found validation error of CNN v.s. run-time (aver- tion is the discounted return, we also record and report the
aged over 30 random initializations). corresponding un-discounted return, which is the ultimate
objective to be maximized. As a result, N and γ should
optimal behavior of Hyperband and BOHB observed here be chosen such that the value of discounted return faith-
can be alleviated if the search space of hyperparameters is fully aligns with its un-discounted counterpart. Fig. 3 plots
chosen to be smaller, in which case the advantage of pure the best (un-discounted) return in an episode against the
exploration can be better manifested. The unsatisfactory total number of steps, in which BO-BOS (with N = 50
performance of BOCA might be explained by the fact that and γ = 0.9) outperforms GP-UCB (for both γ = 0.9 and
it does not make use of the intermediate validation errors γ = 1.0). The observation that the solid red line shows
when selecting the number of epochs. In contrast, BO-BOS better returns than the two dotted red lines might be because
takes into account the observations after each training epoch overly small γ (0.75) and overly large N (100) both enlarge
when choosing the optimal stopping time, and thus is able the disparity between discounted and un-discounted returns
to make better-informed decisions. Moreover, since BOCA since they both downplay the importance of long-term dis-
is designed for general scenarios, the useful assumption counted rewards. Moreover, not discounting the rewards
of monotonic learning curve utilized by BO-BOS is not (γ = 1) leads to poor performance, corroborating the ear-
exploited; therefore, BOCA is expected to perform better lier analysis motivating the use of discounted rewards. The
if the fidelity levels result from data sub-sampling. LC results demonstrate that even though BO-BOS is not imme-
Prediction performs similarly to GP-UCB, which might be diately applicable in the original problem, it can still work
because the predicted final values of the learning curves effectively if the problem is properly transformed, which
are used as real observations in the GP surrogate function substantiates the general applicability of BO-BOS. More-
(Domhan et al., 2015), thus invalidating the theoretical guar- over, we also applied Hyperband in this task, but it failed
antee and deteriorating the convergence of GP-UCB, which to converge to a comparable policy to the other methods
offsets the computational saving provided by early stopping. (achieving an average return of around 2.5), which further
BO-BOS, on the other hand, offers theoretically guaranteed supports our earlier claim stating that Hyperband under-
convergence, thus allowing explicit control over the trade- performs when the search space is large because there are
off between the speed of convergence and the reduction in significantly more parameters (16) than the previous tasks.
the average number of epochs, as discussed in Section 4.
Interestingly, both BO-BOS and GP-UCB use significantly
5.3. Novel Applications of the BO-BOS Algorithm less steps to substantially outperform the benchmarks2
achieved by some recently developed deep RL algorithms,
5.3.1. P OLICY S EARCH FOR RL in which the best-performing algorithm (Deep Deterministic
Thanks to its superb sample efficiency, BO has been found Policy Gradients) achieves an average return of around 120.
effective for policy search in RL (Wilson et al., 2014), espe- Although the simplicity of the task might have contributed
cially when policy evaluation is costly such as gait optimiza- to their overwhelming performances, the results highlight
tion for robots (Lizotte et al., 2007) and vehicle navigation 2
https://spinningup.openai.com/en/latest/spinningup/bench.html
(Brochu et al., 2010). In policy search, the return of a policy
can be naturally solved by BO-BOS by simply replacing the

sequential training of ML models with forward selection.
In this task, we tune four hyperparameters of the gradi-
ent boosting model (XGBoost (Chen & Guestrin, 2016))
trained on an email spam dataset. We additionally compare
with Hyperband since it was previously applied to random
feature approximation in kernel methods (Li et al., 2017).
As shown in Fig. 4, BO-BOS again delivers the best per-
Figure 3. Best-found return (averaged over 5 episodes) v.s. the total formance in this application. Consistent with Fig. 2, the
number of steps of the robot in the environment (averaged over 30 wall-clock time includes all the time incurred during the
random initializations) using the Swimmer-v2 task. The BO-BOS algorithms. The K1,t sequence is made smaller than be-
algorithm is run with different values of N (the maximum number fore: K1,t = K1,t−1 /0.99, because in this setting, more
of epochs) and γ (the discount factor). aggressive early stopping is needed for BO-BOS to show
its advantage. Although Hyperband works well for random
feature approximation (Li et al., 2017), it does not perform
favourably when applied to more structured feature selection
techniques. Both hyperparameter optimization and feature
selection have been shown to be effective for enhancing the
performance of ML models. However, performing either
of them in isolation may lead to sub-optimal performance
since their interaction is un-exploited. Our results suggest
that BO-BOS can effectively improve the efficiency of joint
hyperparameter tuning and feature selection, making the
Figure 4. Best-found validation error of XGBoost v.s. run-time
combined usage of these two pre-processing techniques a
(averaged over 30 random initializations), obtained using joint
hyperparameter tuning and feature selection. more practical choice.
the considerable potential of BO-based policy search algo- 6. Conclusion

rithms for RL. Note that following the common practice
in RL, Fig. 3 is presented in terms of the total number of This paper describes a unifying BO-BOS framework that in-
steps instead of run-time; we believe the additional com- tegrates BOS into BO in a natural way to derive a principled
putation required by BOS can be easily overshadowed in mechanism for optimally stopping hyperparameter evalua-
large-scale RL tasks, as demonstrated in the previous ex- tions during BO. We analyze the regret of the algorithm, and
periments. Most modern RL algorithms rely on enormous derive the BOS parameters that make the resulting BO-BOS
number of samples, making their applications problematic algorithm no-regret. Applications of BO-BOS to hyperpa-
when sample efficiency is of crucial importance (Arulku- rameter tuning of ML models, as well as two other novel
maran et al., 2017). As shown above, BO-BOS achieves applications, demonstrate the practical effectiveness of the
significantly better results than popular deep RL algorithms algorithm. For future work, we plan to use BO-BOS in
while at the same time being far more sample-efficient, thus other applications with iterative function evaluations and
potentially offering practitioners a practically feasible op- generalize BO-BOS to the batch3 (Daxberger & Low, 2017)
tion in solving large-scale RL problems. and high-dimensional (Hoang et al., 2018) BO settings, as
well as to the non-myopic context by appealing to existing
5.3.2. J OINT H YPERPARAMETER T UNING AND F EATURE literature on nonmyopic BO (Ling et al., 2016) and active
S ELECTION learning (Cao et al., 2013; Hoang et al., 2014a;b; Low et al.,
2008; 2009; 2011; 2014a). For applications with a huge
Beside hyperparameter tuning, feature selection is another budget of function evaluations, we like to couple BO-BOS
important pre-processing step in ML, to lower computa- with the use of distributed/decentralized (Chen et al., 2012;
tional cost and boost model performance (Hall, 2000). There 2013a;b; 2015; Hoang et al., 2016; 2019; Low et al., 2015;
are two types of feature selection techniques: filter and Ouyang & Low, 2018) or online/stochastic (Hoang et al.,
wrapper methods; the wrapper methods, such as forward 2015; 2017; Low et al., 2014b; Xu et al., 2014; Yu et al.,
selection (which starts with an empty feature set and in each 2019) sparse GP models to represent the belief of the un-
iteration greedily adds the feature with largest performance known objective function efficiently.
improvement), have been shown to perform well, although 3
computationally costly (Hall & Smith, 1999). Interestingly, A closely related counterpart is batch active learning (Low
et al., 2012; Ouyang et al., 2014).
as a result of the sequential nature of wrapper methods, they
Acknowledgements Chen, T. and Guestrin, C. XGBoost: A scalable tree boost-

ing system. In Proc. KDD, pp. 785–794, 2016.
This research is supported by the National Research Foun-
dation, Prime Minister’s Office, Singapore under its Cam- Davis, G. A. and Cairns, R. D. Good timing: The economics
pus for Research Excellence and Technological Enterprise of optimal stopping. Journal of Economic Dynamics and
(CREATE) programme and Singapore Ministry of Educa- Control, 36(2):255–265, 2012.
tion Academic Research Fund Tier 2, MOE2016-T2-2-156.
Daxberger, E. and Low, K. H. Distributed batch Gaussian
process optimization. In Proc. ICML, pp. 951–960, 2017.
References
Arulkumaran, K., Deisenroth, M. P., Brundage, M., and Dheeru, D. and Karra Taniskidou, E. UCI machine learning
Bharath, A. A. Deep reinforcement learning: A brief repository, 2017. URL http://archive.ics.uci.
survey. IEEE Signal Processing Magazine, 17(6):26–38, edu/ml.
2017. Domhan, T., Springenberg, J. T., and Hutter, F. Speeding
Baker, B., Gupta, O., Raskar, R., and Naik, N. Accelerating up automatic hyperparameter optimization of deep neural
neural architecture search using performance prediction. networks by extrapolation of learning curves. In Proc.
In Proc. NIPS Workshop on Meta-Learning, 2017. IJCAI, pp. 3460–3468, 2015.
Falkner, S., Klein, A., and Hutter, F. BOHB: Robust and

Brochu, E., Cora, V. M., and de Freitas, N. A tutorial
efficient hyperparameter optimization at scale. In Proc.
on Bayesian optimization of expensive cost functions,
ICML, pp. 1436–1445, 2018.
with application to active user modeling and hierarchical
reinforcement learning. arXiv:1012.2599, 2010. Ferguson, T. S. Optimal stopping and applications.
2006. URL https://www.math.ucla.edu/
Brockman, G., Cheung, V., Pettersson, L., Schneider, J.,
˜ tom/Stopping/Contents.html.
Schulman, J., Tang, J., and Zaremba, W. OpenAI Gym.
arXiv:1606.01540, 2016. Friedman, J. H. Greedy function approximation: A gradi-
ent boosting machine. Ann. Statistics, 29(5):1189–1232,
Brockwell, A. E. and Kadane, J. B. A gridding method for
2001.
Bayesian sequential decision problems. J. Computational
and Graphical Statistics, 12(3):566–584, 2003. Hall, M. A. Correlation-based feature selection of discrete
and numeric class machine learning. In Proc. ICML, pp.
Cao, N., Low, K. H., and Dolan, J. M. Multi-robot infor-
359–366, 2000.
mative path planning for active sensing of environmental
phenomena: A tale of two algorithms. In Proc. AAMAS, Hall, M. A. and Smith, L. A. Feature selection for machine
pp. 7–14, 2013. learning: Comparing a correlation-based filter approach
to the wrapper. In Proc. FLAIRS conference, pp. 235–239,
Chen, J., Low, K. H., Tan, C. K.-Y., Oran, A., Jaillet, P., 1999.
Dolan, J. M., and Sukhatme, G. S. Decentralized data
fusion and active sensing with mobile sensors for model- Hoang, Q. M., Hoang, T. N., and Low, K. H. A generalized
ing and predicting spatiotemporal traffic phenomena. In stochastic variational Bayesian hyperparameter learning
Proc. UAI, pp. 163–173, 2012. framework for sparse spectrum Gaussian process regres-
sion. In Proc. AAAI, pp. 2007–2014, 2017.
Chen, J., Cao, N., Low, K. H., Ouyang, R., Tan, C. K.-
Y., and Jaillet, P. Parallel Gaussian process regression Hoang, T. N., Low, K. H., Jaillet, P., and Kankanhalli,
with low-rank covariance matrix approximations. In Proc. M. Active learning is planning: Nonmyopic -Bayes-
UAI, pp. 152–161, 2013a. optimal active learning of Gaussian processes. In Proc.
ECML/PKDD Nectar Track, pp. 494–498, 2014a.
Chen, J., Low, K. H., and Tan, C. K.-Y. Gaussian process-
based decentralized data fusion and active sensing for Hoang, T. N., Low, K. H., Jaillet, P., and Kankanhalli, M.
mobility-on-demand system. In Proc. RSS, 2013b. Nonmyopic -Bayes-optimal active learning of Gaussian
processes. In Proc. ICML, pp. 739–747, 2014b.
Chen, J., Low, K. H., Jaillet, P., and Yao, Y. Gaussian pro-
cess decentralized data fusion and active sensing for spa- Hoang, T. N., Hoang, Q. M., and Low, K. H. A unifying
tiotemporal traffic modeling and prediction in mobility- framework of anytime sparse Gaussian process regression
on-demand systems. IEEE Trans. Autom. Sci. Eng., 12: models with stochastic variational inference for big data.
901–921, 2015. In Proc. ICML, pp. 569–578, 2015.
Hoang, T. N., Hoang, Q. M., and Low, K. H. A dis- Low, K. H., Dolan, J. M., and Khosla, P. Adaptive multi-
tributed variational inference framework for unifying par- robot wide-area exploration and mapping. In Proc. AA-
allel sparse Gaussian process regression models. In Proc. MAS, pp. 23–30, 2008.
ICML, pp. 382–391, 2016.
Low, K. H., Dolan, J. M., and Khosla, P. Information-
Hoang, T. N., Hoang, Q. M., and Low, K. H. Decentral- theoretic approach to efficient adaptive path planning for
ized high-dimensional Bayesian optimization with factor mobile robotic environmental sensing. In Proc. ICAPS,
graphs. In Proc. AAAI, pp. 3231–3238, 2018. pp. 233–240, 2009.
Hoang, T. N., Hoang, Q. M., Low, K. H., and How, J. P. Col- Low, K. H., Dolan, J. M., and Khosla, P. Active Markov
lective online learning of Gaussian processes in massive information-theoretic path planning for robotic environ-
multi-agent systems. In Proc. AAAI, 2019. mental sensing. In Proc. AAMAS, pp. 753–760, 2011.
Jiang, F., Jack Lee, J., and Müller, P. A Bayesian decision-
Low, K. H., Chen, J., Dolan, J. M., Chien, S., and Thomp-
theoretic sequential response-adaptive randomization de-
son, D. R. Decentralized active robotic exploration and
sign. Statistics in Medicine, 32(12):1975–1994, 2013.
mapping for probabilistic field classification in environ-
Kandasamy, K., Dasarathy, G., Oliva, J. B., Schneider, J., mental sensing. In Proc. AAMAS, pp. 105–112, 2012.
and Póczos, B. Gaussian process bandit optimisation with
multi-fidelity evaluations. In Proc. NIPS, pp. 992–1000, Low, K. H., Chen, J., Hoang, T. N., Xu, N., and Jaillet,
2016. P. Recent advances in scaling up Gaussian process pre-
dictive models for large spatiotemporal data. In Ravela,
Kandasamy, K., Dasarathy, G., Schneider, J., and Póczos, S. and Sandu, A. (eds.), Dynamic Data-Driven Environ-
B. Multi-fidelity Bayesian optimisation with continuous mental Systems Science: First International Conference,
approximations. In Proc. ICML, pp. 1799–1808, 2017. DyDESS 2014, pp. 167–181. LNCS 8964, Springer Inter-
national Publishing, 2014a.
Klein, A., Falkner, S., Springenberg, J. T., and Hutter, F.
Learning curve prediction with Bayesian neural networks. Low, K. H., Xu, N., Chen, J., Lim, K. K., and Özgül, E. B.
In Proc. ICLR, 2017. Generalized online sparse Gaussian processes with ap-
Krizhevsky, A. Learning multiple layers of features from plication to persistent mobile robot localization. In Proc.
tiny images. Master’s thesis, Department of Computer ECML/PKDD Nectar Track, pp. 499–503, 2014b.
Science, University of Toronto, 2009.
Low, K. H., Yu, J., Chen, J., and Jaillet, P. Parallel Gaussian
LeCun, Y., Bengio, Y., and Hinton, G. Deep learning. Na- process regression for big data: Low-rank representation
ture, 521(7553):436–444, 2015. meets Markov approximation. In Proc. AAAI, pp. 2821–
2827, 2015.
Lewis, R. J. and Berry, D. A. Group sequential clinical tri-
als: A classical evaluation of Bayesian decision-theoretic Müller, P., Berry, D. A., Grieve, A. P., Smith, M., and
designs. J. American Statistical Association, 89(428): Krams, M. Simulation-based sequential Bayesian design.
1528–1534, 1994. J. Statistical Planning and Inference, 137(10):3140–3150,
2007.
Li, L., Jamieson, K., DeSalvo, G., Rostamizadeh, A., and
Talwalkar, A. Hyperband: A novel bandit-based approach Netzer, Y., Wang, T., Coates, A., Bissacco, A., Wu, B.,
to hyperparameter optimization. JMLR, 18(1):6765–6816, and Ng, A. Y. Reading digits in natural images with
2017. unsupervised feature learning. In Proc. NIPS Workshop
Ling, C. K., Low, K. H., and Jaillet, P. Gaussian process on Deep Learning and Unsupervised Feature Learning,
planning with Lipschitz continuous reward functions: To- 2011.
wards unifying Bayesian optimization, active learning,
Ouyang, R. and Low, K. H. Gaussian process decentral-
and beyond. In Proc. AAAI, pp. 1860–1866, 2016.
ized data fusion meets transfer learning in large-scale
Lizotte, D. J., Wang, T., Bowling, M. H., and Schuurmans, distributed cooperative perception. In Proc. AAAI, pp.
D. Automatic gait optimization with Gaussian process 3876–3883, 2018.
regression. In IJCAI, pp. 944–949, 2007.
Ouyang, R., Low, K. H., Chen, J., and Jaillet, P. Multi-robot
Longstaff, F. A. and Schwartz, E. S. Valuing American active sensing of non-stationary Gaussian process-based
options by simulation: A simple least-squares approach. environmental phenomena. In Proc. AAMAS, pp. 573–
The Review of Financial Studies, 14(1):113–147, 2001. 580, 2014.
Powell, W. B. and Ryzhov, I. O. Optimal learning. John

Wiley & Sons, 2012.
Rasmussen, C. E. and Williams, C. K. I. Gaussian processes
for machine learning. MIT Press, 2006.
Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and
De Freitas, N. Taking the human out of the loop: A
review of Bayesian optimization. Proceedings of the
IEEE, 104(1):148–175, 2016.
Silver, D., Huang, A., Maddison, C. J., Guez, A., Sifre, L.,
van den Driessche, G., Schrittwieser, J., Antonoglou, I.,
Panneershelvam, V., Lanctot, M., Dieleman, S., Grewe,
D., Nham, J., Kalchbrenner, N., Sutskever, I., Lillicrap, T.,
Leach, M., Kavukcuoglu, K., Graepel, T., and Hassabis,
D. Mastering the game of Go with deep neural networks
and tree search. Nature, 529(7587):484–489, 2016.
Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M.
Gaussian process optimization in the bandit setting: No
regret and experimental design. In Proc. ICML, pp. 1015–
1022, 2010.
Swersky, K., Snoek, J., and Adams, R. P. Freeze-thaw
Bayesian optimization. arXiv:1406.3896, 2014.
Todorov, E., Erez, T., and Tassa, Y. MuJoCo: A physics
engine for model-based control. In Proc. IEEE/RSJ IROS,
pp. 5026–5033, 2012.
Wathen, J. K. and Thall, P. F. Bayesian adaptive model
selection for optimizing group sequential clinical trials.
Statistics in Medicine, 27(27):5586–5604, 2008.
Wilson, A., Fern, A., and Tadepalli, P. Using trajectory
data to improve Bayesian optimization for reinforcement
learning. Journal of Machine Learning Research, 15(1):
253–282, 2014.
Wu, J. and Frazier, P. I. Continuous-fidelity Bayesian opti-
mization with knowledge gradient. In Proc. NIPS Work-
shop on Bayesian Optimization, 2017.
Xu, N., Low, K. H., Chen, J., Lim, K. K., and Özgül, E. B.
GP-Localize: Persistent mobile robot localization using
online sparse Gaussian process observation model. In
Proc. AAAI, pp. 2585–2592, 2014.
Yu, H., Hoang, T. N., Low, K. H., and Jaillet, P. Stochas-
tic variational inference for Bayesian sparse Gaussian
process regression. In Proc. IJCNN, 2019.
Zhang, Y., Hoang, T. N., Low, K. H., and Kankanhalli, M.
Near-optimal active learning of multi-output Gaussian
processes. In Proc. AAAI, pp. 2351–2357, 2016.
Zhang, Y., Hoang, T. N., Low, K. H., and Kankanhalli, M.
Information-based multi-fidelity Bayesian optimization.
In Proc. NIPS Workshop on Bayesian Optimization, 2017.
A. Approximate Backward Induction for Bayesian Optimal Stopping

In this section, we will present a commonly-used approximate backward induction algorithm for solving the BOS problem.
The algorithm uses summary statistics to compactly represent the posterior belief Pr(θt |yt,n ) which is computed from the
prior belief Pr(θt ) and the noisy outputs yt,n observed up till epoch n in iteration t.
In the approximate backward induction algorithm of Müller et al. (2007), the entire space of summary statistics is firstly
partitioned into a number of discrete intervals in each epoch, which results in a two-dimensional domain with one axis
being the number of epochs and the other axis representing the discretized intervals of the summary statistic (i.e., assuming
the summary statistic is one-dimensional). In the beginning, a number of sample paths are generated through forward
simulation: Firstly, a large number of samples are drawn from the prior belief Pr(θt ). Then, for each sample drawn from
Pr(θt ), an entire sample path is generated from epochs 1 to N through repeated sampling. In this manner, each sample path
leads to a curve in the 2-D domain and fully defines N posterior beliefs with one in each epoch. Starting from the last epoch
N , for each interval, the expected loss of a terminal decision d1 or d2 is evaluated for every sample path ending in this
interval (since each such sample path ends with a particular posterior belief in epoch N ), and their empirical average is
used to approximate the expected loss of the particular terminal decision for this interval. The minimum of the expected
losses among the two terminal decisions is the expected loss for this particular interval, which is equivalent to (2) except that
decision d0 is not available in the last epoch N .
Next, the algorithm proceeds backwards from epoch n = N − 1 all the way to epoch n = 1. In each epoch n, the expected
loss of each terminal decision is evaluated in the same way as that in the last epoch N , as described above. To evaluate the
expected loss of the continuation decision for an interval, for each sample path passing through this interval, the expected loss
for the interval that it passes through in the next epoch n + 1 is recorded and an average of all the recorded expected losses
in the next epoch n + 1 is summed with the cost cd0 of observing the noisy output yt,n+1 to yield the expected loss of the
continuation decision d0 for this particular interval; this is equivalent to approximating the Eyt,n+1 |yt,n [ρt,n+1 (yt,n+1 )]+cd0
term in (2) via Monte Carlo sampling of the posterior belief Pr(yt,n+1 |yt,n ). Following (2), the minimum of the expected
losses among all terminal and continuation decisions is the expected loss for this particular interval and the corresponding
decision is recorded as the optimal decision when the summary statistic falls into this interval. Then, the algorithm continues
backwards until epoch n = 1 is reached. After the algorithm has finished running, the optimal decision computed in every
pair of epoch and interval will form the optimal decision rules which serve as the output of the approximate backward
induction algorithm.
B. Approximate Backward Induction Algorithm for Solving BOS Problem in BO-BOS

In this section, we will describe the approximate backward induction algorithm for solving the BOS problem (line 5) in each
iteration of BO-BOS (Algorithm 1), which is adapted from the algorithm introduced in Appendix A.
To account for Assumption 1b in the approximate backward induction algorithm, we adopt the kernel k introduced
in (Swersky et al., 2014) to incorporate the inductive bias that the learning curve (in the form of validation error) of the ML
model is approximately exponentially decreasing in the number of training epochs, which can be expressed as
Z ∞
0 βα
k(n, n ) , exp(−λn) exp(−λn0 ) φ(λ) dλ = α (5)
0 (n + n0 + β)
for all epochs n, n0 = 1, . . . , N where φ is a probability measure over λ that is chosen to be a Gamma prior with parameters
α and β. The above kernel (5) is used to fit a GP model to the validation errors 1 − yt,N0 of the ML model trained using xt
for a fixed number N0 of initial epochs (e.g., N0 = 8 in all our experiments when N = 50), specifically, by computing
the values of parameters α and β in (5) via Bayesian update (i.e., assuming that the validation errors follow the Gamma
conjugate prior with respect to an exponential likelihood). Samples are then drawn from the resulting GP posterior belief for
forward simulation of sample paths from epochs N0 + 1 to N , which are used to estimate the Pr(θt |yt,n ) and Pr(yt,n+1 |yt,n )
terms necessary for approximate backward induction. Fig. 5 plots some of such sample paths and demonstrates that the GP
kernel in (5) can characterize a monotonic learning curve (Assumption 1b) well.
Following the practices in related applications of BOS (Brockwell & Kadane, 2003; Jiang et al., 2013; Müller et al., 2007),
the average validation error (or, equivalently, average validation accuracy) over epochs 1 to n is used as the summary
statistics. Firstly, the entire space of summary statistics is partitioned into a number of discrete intervals in each epoch,
which results in a two-dimensional domain with one axis being the number of epochs and the other axis representing the
0.3
Validation Error
0.2
0.1
0.0
5 10 15 20 25 30 35 40 45
Epochs
Figure 5. Forward simulation of some sample paths drawn from a GP posterior belief based on the kernel in (5).
discretized intervals of the summary statistic (i.e., average validation error). Next, a forward simulation of a large number
(i.e., 100, 000 in all our experiments) of sample paths is performed using the GP kernel in (5), as described above. Each
sample path corresponds to a curve in the 2-D domain. Starting from the last epoch N , for each interval, we consider all
sample paths ending in this interval and use the proportion of such sample paths with a validation accuracy (from model
training for N epochs) larger than the currently found maximum (offset by a noise correction term) to estimate the posterior
∗
probability Pr(θt = θt,2 |yt,N ) = Pr(f ([xt , N ]) > yt−1 − ξt |yt,N ), which is in turn used to evaluate the expected losses of
4
the terminal decisions d1 and d2 for this interval. The minimum of the expected losses among the two terminal decisions is
the expected loss for this particular interval.
Next, the algorithm proceeds backwards from epoch n = N − 1 all the way to epoch n = N0 + 1. In each epoch n, the
expected loss of each terminal decision is evaluated in the same way as that in the last epoch N , as described above. The
expected loss of the continuation decision d0 is evaluated in the same way as that in Appendix A: For each sample path
passing through an interval in epoch n, the expected loss for the interval that it passes through in the next epoch n + 1
is recorded and an average of all the recorded expected losses in the next epoch n + 1 is summed with the cost cd0 of
observing the validation accuracy yt,n+1 to yield the expected loss of the continuation decision d0 for this particular interval.
Note that this step is equivalent to approximating the Eyt,n+1 |yt,n [ρt,n+1 (yt,n )] term in (2) via Monte Carlo sampling of
the posterior belief Pr(yt,n+1 |yt,n ). Following (2), the minimum of expected losses among all terminal and continuation
decisions is the expected loss for this particular interval and the corresponding decision is recorded as the optimal decision
to be recommended when the summary statistic falls into this particular interval. Then, the algorithm continues backwards
until epoch n = N0 + 1 is reached. We present in Algorithm 2 the pseudocode for the above-mentioned approximate
backward induction algorithm for ease of understanding.
After solving our BOS problem for early stopping in BO using the approximate backward induction algorithm described
above, Bayes-optimal decision rules are obtained in every pair of epoch and interval. Fig. 6 shows an example of optimal
decision rules obtained from solving an instance of our BOS problem where the white, yellow, and red regions correspond
to recommending optimal continuation decision d0 and terminal decisions d1 and d2 , respectively. In particular, after model
training under xt to yield the validation error 1 − yt,n in epoch n, the summary statistic is updated to the average validation
error over epochs 1 to n. The updated summary statistic falls into an interval with a corresponding optimal decision to be
recommended. For example, Fig. 6 shows that if the summary statistic falls into the yellow region in any epoch n, then the
optimal terminal decision d1 is recommended to early-stop model training under xt (assuming that C2 is satisfied). If the
summary statistic falls into any other region, then model training continues under xt for one more epoch and the above
procedure is repeated in epoch n + 1 until the last epoch n = N is reached. This procedure, together with C2, constitutes
lines 6 to 9 in Algorithm 1.
4
In contrast to the approximate backward induction algorithm of Müller et al. (2007) (Appendix A), we employ a computationally
cheaper way to approximate the expected losses of the terminal decisions for an interval.
Algorithm 2 Approximate Backward Induction Algorithm for Solving BOS Problem in BO-BOS
1: Partition the domain of summary statistics into M discrete intervals
2: Train the ML model using xt for N0 epochs
3: Generate a large number of forward simulation samples using kernel (5)
4: Let n = N
5: for m = 1, 2, . . . , M do
6: Find all sample paths ending in interval m at epoch n, denoted as S
∗
7: Estimate Pr(θt = θt,2 |yt,n ) = Pr(f ([xt , N ]) > yt−1 − ξt |yt,n ) by the proportion of S that end up (after N epochs)
∗
having larger validation accuracy than yt−1 − ξt
8: Calculate the expected losses of the terminal decisions d1 and d2 using (4)
9: Use the minimum of these two expected losses as the expected loss of epoch n and interval m
10: for n = N − 1, N − 2, . . . , N0 + 1 do
11: for m = 1, 2, . . . , M do
12: Find all sample paths passing through interval m at epoch n, denoted as S
∗
13: Estimate Pr(θt = θt,2 |yt,n ) = Pr(f ([xt , N ]) > yt−1 − ξt |yt,n ) by the proportion of S that end up (after N
∗
epochs) having larger validation accuracy than yt−1 − ξt
14: Calculate the expected losses of the terminal decisions d1 and d2 using (4)
15: ld0 ,n+1,m = 0
16: for each sample path s in S do
17: ld0 ,n+1,m = ld0 ,n+1,m + the expected loss of the interval reached by s at epoch n + 1
18: ld0 ,n+1,m = ld0 ,n+1,m /|S|
19: Calculate the expected loss of the continuation decision d0 as: Eyt,n+1 |yt,n [ρt,n+1 (yt,n )] + cd0 = ld0 ,n+1,m + cd0
20: Use the minimum expected losses among d1 , d2 and d0 as the expected loss of epoch n and interval m (follow-
ing (2)), and record the corresponding decision as the optimal decision
1.0
Summary Statistics
0.75
0.5
0.25
0.0
1 6 11 16 21 26 31 36
Epochs
Figure 6. An example of optimal decision rules obtained from solving an instance of our BOS problem: White, yellow, and red regions
correspond to recommending optimal continuation decision d0 and terminal decisions d1 and d2 , respectively. The sample paths cannot
reach the black regions due to the use of the GP kernel in (5) for characterizing a monotonic learning curve (Assumption 1b).
C. Proof of Theorems 1 and 2

In this section, we prove the theoretical results in this paper.
C.1. Regret Decomposition

In this work, it is natural and convenient to define the instantaneous regret at step t as rt = f (z∗ ) − ft∗ , in which z∗ is
the location of the global maximum: z∗ = argmaxz f (z), and ft∗ is the maximum observed function value from iterations
1 to t: ft∗ = maxt0 =1,...,t f (zt0 ). Subsequently, the cumulative regret and simple regret after T iterations are defined as
PT
RT = t=1 rt and ST = mint=1,...,T rt respectively. As a result, as long as we can show that RT grows sub-linearly in
T , then we can conclude that the average regret RTT asymptotically goes to 0; therefore, ST vanishes asymptotically since
it is upper-bounded by the average regret: ST ≤ RTT . In contrast to the more commonly used definition of instantaneous
regret: rt = f (z∗ ) − f (zt ), the slightly modified definition introduced here is justified in the sense that the induced
definition of simple regrets, which is the ultimate goal of the theoretical analysis, obtained in both cases are equivalent, i.e.,
mint=1,...,T f (z∗ ) − ft∗ = mint=1,...,T f (z∗ ) − f (zt ).
The instantaneous regret defined above can be further decomposed as
rt =f (z∗ ) − ft∗ = f (z∗ ) − 0 max f (zt0 )

t =1,...,t
(6)
=f (z∗ ) − max{ft−1
∗
, f (zt )}
Note that in our algorithm, the BO iterations can be divided into two types: 1) t+ such that nt+ = N : those iterations that
are not early-stopped; and 2) t− such that nt− < N : those that are early-stopped. For all t+ , it follows from Equation 6 that
rt = f (z∗ ) − max{ft−1∗
, f (zt )} ≤ f (z∗ ) − f (zt ) = f (z∗ ) − f ([xt , nt ]) = f (z∗ ) − f ([zt , N ]) , rt+ ; for all t− , from
Equation 6, we have that rt = f (z∗ ) − max{ft−1 ∗
P, f (zt )} ≤Pf (z∗ ) − ft−1
∗
, rt− . In the following, we will focus on the
analysis of the sum of all rt+ and all rt− : RT = t+ rt+ + t− rt− . As a result of the definition, RT0 is an upper bound
0
of RT , therefore, sub-linear growth of RT0 implies that RT also grows sub-linearly.
Next, note that for all t− such that nt− < N (when xt is early-stopped),
rt− = f (z∗ ) − ft−1

∗
= f (z∗ ) − f ([xt , N ]) + f ([xt , N ]) − ft−1
∗
(1) (7)
≤ f (z∗ ) − f ([xt , nt ]) + f ([xt , N ]) − ft−1
∗
in which (1) results from Assumption 1. As a result, RT0 can be re-written as

(1) X X
RT0 = [f (z∗ ) − f ([xt , N ])] + [f (z∗ ) − ft−1
∗
]
{t|nt =N } {t|nt <N }
X X X
∗
= [f (z ) − f ([xt , N ])] + [f (z∗ ) − f ([xt , N ])] + ∗
[f ([xt , N ]) − ft−1 ]
{t|nt =N } {t|nt <N } {t|nt <N }
(2) X X X
≤ [f (z∗ ) − f ([xt , N ])] + [f (z∗ ) − f ([xt , nt ])] + ∗
[f ([xt , N ]) − ft−1 ]
{t|nt =N } {t|nt <N } {t|nt <N }
(8)
T
(3) X X
= [f (z∗ ) − f ([xt , nt ])] + ∗
[f ([xt , N ]) − ft−1 ]
t=1 {t|nt <N }
T
X X
, rt,1 + rt,2
t=1 {t|nt <N }
, RT,1 + RT,2
in which (1) makes use of the definition of RT0 , (2) results from Equation 7 and (3) follows by combining the first two terms
on the previous line. The first term following (3) of Equation 8 is summed over all time steps, whereas the second term is
only summed over those time steps that are early-stopped (nt < N ). As mentioned earlier, in the sequel, we will attempt to
prove an upper bound on the expected value of RT0 ,
E[RT0 ] ≤ E[RT,1 ] + E[RT,2 ] (9)
in which the expectation is taken with respect to the posterior probabilities used in the BOS problems, corresponding to
∗
those iterations that are early-stopped: Πt∈{t0 |t0 =1,...,T,nt0 <N } Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt ). Note that the probability
distributions are independent across all the early-stopped iterations, therefore, for each early-stopped iteration t, the
∗
expectations of both rt,1 and rt,2 are only taken over the specific distribution: Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt ); whereas
for each not-early-stopped iteration t, E[rt,1 ] = rt,1 (whereas rt,2 is absent). In the next two sections, we will prove upper
bounds on E[RT,1 ] and E[RT,2 ] respectively.
C.2. Upper Bound on E[RT,1 ]

In this section, we will upper-bound the term E[RT,1 ]. As mentioned in the main text, for simplicity, we will focus on the
case in which the underlying domain D is discrete, i.e., |D| < ∞. To begin with, we will need a supporting lemma showing
a uniform upper bound over the entire domain.
Lemma 1. Suppose that δ ∈ (0, 1) and βt , 2 log(|D|t2 π 2 /6δ). Then, with probability ≥ 1 − δ
1/2
|f (z) − µt−1 (z)| ≤ βt σt−1 (z) ∀z ∈ D, t ≥ 1 .
The proof of lemma 1 makes use of standard Gaussian tail bounds and a number of union bounds, and the proof is identical
to the proof of lemma 5.1 in (Srinivas et al., 2010). The next supporting lemma makes use of the Lipschitz continuity of f to
bound the differences between function values whose inputs only differ by the dimension corresponding to the number of
training epochs.
Lemma 2. Suppose that Assumption 2 holds and let δ 0 ∈ (0, 1). Then, with probability ≥ 1 − δ 0 ,
r
da
|f ([x, N ]) − f ([x, n])| ≤ N b log 0 ∀x, n = 1, . . . , N .
δ
Proof. Let z = [x, n] denote the input to the objective function f . Assumption 2, together with a union bound over
L 2
j = 1, . . . , d, implies that with probability ≥ 1 − dae−( b ) ,
|f (z) − f (z0 )| ≤ L||z − z0 ||1 ∀z ∈ D
Since [x, N ] and [x, n] differ only by the dimension corresponding to the number of training epochs, we have that
|f ([x, N ]) − f ([x, n])| ≤ LN

L 2
Then, the lemma follows by letting δ 0 = dae−( b ) .
The next lemma bounds E[rt,1 ] by the Gaussian process posterior standard deviation with some scaling constants.
Lemma 3. Let δ, δ 0 ∈ (0, 1) and κ ≥ 1 be the constant used in C2 in the BO-BOS algorithm. Then, at iteration t of the
BO-BOS algorithm, we have that, with probability ≥ 1 − δ − δ 0 ,
r
da
E[rt,1 ] ≤ 2κβt σt−1 ([xt , nt ]) + N b log 0 1nt <N .
1/2
δ
Proof. Firstly, with probability ≥ 1 − δ,
(1) (2) (3)

f (z∗ ) = f ([x∗ , N ]) ≤ µt−1 ([x∗ , N ]) + βt
1/2 1/2
σt−1 ([x∗ , N ]) ≤ µt−1 ([xt , N ]) + βt σt−1 ([xt , N ]) (10)
in which (1) follows from Assumption 1 which states that, for each x, the function value is monotonically non-decreasing in
the number of training epochs, which implies that at the (unknown) global maximum z∗ , the dimension corresponding to the
number of epochs is equal to N . (2)

√makes use of Lemma 1, whereas (3) is due to the way xt is selected in the algorithm,
i.e., xt = argmaxx µt−1 ([x, N ]) + βt σt−1 ([x, N ]). As a result, we have that with probability ≥ 1 − δ − δ 0
(1)
1/2
E[rt,1 ] = E[f (z∗ ) − f ([xt , nt ])] ≤ E[βt σt−1 ([xt , N ]) + µt−1 ([xt , N ]) − f ([xt , nt ])]
1/2
= E[βt σt−1 ([xt , N ]) + µt−1 ([xt , N ]) − f ([xt , N ]) + f ([xt , N ]) − f ([xt , nt ])]
(2)
1/2
≤ E[2βt σt−1 ([xt , N ])] + E[f ([xt , N ]) − f ([xt , nt ])]
r
(3) da (11)
≤ E[2βt σt−1 ([xt , N ])] + N b log 0 1nt <N
1/2
δ
r
(4) da
≤ 2βt σt−1 ([xt , N ]) + N b log 0 1nt <N
1/2
δ
r
(5) da
≤ 2κβt σt−1 ([xt , nt ]) + N b log 0 1nt <N
1/2
δ
in which (1) follows from Equation 10, and (2) results from Lemma 1 and the linearity of the expectation operator. 1nt <N
in (3) is the indicator function, which takes the value of 1 if the event nt < N is true and 0 otherwise. (3) is obtained
by analyzing two different cases separately: if nt = N (xt is not early-stopped),
q then E[f ([xq t , N ]) − f ([xt , nt ])] = 0;
if nt < N (xt is early-stopped), then E[f ([xt , N ]) − f ([xt , nt ])] ≤ E[N b log da da
δ 0 ] = N b log δ 0 with probability
0
≥ 1 − δ following Lemma 2. (4) is due to the fact that σt−1 ([xt , N ]) only depends on the observations up to step
∗
t − 1 and is not dependent on the probability Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt ). (5) follows from the design of the
algorithm; in particular, if nt < N , then κσt−1 ([xt , nt ]) ≥ σt−1 ([xt , N ]) is guaranteed by C2; otherwise, if nt = N , then
κσt−1 ([xt , nt ]) ≥ σt−1 ([xt , nt ]) = σt−1 ([xt , N ]) since κ ≥ 1.
PT
Subsequently, we can upper bound E[RT,1 ] = t=1 E[rt,1 ] by extensions of Lemma 5.3 and 5.4 from (Srinivas et al., 2010),
which are presented here for completeness. The following lemma connects the information gain about the objective function
with the posterior predictive variance, whose proof results from straightforward extension of Lemma 5.3 of (Srinivas et al.,
2010).
Lemma 4. Let yT be a set of observations of size T , and let fT be the corresponding function values. The information gain
about fT from observing yT is
T
1X
I(yT ; fT ) = log[1 + σ −2 σt−1
2
([xt , nt ])] .
2 t=1
Next, we use the following lemma to bound the sum of the first term of the expected instantaneous regret as given in Lemma
3.
8 2 2
Lemma 5. Let δ ∈ (0, 1), C1 , log(1+σ −2 ) , βt , 2 log(|D|t π /6δ), and γT , maxA∈D,|A|=T I(yA ; fA ) is the maximum
information gain about f from any subset of size T . Then,
T
1/2
X p p
2βt σt−1 ([xt , nt ]) ≤ T C1 βT I(yT ; fT ) ≤ T C1 βT γT .
t=1
Proof. Firstly, we have that
(1)
1/2
(2βt σt−1 ([xt , nt ]))2 = 4βt σt−1
2
([xt , nt ]) ≤ 4βT σ 2 [σ −2 σt−1
2
([xt , nt ])]
(2) σ −2
≤ 4βT σ 2 log[1 + σ −2 σt−1
2
([xt , nt ])] (12)
log(1 + σ −2 )
8 1
≤ βT −2
log[1 + σ −2 σt−1
2
([xt , nt ])]
log(1 + σ ) 2
σ −2
in which (1) holds since βt is monotonically increasing in t; (2) results from the fact that σ −2 x ≤ log(1+σ −2 ) log[1 + σ −2 x]
2
for x ∈ (0, 1], whereas 0 < σt−1 ([xt , nt ]) ≤ 1. Next, summing over t = 1, . . . , T , we get
T T
X 1/2 8 1X
(2βt σt−1 ([xt , nt ]))2 ≤ βT −2
log[1 + σ −2 σt−1
2
([xt , nt ])]
t=1
log(1 + σ ) 2 t=1
(13)
(1) 8 (2)
= βT −2
I(yT ; fT ) ≤ C1 βT γT
log(1 + σ )
in which (1) results from Lemma 4, and (2) follows from the definitions of C1 and γT . Next, making use of the Cauchy-
Schwarz inequality, we get
v
T u T
X 1/2
√ uX 1/2
2βt σt−1 ([xt , nt ]) ≤ T t (2βt σt−1 ([xt , nt ]))2
t=1 t=1
(14)
p
≤ C1 T βT γT
which completes the proof.
Next, putting everything together, we get the follow lemma on the upper bound on E[RT,1 ].
Lemma 6. Suppose that Assumptions 1 and 2 hold. Let δ, δ 0 ∈ (0, 1), C1 = 8/ log(1 + σ −2 ), βt = 2 log(|D|t2 π 2 /6δ),
κ ≥ 1 be the constant used in C2 in the BO-BOS algorithm, and γT = maxA∈D,|A|=T I(yA ; fA ). Let τT = t=1 1nt <N
PT
be the number of BO iterations in which early stopping happens from iterations 1 to T . Assume that f is a sample from a
GP, and y(z) = f (z) + ∀z ∈ D in which ∼ N (0, σ 2 ). Then, with probability ≥ 1 − δ − δ 0 ,
T r
X p da
E[RT,1 ] = E[rt,1 ] ≤ κ T C1 βT γT + N b log 0 τT ∀T ≥ 1 .
t=1
δ
Proof.
T (2) T r
da
1n <N ]
(1) X X 1/2
E[RT,1 ] = E[rt,1 ] ≤ [2κβt σt−1 ([xt , nt ]) + N b log
t=1 t=1
δ0 t
(3) T r
da
N b log 0 1nt <N
p X
≤ κ T C1 βT γT +
t=1
δ
(15)
r T
da X
1n <N
p
= κ T C1 βT γT + N b log 0
δ t=1 t
r
p da
= κ T C1 βT γT + N b log 0 τT
δ
in which (1) follows from the linearity of the expectation operator, (2) results from Lemma 3, and (3) follows from Lemma
5.
C.3. Upper Bound on E[RT,2 ]

In this section, we prove an upper bound on E[RT,2 ]. A few supporting lemmas will be presented and proved first. To begin
with, the next lemma derives the appropriate choice of the incumbent values used in the BOS problems in different iterations
of the BO-BOS algorithm.
Lemma 7. Let the objective function f be a sample from a GP and y(z) = f (z) + ∀z ∈ D in which ∼ N (0, σ 2 ). Let
δ 00 ∈ (0, 1). At iteration t > 1, define ft−1
∗ ∗
, maxt0 =1,...,t−1 f (zt0 ) and yt−1 , maxt0 =1,...,t−1 yt0 ; for iteration t = 1,
∗ ∗ 00
define f0 , 0 and y0 , 0. Then with probability ≥ 1 − δ ,
∗ ∗
ft−1 ≥ yt−1 − ξt ∀t ≥ 1
in which r
π 2 t2 (t − 1)
ξt = 2σ 2 log ∀t > 1
6δ 00
and ξ1 = 0.
Proof. The lemma trivially holds for t = 1. Assume we are at iteration t > 1 of the BO-BOS algorithm, and let
t0 ∈ {1, 2, . . . , t − 1}. Since yt0 = f (zt0 ) + , in which ∼ N (0, σ 2 ), we have that yt0 ∼ N (f (zt0 ), σ 2 ). Making use of
the upper deviation inequality for Gaussian distribution and the definition of ξt , we get
ξt 2 6δ 00
Pr[yt0 ≥ f (zt0 ) + ξt ] ≤ e− 2σ2 = (16)
π 2 t2 (t− 1)
Denote the event that {∃ t0 ∈ {1, 2, . . . , t − 1} s.t. yt0 ≥ f (zt0 ) + ξt } as At . Next, taking a union bound over the entire
observation history t0 ∈ {1, 2, . . . , t − 1}, we get
t−1
X
Pr[At ] ≤ Pr[yt0 ≥ f (zt0 ) + ξt ]
t0 =1 (17)
6δ 00 6δ 00
≤(t − 1) 2 2 = 2 2
π t (t − 1) π t
00
which implies that at iteration t, with probability ≥ 1 − π6δ2 t2 , yt0 − f (zt0 ) < ξt ∀t0 ∈ {1, 2, . . . , t − 1}, which further
∗ ∗
suggests that yt−1 − ft−1 ≤ ξt at iteration t. Next, taking a union bound over t ≥ 1, we get
X X 6δ 00
Pr[∃t ≥ 1 s.t. At holds] ≤ Pr[At ] ≤ = δ 00 (18)
π 2 t2
t≥1 t≥1
which suggests that, with probability ≥ 1 − δ 00 , yt−1

∗ ∗
− ft−1 ≤ ξt ∀t ≥ 1, and thus completes the proof.
The next lemma shows that, with appropriate choices of the incumbent value, the posterior probability used in Bayesian
optimal stopping is upper-bounded.
∗
Lemma 8. If in iteration t of the BO-BOS algorithm, the BOS algorithm is run with the incumbent value yt−1 − γt and the
corresponding cost parameters K1 , K2 and cd0 , and the algorithm early-stops after nt < N epochs, then with probability
≥ 1 − δ 00 ,
∗ K2 + cd0
Pr(f ([xt , N ]) > ft−1 |yt,nt ) ≤ ∀t ≥ 1 . (19)
K1
Proof. Recall that when running the Bayesian optimal stopping algorithm in iteration t of BO-BOS, we only early-stop the
experiment (nt < N ) when we can safely conclude that the performance of the currently evaluated hyperparameter xt will
end up having smaller (or equal) validation accuracy than the currently observed optimum offset by a noise correction term:
∗
yt−1 − ξt ; i.e., when the expected loss of decision d1 is the smallest among all decisions. Therefore, when the evaluation of
xt is early-stopped after nt < N epochs, we can conclude that
∗
K1 Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt )
h
∗ ∗
≤Eyt,nt +1 |yt,nt min{K1 Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt +1 ), K2 Pr(f ([xt , N ]) ≤ yt−1 − ξt |yt,nt +1 ),
i
Eyt,nt +2 |yt,nt +1 [ρt,nt +2 (yt,nt +2 )] + cd0 } + cd0
(20)
∗
≤Eyt,nt +1 |yt,nt [K2 Pr(f ([xt , N ]) ≤ yt−1 − ξt |yt,nt +1 )] + cd0
∗
≤K2 Eyt,nt +1 |yt,nt [Pr(f ([xt , N ]) ≤ yt−1 − ξt |yt,nt +1 )] + cd0
≤K2 + cd0
Equation 20, together with Lemma 7, implies that
∗ ∗
Pr(f ([xt , N ]) > ft−1 |yt,nt ) ≤Pr(f ([xt , N ]) > yt−1 − ξt |yt,nt )
K2 + cd0 (21)
≤
K1
which holds uniformly for all t ≥ 1 with probability ≥ 1 − δ 00 .
Subsequently, we use the next Lemma to upper-bound E[RT2 ] by the BOS cost parameters. We set K2 and cd0 as constants,
and use different values of K1 in different iterations t of the BO-BOS algorithm, which is represented by K1,t .
K2 +cd0
Lemma 9. In iteration t of the BO-BOS algorithm, define K1,t , ηt . Then, with probability ≥ 1 − δ 00 ,
T
X
E[RT,2 ] ≤ ηt ∀T ≥ 1 .
t=1
Proof. Recall that according to Assumption 1, the value of the objective function f is bounded in the range [0, 1]. In iteration
t, assume we early-stop the evaluation of xt after nt < N epochs, then
(1)
∗
E[f ([xt , N ]) − ft−1 |yt,nt ] ≤ E[1f ([xt ,N ])−ft−1
∗
∗
>0 |yt,nt ] = Pr(f ([xt , N ]) > ft−1 |yt,nt )
(22)
Step (1) in Equation 22 is because x ≤ 1x>0 ∀x ∈ [−1, 1] and substituting x = f ([xt , N ]) − ft−1
∗
. As a result, with
probability ≥ 1 − δ 00
(1) X (2) X (3) X

∗ ∗
E[RT,2 ] = E[rt,2 ] = E[f ([xt , N ]) − ft−1 |yt,nt ] ≤ Pr(f ([xt , N ]) > ft−1 |yt,nt )
{t|nt <N } {t|nt <N } {t|nt <N }
T
(23)
(4) X X
≤ ηt ≤ ηt
{t|nt <N } t=1
in which (1) follows from the linearity of expectation, (2) holds because the Expectation of rt,2 is taken over Pr(f ([xt , N ]) >
∗
yt−1 − ξt |yt,nt ), (3) results from Equation 22, and (4) follows from Lemma 8. This completes the proof.
C.4. Putting Things Together

In this section, we put everything from the previous two sections together to prove the main theorems.
C.4.1. P ROOF OF T HEOREM 1

RT
Theorem 1 can be proven by combining Lemmas 6 and 9, and making use of the fact that ST ≤ T .
C.4.2. P ROOF OF T HEOREM 2

Below we analyze the asymptotic behavior of each of the three terms in the upper bound of E[ST ] in Theorem 1, which is
re-presented here for ease of reference.
√ PT r
κ T C1 βT γT t=1 ηt 1 da
E[ST ] ≤ + + N b log 0 τT . (24)
T T T δ
The first term in the upper bound of E[ST ] Firstly, the first term in the upper bound matches the upper bound on the
simple regret of the GP-UCB algorithm (Srinivas et al., 2010) (up to the constant κ). The maximum information gain,
γT , has been analyzed for a few of the commonly used kernels in GP (Srinivas et al., 2010). For example, for the Square
Exponential kernel, γT = O((log T )d+1 ), whereas for the Matérn kernel with ν > 1, γT = O(T d(d+1)/(2ν+d(d+1)) log T ).
Plugging both expressions of γT into Theorem √
1, together with the expression of βT as given in Theorem 1, shows that both
kernels lead to sub-linear growth of the term T C1 βT γT , which implies that the first term in the upper bound of E[ST ]
asymptotically goes to 0.
The second term in the upper bound of E[ST ] Given that K1,t is an increasing sequence with K1,1 ≥ K2 + cd0 , the
PT PT K2 +cd0
series t=1 ηt = t=1 K 1,t
grows sub-linearly, thus making the second term in the upper bound of E[ST ] given in
PT
ηt
Theorem 1, t=1
T , asymptotically go to 0.
The third term in the upper bound of E[ST ] Next, suppose that K1,t becomes +∞ for the first time at iteration T0 .
Since K1,t is a non-decreasing sequence, K1,t = +∞ for all t ≥ T0 . Therefore, for t ≥ T0 , decision d1 will never be taken
and the algorithm will never early-stop. In other words, nt = N for all t ≥ T0 .
Therefore, we can conclude that τT ≤ T0 for all T ≥ 1. As a result, the last term in the upper bound on E[ST ] in Theorem 1
can be upper-bounded by
q
T N b log da
r
τT da 0 δ0 1
N b log 0 ≤ =O (25)
T δ T T
which asymptotically goes to 0 as T goes to +∞, because the numerator term is a constant. Therefore, this term also
asymptotically vanishes in the upper bound.
To summarize, if the BOS parameters are selected according to Theorem 2, we have that
√ PT !
T βT γT t=1 ηt 1
E[ST ] = O + + (26)
T T T
and E[ST ] goes to zero asymptotically.
D. Additional Experimental Details

In each experiment, the same initializations (6 initial points if not further specified) are used for all BO-based methods:
GP-UCB, BOCA, LC Prediction, and BO-BOS. The Square Exponential kernel is used for BOCA since the algorithm
is only given for this kernel (Kandasamy et al., 2017), the other BO-based algorithms use the Matérn kernel; the kernel
hyperparameters are updated by maximizing the Gaussian process marginal likelihood after every 10 BO iterations. In the
BO-BOS algorithm, since the number of training epochs is an input to the GP surrogate function, some of the intermediate
observations (n < N ) can be used as additional input to GP to improve the modeling of the objective function. However,
using the observation after every epoch as input leads to poor scalability. Therefore, for all experiments with N = 50 (which
include most of the experiments), we use the observations after first, 10-th, 20-th, 30-th and 40-th epochs as additional
inputs to the GP surrogate function; whereas for the RL experiment with N = 100 in section 5.3.1, we use the 1-th, 20-th,
40-th, 60-th and 80-th intermediate observations as additional inputs. 100,000 forward simulation samples are used for
each BOS algorithm; the grid size of the discretized summary statistics is set to 100; for simplicity, the incumbent value
∗
at iteration t is chosen as yt−1 = maxt0 =1,...,t−1 yt0 , thus ignoring the observation noise. In the LC Prediction algorithm
(Domhan et al., 2015), learning curve prediction is performed after every 2 epochs. In Hyperband (Li et al., 2017), the
successive halving parameter η is set to 3 as recommended by the original authors, and the maximum number of epochs
is set to N = 80 (we observed that setting N = 80 led to better performance than N = 50 since it allows the Hyperband
algorithm to run for more epochs overall).
D.1. Hyperparameter Tuning for Logistic Regression

In the first set of experiments, we perform hyperparameter tuning for a simple ML model, logistic regression (LR). The
LR model is trained using the MNIST image dataset, which consists of 70,000 images of the 10 digits, corresponding to a
10-class classification problem. Three hyperparameters are tuned: the batch size (20 to 500), L2 regularization parameter
(10−6 to 1.0), and learning rate (10−3 to 0.1). We use 80% of the images as the training set and the remaining 20% as the
validation set.
Some of the learning curves during a particular run of the BO-BOS algorithm is shown in Fig. 7. It can be observed
that the learning curves that show minimal potential in achieving small validation errors are early-stopped, whereas the
promising hyperparameter settings are run for larger number of epochs. The reliability of the early stopping achieved by the
BO-BOS algorithm is demonstrated in Fig. 8. In this figure, the green triangles correspond to the learning curves that are not
early-stopped (nt = N ), and the red circles represent the final validation errors (after training for the maximum number of
epochs N ) that could have been reached by the early-stopped learning curves (nt < N ). Note that the red circles are shown
only for the purpose of illustration and are not observed in practice. As displayed in the figure, the early stopping decisions
made during the BO-BOS algorithm are reliable, since those early-stopped learning curves all end up having large validation
errors.
0.35
0.30
validation error
0.25
0.20
0.15
0.10
0 10 20 30 40 50
epochs
Figure 7. Some learning curves during the BO-BOS algorithm.
0.30 Validation Error

Final Value of Early-Stopped Iterations
Fully-Evaluated Iterations
0.25
Validation Error
0.20
0.15
0.10
0.05
0 5 10 15 20 25
BO iterations
Figure 8. Illustration of the effectiveness of the early stopping decisions made during the BO-BOS algorithm.
D.2. Hyperparameter Tuning for Convolutional Neural Networks

Next, we tune the hyperparameters of convolutional neural networks (CNN) using the CIFAR-10 (Krizhevsky, 2009) and
Street View House Numbers (SVHN) dataset (Netzer et al., 2011). Both tasks correspond to 10-class classification problems.
For CIFAR-10, 50, 000 images are used as the training set and 10,000 images are used as the validation set; for SVHN,
73,257 and 26032 images are used as the training and validation sets respectively following the original dataset partition. The
CNN model consists of three convolutional layers (each followed by a max-pooling layer) followed by one fully-connected
layer. We tune six hyperparameters in both experiments: the batch size (32 to 512), learning rate (10−7 to 0.1), learning
rate decay (10−7 to 10−3 ), L2 regularization parameter (10−7 to 10−3 ), the number of convolutional filters in each layer
(128 to 256), and the number of units in the fully-connected layer (256 to 512). In addition to the results in Figure 2, the
corresponding figures with standard error is presented below in Figures 9 and 10, which demonstrate the robustness of the
performance advantages of BO-BOS.
0.32
0.31
0.30
0.29
Validation Error
0.28
0.27
0.26 GP-UCB
Hyperband
0.25 BOCA
LC Prediction
0.24 BO-BOS
0 1 2 3 4 5 6 7 8
Wall-Clock Time (Hours)
Figure 9. Best-found validation error of CNN v.s. run-time using the CIFAR-10 dataset, with standard error (averaged over 30 random
initializations).
0.16
GP-UCB
0.15 Hyperband
BOCA
0.14 LC Prediction
BO-BOS
Validation Error
0.13
0.12
0.11
0.10
0.09
0 1 2 3 4 5 6 7 8
Figure 10. Best-found validation error of CNN v.s. run-time using the SVHN dataset, with standard error (averaged over 30 random
initializations).
D.3. Policy Search for Reinforcement Learning

We apply our algorithm to a continuous control task: the Swimmer-v2 environment from OpenAI Gym, MuJoCo (Brockman
et al., 2016; Todorov et al., 2012). The task involves controlling two joints of a swimming robot to make it swim forward
as fast as possible. The state of the robot is represented by an 8-dimensional feature vector, and the action space is
2-dimensional corresponding to the two joints. We use a linear policy, in which the policy is represented by an 8 × 2 matrix
that maps each state vector to the corresponding action vector. In this setting, the input parameters, x, to the GP-UCB and
BO-BOS algorithms are the 16 parameters of the policy matrix, and the objective function is the discounted cumulative
rewards in an episode. Each episode of the task consists of 1,000 steps. We set N (the maximum number of epochs) to
be smaller than 1, 000 by treating a fixed number of consecutive steps as one single epoch. E.g., we can set N = 50 or
N = 100 by treating every 20 or 10 consecutive steps as one epoch respectively. The rewards are clipped, scaled, and
normalized such that the discounted cumulative rewards of each episode is bounded in the range [0, 1]; for each evaluated
policy, we also record the un-discounted and un-scaled cumulative rewards, which are the ultimate objective to be maximized
and reported in Fig. 3 in the main text. Each policy evaluation consists of running 5 independent episodes with the given
policy, and returning the average discounted cumulative rewards, i.e., average return, as the observed function value.
As mentioned in the main text, the rewards are discounted in order to make the objective function, the discounted cumulative
rewards, resemble the learning curves of ML models, such that the BO-BOS algorithm can be naturally applied. This
rationale is illustrated in Fig. 11, which plots some example un-discounted (γ = 1.0) and discounted (γ = 0.9) cumulative
rewards respectively. The figures indicate that, compared with un-discounted cumulative rewards, discounted cumulative
rewards bear significantly closer resemblance to the learning curves of ML models, thus supporting the claim made in the
main text motivating the use of discounted rewards, as well as the experimental results shown in Fig. 3 (specifically, the
poor performance of the curve corresponding to N = 50 and γ = 1.0). In addition to the results presented in the main text
in section 5.3.1, we further present the results with standard errors in Fig. 12, to emphasize the significant performance
advantage offered by BO-BOS compared with GP-UCB. To avoid clutter, we only present the results with error bar for
GP-UCB with γ = 1.0 and BO-BOS with N = 50 and γ = 0.9, which are best-performing settings for GP-UCB and
BO-BOS respectively.
0.65
Un-discounted Cumulative Reward
0.60
0.55
0.50
0.45
0.40
0 10 20 30 40 50
Epochs
(a) Un-discounted (γ = 1.0).
0.85
0.80
Discounted Cumulative Reward
0.75
0.70
0.65
0.60
0.55
0.50
0.45
0 10 20 30 40 50
Epochs
(b) Discounted (γ = 0.9).
Figure 11. Example curves of un-discounted and discounted cumulative rewards
D.4. Joint Hyperparameter Tuning and Feature Selection

In this set of experiments, we use the gradient boosting model (XGBoost (Chen & Guestrin, 2016)), tuning four hyperparam-
eters: the learning rate (10−3 to 0.5), maximum depth of each decision tree (2 to 15), feature sub-sampling ratio for each tree
(0.3 to 1.0), and L1 regularization parameter (0.0 to 5.0). We use the email spam dataset from the UCI Machine Learning
200
Average Return
150 GP-UCB ( = 1.0)
BO-BOS (N = 50, = 0.9)
100
50
0 1 2 3 4 5 6
Steps 1e5
Figure 12. Best-found return (averaged over 5 episodes) v.s. the total number of steps of the robot in the environment (averaged over 30
random initializations) using the Swimmer-v2 task, with standard error.
0.038
GP-UCB
Hyperband
0.036 BO-BOS
Validation Error
0.034
0.032
0.030
0.028
0 1 2 3 4 5 6
Figure 13. Best-found validation error of XGBoost v.s. run-time with standard error (averaged over 30 random initializations), obtained
using joint hyperparameter tuning and feature selection.
Repository (Dheeru & Karra Taniskidou, 2017), which represents a binary classification problem: whether the email is a
spam or not. We use 3065 emails as the training set and the remaining 1536 emails as the validation set; each email consists
of 57 features. The maximum number of features for each hyperparameter setting is set as N = 50. In addition to Figure 4
in the main text, the same plot with error bar (standard error) is shown in Figure 13.

Bayesian Optimization Meets Bayesian Optimal Stopping

Uploaded by

Document Informationclick to expand document information

Document Informationclick to expand document information

Copyright:

Available Formats

Bayesian Optimization Meets Bayesian Optimal Stopping

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Bayesian Optimization Meets Bayesian Optimal Stopping

Uploaded by

Copyright:

Available Formats

Bayesian Optimization Meets Bayesian Optimal Stopping

Zhongxiang Dai 1 Haibin Yu 1 Bryan Kian Hsiang Low 1 Patrick Jaillet 2

Abstract Silver et al., 2016). However, a major difficulty faced by ML

4. Theoretical Analysis The proof of both theorems is presented in Appendix C.

algorithm early-stops very often, thus leading to smaller

is usually estimated by running the agent in the environment

can be naturally solved by BO-BOS by simply replacing the

the considerable potential of BO-based policy search algo- 6. Conclusion

Acknowledgements Chen, T. and Guestrin, C. XGBoost: A scalable tree boost-

Falkner, S., Klein, A., and Hutter, F. BOHB: Robust and

Powell, W. B. and Ryzhov, I. O. Optimal learning. John

A. Approximate Backward Induction for Bayesian Optimal Stopping

B. Approximate Backward Induction Algorithm for Solving BOS Problem in BO-BOS

C. Proof of Theorems 1 and 2

C.1. Regret Decomposition

rt =f (z∗ ) − ft∗ = f (z∗ ) − 0 max f (zt0 )

of RT , therefore, sub-linear growth of RT0 implies that RT also grows sub-linearly.

rt− = f (z∗ ) − ft−1

in which (1) results from Assumption 1. As a result, RT0 can be re-written as

prove an upper bound on the expected value of RT0 ,

E[RT0 ] ≤ E[RT,1 ] + E[RT,2 ] (9)

C.2. Upper Bound on E[RT,1 ]

|f (z) − f (z0 )| ≤ L||z − z0 ||1 ∀z ∈ D

|f ([x, N ]) − f ([x, n])| ≤ LN

Proof. Firstly, with probability ≥ 1 − δ,

(1) (2) (3)

number of epochs is equal to N . (2)

information gain about f from any subset of size T . Then,

Proof. Firstly, we have that

which completes the proof.

C.3. Upper Bound on E[RT,2 ]

which suggests that, with probability ≥ 1 − δ 00 , yt−1

which holds uniformly for all t ≥ 1 with probability ≥ 1 − δ 00 .

(1) X (2) X (3) X

C.4. Putting Things Together

C.4.1. P ROOF OF T HEOREM 1

C.4.2. P ROOF OF T HEOREM 2

and E[ST ] goes to zero asymptotically.

D. Additional Experimental Details

D.1. Hyperparameter Tuning for Logistic Regression

Figure 7. Some learning curves during the BO-BOS algorithm.

0.30 Validation Error

D.2. Hyperparameter Tuning for Convolutional Neural Networks

D.3. Policy Search for Reinforcement Learning

(a) Un-discounted (γ = 1.0).

(b) Discounted (γ = 0.9).

Figure 11. Example curves of un-discounted and discounted cumulative rewards

D.4. Joint Hyperparameter Tuning and Feature Selection

You might also like