BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization

Introduction

Computational modeling and machine learning (ML) have led to an acceleration of scientific innovation in diverse areas, ranging from drug design to robotics to material science. These tasks often involve solving time- and resource-intensive global optimization problems to achieve optimal performance. Bayesian optimization (BO) , an established methodology for sample-efficient sequential optimization, has been proposed as an effective solution to such problems, and has been applied successfully to tasks ranging from hyperparameter optimization , robotic control , chemical design , and tuning and policy search for internet-scale software systems . Meanwhile, ML research has been undergoing a revolution driven largely by new programming frameworks and hardware that reduce the time from ideation to execution . While BO has become rich with new methodologies, today there is no coherent framework that leverages these computational advances to simplify and accelerate BO research in the same way that modern frameworks have for deep learning. In this paper, we address this gap by introducing BoTorch, a modular and scalable Monte Carlo (MC) framework for BO that is built around modern paradigms of computation, and theoretically grounded in novel convergence results. Our contributions include:

A novel approach to optimizing MC acquisition functions that effectively combines with deterministic higher-order optimization algorithms and variance reduction techniques.

The first convergence results for sample average approximation (SAA) of MC acquisition functions, including novel general convergence results for SAA via randomized quasi-MC.

A new, SAA-based “one-shot” formulation of the Knowledge Gradient, a look-ahead acquisition function, with improved performance over the state-of-the-art.

Composable model-agnostic abstractions for MC BO that leverage modern computational technologies, including auto-differentiation and scalable parallel computation on CPUs and GPUs.

We discuss related work in Section 2 and then present the methodology underlying BoTorch in Sections 3 and 4. Details of the BoTorch framework, including its modular abstractions and implementation examples, are given in Section 5. Numerical results are provided in Section 6.

Background and Related Work

Popular libraries for BO include Spearmint , GPyOpt , Cornell-MOE , RoBO , Emukit , and Dragonfly . We provide further discussion of these packages in Appendix A. Two other libraries, ProBO and GPFlowOpt , are of particular relevance. ProBO is a recently suggested frameworkNo implementation of ProBO is available at the time of this writing. for using general probabilistic programming in BO. While its model-agnostic approach is similar to ours, ProBO, unlike BoTorch, does not benefit from gradient-based optimization provided by differentiable programming, or algebraic methods designed to exploit GPU acceleration. GPFlowOpt inherits support for auto-differentiation and hardware acceleration from TensorFlow [via GPFlow, 64], but unlike BoTorch, it does not use algorithms designed to specifically exploit this potential. Neither ProBO nor GPFlowOpt naturally support MC acquisition functions. In contrast to all existing libraries, BoTorch is a modular programming framework and employs novel algorithmic approaches that achieve a high degree of flexibility and performance.

The MC approach to optimizing acquisition functions has been considered in the BO literature to an extent, typically using stochastic methods for optimization . Our work takes the distinctive view of sample average approximation (SAA), an approach that combines sampling with deterministic optimization and variance reduction techniques. To our knowledge, we provide the first theoretical analysis and systematic implementation of this approach in the BO setting.

Monte-Carlo Acquisition Functions

The next step in BO is to optimize an acquisition function evaluated on $f_{\mathcal{D}}(\mathbf{x})$ over the candidate set $\mathbf{x}$. Following , many acquisition functions can be written as

In some situations, the expectation over $f_{\mathcal{D}}(\mathbf{x})$ in (1) and its gradient $\nabla_{\mathbf{x}}\alpha(\mathbf{x};\Phi,\mathcal{D})$ can be computed analytically, e.g. if one considers a single-output ($m\!=\!1$) model, a single candidate ($q\!=\!1$) point $x$, a Gaussian posterior $f_{\mathcal{D}}(x)=\mathcal{N}(\mu_{x},\sigma_{x}^{2})$, and the identity objective $g(f)\equiv f$. Expected Improvement (EI) is a popular acquisition function that maximizes the expected difference between the currently observed best value $f^{*}$ (assuming noiseless observations) and $f$ at the next query point, through the utility $a(f,f^{*})=\max(f-f^{*},0)$. EI and its gradient have a well-known analytic form .

In general, analytic expressions are not available for arbitrary objective functions $g(\cdot)$, utility functions $a(\cdot\,,\cdot)$, non-Gaussian model posteriors, or collections of points $\mathbf{x}$ which are to be evaluated in a parallel or asynchronous fashion . Instead, MC integration can be used to approximate the expectation (1) using samples from the posterior. An MC approximation $\hat{\alpha}_{\!N}(\mathbf{x};\Phi,\mathcal{D})$ of (1) using $N$ samples $\xi_{\mathcal{D}}^{i}(\mathbf{x})\sim f_{\mathcal{D}}(\mathbf{x})$ is straightforward:

The obvious way to evaluate (2) is to draw i.i.d. samples $\xi_{\mathcal{D}}^{i}(\mathbf{x})$. Alternatively, randomized quasi-Monte Carlo (RQMC) techniques can be used to significantly reduce the variance of the estimate and its gradient (see Appendix E for additional details).

MC Bayesian Optimization via Sample Average Approximation

To generate a new candidate set $\mathbf{x}$, one must optimize the acquisition function $\alpha$. Doing this effectively, especially in higher dimensions, typically requires using gradient information. For differentiable analytic acquisition functions (e.g. EI, UCB), one can either manually implement gradients, or use auto-differentiation to compute $\nabla_{\!x}\alpha(x;\Phi,\mathcal{D})$, provided one can differentiate through the posterior parameters (as is the case for Gaussian posteriors).

An unbiased estimate of the MC acquisition function gradient $\nabla_{\!\mathbf{x}}\alpha(\mathbf{x};\Phi,\mathcal{D})$ can often be obtained from (2) via the reparameterization trick . The basic idea is that $\xi\sim f_{\mathcal{D}}(\mathbf{x})$ can be expressed as a suitable (differentiable) deterministic transformation $\xi=h_{\mathcal{D}}(\mathbf{x},\epsilon)$ of an auxiliary random variable $\epsilon$ independent of $\mathbf{x}$. For instance, if $f_{\mathcal{D}}(\mathbf{x})\sim\mathcal{N}(\mu_{\mathbf{x}},\Sigma_{\mathbf{x}})$, then $h_{\mathcal{D}}(\mathbf{x},\epsilon)=\mu_{\mathbf{x}}+L_{\mathbf{x}}\epsilon$, with $\epsilon\sim\mathcal{N}(0,I)$ and $L_{\mathbf{x}}L_{\mathbf{x}}^{T}=\Sigma_{\mathbf{x}}$. If $a(\cdot,\Phi)$ and $g(\cdot)$ are differentiable, then $\nabla_{\!\mathbf{x}}a(g(\xi),\Phi)=\nabla_{\!g}a\nabla_{\!\xi}g\nabla_{\!\mathbf{x}}h_{\mathcal{D}}(\mathbf{x},\epsilon)$.

Our primary methodological contribution is to take a sample average approximation approach to BO. The conventional way of optimizing MC acquisition functions of the form (2) is to re-draw samples from $\epsilon$ for each evaluation and apply stochastic first-order methods such as Stochastic Gradient Descent (SGD) . In our SAA approach, rather than re-drawing samples from $\epsilon$ for each evaluation of the acquisition function, we draw a set of base samples $E:=\{\epsilon^{i}\}_{i=1}^{N}$ once, and hold it fixed between evaluations throughout the course of optimization (this can be seen as a specific incarnation of the method of common random numbers). Conditioned on $E$, the resulting MC estimate $\hat{\alpha}_{\!N}(\mathbf{x};\Phi,\mathcal{D})$ is deterministic. We then obtain the candidate set $\hat{\mathbf{x}}_{\!N}^{*}$ as

The gradient $\nabla_{\mathbf{x}}\hat{\alpha}_{\!N}(\mathbf{x};\Phi,\mathcal{D})$ can be computed as the average of the sample-level gradients, exploiting auto-differentiation. We emphasize that whether this average is a “proper” (i.e., unbiased, consistent) estimator of $\nabla_{\mathbf{x}}\alpha(\mathbf{x};\Phi,\mathcal{D})$ is irrelevant for the convergence results we will derive below.

Under relatively weak conditions,Many utility functions $a$ are Lipschitz, including those representing (parallel) EI and UCB . Lipschitzness is a sufficient condition, and convergence can also be shown in less restrictive settings (see Appendix D). Theorem 1 ensures not only that the optimizer $\hat{\mathbf{x}}_{\!N}^{*}$ of $\hat{\alpha}_{\!N}$ converges to an optimizer of the true $\alpha$ with probability one, but also that the convergence (in probability) happens at an exponential rate. We stated Theorem 1 informally and for i.i.d. base samples for simplicity. In Appendix D.3 we give a formal statement, and extend it to base samples generated by a family of RQMC methods, leveraging recent theoretical advances . While at this point we do not characterize improvements in theoretical convergence rates of RQMC over MC for SAA, we observe empirically that RQMC methods work remarkably well in practice (see Figures 3 and 3).

The primary benefit from SAA comes from the fact that in order to optimize $\hat{\alpha}_{\!N}(\mathbf{x};\Phi,\mathcal{D})$ for fixed base samples $E$, one can now employ the full toolbox of deterministic optimization, including quasi-Newton methods that provide faster convergence speeds and are generally less sensitive to optimization hyperparameters than stochastic first-order methods. By default, we use multi-start optimization via L-BFGS-B in conjunction with an initialization heuristic that exploits fast batch evaluation of acquisition functions (see Appendix F.1). We find that in practice the bias from using SAA only has a minor effect on the performance relative to using the analytic ground truth, and often improves performance relative to stochastic approaches (see Appendix E), while avoiding tedious tuning of optimization hyperparameters such as learning rates.

2 One-Shot Formulation of the Knowledge Gradient using SAA

The acquisition functions mentioned above, such as EI and UCB, are myopic, that is, they do not take into account the effect of an observation on the model in future iterations. In contrast, look-ahead methods do. Our SAA approach enables a novel formulation of a class of look-ahead acquisition functions. For the purpose of this paper we focus on the Knowledge Gradient (KG) , but our methods extend to other look-ahead acquisition functions such as two-step EI .

KG quantifies the expected increase in the maximum of $f$ from obtaining the additional (random) observation data $\{\mathbf{x},y_{\mathcal{D}}(\mathbf{x})\}$. KG often shows improved BO performance relative to simpler, myopic acquisition functions such as EI , but in its traditional form it is computationally expensive and hard to implement, two challenges that we address in this work. Writing $\mathcal{D}_{\mathbf{x}}:=\mathcal{D}\cup\{\mathbf{x},\mathbf{y}_{\mathcal{D}}(\mathbf{x})\}$, we introduce a generalized variant of parallel KG (qKG) :

Conditional on the fixed base samples, (5) does not exhibit the nested structure used in the conventional formulation (which requires solving an optimization problem to get a noisy gradient estimate). Moving the maximization outside of the sample average yields the equivalent problem

Programmable Bayesian Optimization with BoTorch

SAA provides an efficient and robust approach to optimizing MC acquisition functions through the use of deterministic gradient-based optimization. In this section, we introduce BoTorch, a complementary differentiable programming framework for Bayesian optimization research. Following the conceptual framework outlined in Figure 1, BoTorch provides modular abstractions for representing and implementing sophisticated BO procedures. Operations are implemented as PyTorch modules that are highly parallelizable on modern hardware and end-to-end differentiable, which allows for efficient optimization of acquisition functions. Since the chain of evaluations on the sample level does not make any assumptions about the form of the posterior, BoTorch’s primitives can be directly used with any model from which re-parameterized posterior samples can be drawn, including probabilistic programs , Bayesian neural networks , and more general types of GPs . In this paper, we focus on an efficient and scalable implementation of GPs, GPyTorch .

To illustrate the core components of BoTorch, we demonstrate how both known and novel acquisition functions can readily be implemented. For the purposes of exposition, we show a set of simplified implementations here; details and additional examples are given in Appendices G and H.

In our first example, we consider $q$ParEGO , a variant of ParEGO , a method for multi-objective optimziation.

Code Example LABEL:codex:moo implements the inner loop of $q$ParEGO. We begin by instantiating a GenericMCObjective module that defines an augmented Chebyshev scalarization. This is an instance of BoTorch’s abstract MCObjective, which applies a transformation $g(\cdot)$ to samples $\xi$ from a posterior in its forward($\xi$) pass. In line 5, we instantiate an MCAcquisitionFunction module, in this case, qExpectedImprovement, parallel EI. Acquisition functions combine a model and the objective into a single module that assigns a utility $\alpha(\mathbf{x})$ to a candidate set $\mathbf{x}$ in its forward pass. Models can be any PyTorch module implementing a probabilistic model conforming to BoTorch’s basic Model API. Finally, candidate points are selected by optimizing the acquisition function, through the use of the optimize_acqf() utility function, which finds the candidates $\mathbf{x}^{*}\in\operatorname*{arg\,max}_{\mathbf{x}}\alpha(\mathbf{x})$. Auto-differentiation makes it straightforward to use gradient-based optimization even for complex acquisition functions and objectives. Our SAA approach permits the use of deterministic higher-order optimization to efficiently and reliably find $\mathbf{x}^{*}$.

In it is shown how performing operations on independently modeled objectives yields better optimization performance when compared to modeling combined outcomes directly (e.g., for the case of calibrating the outputs of a simulator). MCObjective is a powerful abstraction that makes this straightforward. It can also be used to implement unknown (i.e. modeled) outcome constraints: BoTorch implements a ConstrainedMCObjective to compute a feasibility-weighted objective using a sample-level differentiable relaxation of the feasibility .

2 Implementing Parallel, Asynchronous Noisy Expected Improvement

Code Example LABEL:codex:Abstractions:ImplementationExamples:NEI:qNEI provides an implementation of qNEI as formulated in (7). New MC acquisition functions are defined by extending an MCAcquisitionFunction base class and defining a forward pass that compute the utility of a candidate $\mathbf{x}$. In the constructor (not shown), the programmer sets X_baseline to an appropriate subset of the points at which the function was observed.

3 Look-ahead Bayesian Optimization with One-Shot KG

Code Example LABEL:codex:Abstractions:ImplementationExamples:OKG shows a simplified OKG implementation, as discussed in Section 4.2.

Experiments

BoTorch utilizes inference and optimization methods designed to exploit parallelization via batched computation, and integrates closely with GPyTorch . These model have fast test-time (predictive) distributions and sampling. This is crucial for BO, where the same models are evaluated many times in order to optimize the acquisition function. GPyTorch makes use of structure-exploiting algebra and local interpolation for $\mathcal{O}(1)$ computations in querying the predictive distribution, and $\mathcal{O}(T)$ for drawing a posterior sample at $T$ points, compared to the standard $\mathcal{O}(n^{2})$ and $\mathcal{O}(T^{3}n^{3})$ computations .

Figure 5 reports wall times for batch evaluation of qExpectedImprovement at multiple candidate sets $\{\mathbf{x}^{i}\}_{i=1}^{b}$ for different MC samples sizes $N$, on both CPU and GPU for a GPyTorch GP. We observe significant speedups from running on the GPU, with scaling essentially linear in the batch size $b$, except for very large $b$ and $N$. Figure 5 shows between 10–40X speedups when using fast predictive covariance estimates over standard posterior inference in the same setting. The speedups grow slower on the GPU, whose cores do not saturate as quickly as on the CPU when doing standard posterior inference (for additional details see Appendix B). Together, batch evaluation and fast predictive distributions enable efficient, parallelized acquisition function evaluation for a very large number (tens of thousands) of points. This scalability allows us to implement and exploit novel highly parallelized initialization and optimization techniques.

2 Bayesian Optimization Performance Comparisons

We compare (i) the empirical performance of standard algorithms implemented in BoTorch with those from other popular BO libraries, and (ii) our novel acquisition function, OKG, against other acquisition functions, both within BoTorch and in other packages. We isolate three key frameworks—GPyOpt, Cornell MOE (MOE EI, MOE KG), and Dragonfly—because they are the most popular libraries with ongoing supportWe were unable to install GPFlowOpt due to its incompatibility with current versions of GPFlow/TensorFlow. and are most closely related to BoTorch in terms of state-of-the-art acquisition functions. GPyOpt uses an extension of EI with a local penalization heuristic (henceforth GPyOpt LP-EI) for parallel optimization . For Dragonfly, we consider its default ensemble heuristic (henceforth Dragonfly GP Bandit) .

Our results provide three main takeaways. First, we find that BoTorch’s algorithms tend to achieve greater sample efficiency compared to those of other packages (all packages use their default models and settings). Second, we find that OKG often outperforms all other acquisition functions. Finally, OKG is more computationally scalable than MOE KG (the gold-standard implementation of KG), showing significant reductions in wall time (up to 6X, see Appendix C.2) while simultaneously achieving improved optimization performance (Figure 7).

Synthetic Test Functions: We consider BO for parallel optimization of $q=4$ design points, on four noisy synthetic functions used in Wang et al. : Branin, Rosenbrock, Ackley, and Hartmann. Figure 7 reports means and 95% confidence intervals over 100 trials for Hartmann; results for the other functions are qualitatively similar and are provided in Appendix C.1, together with details on the evaluation. Results for constrained BO using a differentiable relaxation of the feasibility indicator on the sample level are provided in Appendix C.3.

Hyperparameter Optimization: We illustrate the performance of BoTorch on real-world applications, represented by three hyperparameter optimization (HPO) experiments: (1) Tuning 5 parameters of a deep Q-network (DQN) learning algorithm on the Cartpole task from OpenAI gym and the default DQN agent implemented in Horizon , Figure 7; (2) Tuning 6 parameters of a neural network surrogate model for the UCI Adult data set introduced by Falkner et al. , available as part of HPOlib2 , Figure 17 in Appendix C.4; (3) Tuning 3 parameters of the recently proposed Stochastic Weight Averaging (SWA) procedure of Izmailov et al. on the VGG-16 architecture for CIFAR-10, which achieves superior accuracy compared to previously reported results. A more detailed description of these experiments is given in Appendix C.4.

Discussion and Outlook

We presented a novel strategy for effectively optimizing MC acquisition functions using SAA, and established strong theoretical convergence guarantees (in fact, our RQMC convergence results are novel more generally, and of independent interest). Our proposed OKG method, an extension of this approach to “one-shot” optimization of look-ahead acquisition functions, constitutes a significant development of KG, improving scalability and allowing for generic composite objectives and outcome constraints. This approach can naturally be extended to multi-step and other look-ahead approaches .

We make these methodological and theoretical contributions available in our open-source library BoTorch (https://botorch.org), a modern programming framework for BO that features a modular design and flexible API, our distinct SAA approach, and algorithms specifically designed to exploit modern computing paradigms such as parallelization and auto-differentiation. BoTorch is particularly valuable in helping researchers to rapidly assemble novel BO techniques. Specifically, the basic MC acquisition function abstraction provides generic support for batch optimization, asynchronous evaluation, RQMC integration, and composite objectives (including outcome constraints).

Our empirical results show that besides increased flexibility, our advancements in both methodology and computational efficiency translate into significantly faster and more accurate closed-loop optimization performance on a range of standard problems. While other settings such as high-dimensional , multi-fidelity , or multi-objective BO, and non-MC acquisition functions such as Max-Value Entropy Search , are outside the scope of this paper, these approaches can readily be realized in BoTorch and are included in the open-source software package. One can also naturally generalize BO procedures to incorporate neural architectures in BoTorch using standard PyTorch models. In particular, deep kernel architectures , deep Gaussian processes , and variational auto-encoders can easily be incorporated into BoTorch’s primitives, and can be used for more expressive kernels in high-dimensions.

In summary, BoTorch provides the research community with a robust and extensible basis for implementing new ideas and algorithms in a modern computational paradigm, theoretically backed by our novel SAA convergence results.

Bayesian optimization is a generic methodology for optimizing black-box functions, and therefore, by its very nature, not tied to any particular application domain. As mentioned earlier in the paper, Bayesian optimization has been used for various arguably good causes, including drug discovery or reducing the energy footprint of ML applications by reducing the computational cost of tuning hyperparameters. In the Appendix, we give an specific example for how our work can be applied in a public health context, namely to efficiently distribute survey locations for estimating malaria prevalence. BoTorch as a tool specifically has been used in various applications, including transfer learning for neural networks , high-dimensional Bayesian optimization , drug discovery , sim-to-real transfer , trajectory optimization , and nano-material design . However, there is nothing inherent to this work and Bayesian optimization as a field more broadly that would preclude it from being abused in some way, as is the case with any general methodology.

Acknowledgments and Disclosure of Funding

We wish to thank Art Owen for insightful conversations on quasi-Monte-Carlo methods. We also express our appreciation to Peter Frazier and Javier Gonzalez for their helpful feedback on earlier versions of this paper.

Andrew Gordon Wilson is supported by NSF I-DISRE 193471, NIH R01 DA048764-01A1, NSF IIS-1910266, and NSF 1922658 NRT-HDR: FUTURE Foundations, Translation, and Responsibility for Data Science.

References

Appendix A Brief Overview of Other Software Packages for BO

One of the earliest commonly-used packages is Spearmint , which implements a variety of modeling techniques such as MCMC hyperparameter sampling and input warping . Spearmint also supports parallel optimization via fantasies, and constrained optimization with the expected improvement and predictive entropy search acquisition functions . Spearmint was among the first libraries to make BO easily accessible to the end user.

GPyOpt builds on the popular GP regression framework GPy . It supports a similar set of features as Spearmint, along with a local penalization-based approach for parallel optimization . It also provides the ability to customize different components through an alternative, more modular API.

Cornell-MOE implements the Knowledge Gradient (KG) acquisition function, which allows for parallel optimization, and includes recent advances such as large-scale models incorporating gradient evaluations and multi-fidelity optimization . Its core is implemented in C++, which provides performance benefits but renders it hard to modify and extend.

RoBO implements a collection of models and acquisition functions, including Bayesian neural nets and multi-fidelity optimization .

Emukit is a Bayesian optimization and active learning toolkit with a collection of acquisition functions, including for parallel and multi-fidelity optimization. It does not provide specific abstractions for implementing new algorithms, but rather specifies a model API that allows it to be used with the other toolkit components.

The recent Dragonfly library supports parallel optimization, multi-fidelity optimization , and high-dimensional optimization with additive kernels . It takes an ensemble approach and aims to work out-of-the-box across a wide range of problems, a design choice that makes it relatively hard to extend.

Appendix B Parallelism and Hardware Acceleration

Batch evaluation, an important element of modern computing, enables automatic dispatch of independent operations across multiple computational resources (e.g. CPU and GPU cores) for parallelization and memory sharing. All BoTorch components support batch evaluation, which makes it easy to write concise and highly efficient code in a platform-agnostic fashion. Batch evaluation enables fast queries of acquisition functions at a large number of candidate sets in parallel, facilitating novel initialization heuristics and optimization techniques.

Appendix C Additional Empirical Results

This section describes a number of empirical results that were omitted from the main paper due to space constraints.

Algorithms start from the same set of $2d+2$ QMC sampled initial points for each trial, with $d$ the dimension of the design space. We evaluate based on the true noiseless function value at the “suggested point” (i.e., the point to be chosen if BO were to end at this batch). OKG, MOE KG, and NEI use “out-of-sample” suggestions (introduced as $\chi_{n}$ in Section D.4), while the others use “in-sample” suggestions .

All functions are evaluated with noise generated from a $\mathcal{N}(0,.25)$ distribution. Figures 9-11 give the results for all synthetic functions from Section 6. The results show that BoTorch’s NEI and OKG acquisition functions provide highly competitive performance in all cases.

C.2 One-Shot KG Computational Scaling

Figure 13 shows the wall time for generating a set of $q=8$ candidates as a function of the number of total data points $n$ for both standard (Cholesky-based) as well as scalable (Linear CG) posterior inference methods, on both CPU and GPU. While the GPU variants have a significant overhead for small models, they are significantly faster for larger models. Notably, our SAA based OKG is significantly faster than MOE KG, while at the same time achieving much better optimization performance (Figure 13).

C.3 Constrained Bayesian Optimization

We present results for constrained BO on a synthetic function. We consider a multi-output function $f=(f_{1},f_{2})$ and the optimization problem:

Both $f_{1}$ and $f_{2}$ are observed with $\mathcal{N}(0,0.5^{2})$ noise and we model the two components using independent GP models. A constraint-weighted composite objective is used in each of the BoTorch acquisition functions EI, NEI, and OKG.

Results for the case of a Hartmann6 objective and two types of constraints are given in Figures 15-15 (we only show results for BoTorch’s algorithms, since the other packages do not natively support optimization subject to unknown constraints).

The regret values are computed using a feasibility-weighted objective, where “infeasible” is assigned an objective value of zero. For random search and EI, the suggested point is taken to be the best feasible noisily observed point, and for NEI and OKG, we use out-of-sample suggestions by optimizing the feasibility-weighted version of the posterior mean. The results displayed in Figure 15 are for the constrained Hartmann6 benchmark from . Note, however, that the results here are not directly comparable to the figures in because (1) we use feasibility-weighted objectives to compute regret and (2) they follow a different convention for suggested points. We emphasize that our contribution of outcome constraints for the case of KG has not been shown before in the literature.

C.4 Hyperparameter Optimization Details

This section gives further detail on the experimental settings used in each of the hyperparameter optimization problems. As HPO typically involves long and resource intensive training jobs, it is standard to select the configuration with the best observed performance, rather than to evaluate a “suggested” configuration (we cannot perform noiseless function evaluations).

DQN and Cartpole: We consider the case of tuning a deep Q-network (DQN) learning algorithm on the Cartpole task from OpenAI gym and the default DQN agent implemented in Horizon . Figure 17 shows the results of tuning five hyperparameters, exploration parameter (“epsilon”), the target update rate, the discount factor, the learning rate, and the learning rate decay. We allow for a maximum of 60 training episodes or 2000 training steps, whichever occurs first. To reduce noise, each “function evaluation” is taken to be an average of 10 independent training runs of DQN. Figure 17 presents the optimization performance of various acquisition functions from the different packages, using 15 rounds of parallel evaluations of size $q=4$, over 100 trials. While in later iterations all algorithms achieve reasonable performance, BoTorch OKG, EI, NEI, and GPyOpt LP-EI show faster learning early on.

Neural Network Surrogate: We consider the neural network surrogate model for the UCI Adult data set introduced by Falkner et al. , which is available as part of HPOlib2 . We use a surrogate model to achieve a high level of precision in comparing the performance of the algorithms without incurring excessive computational training costs. This is a six-dimensional problem over network parameters (number of layers, units per layer) and training parameters (initial learning rate, batch size, dropout, exponential decay factor for learning rate). Figure 17 shows optimization performance in terms of best observed classification accuracy. Results are means and 95% confidence intervals computed from 200 trials with 75 iterations of size $q=1$. All BoTorch algorithms perform quite similarly here, with OKG doing slightly better in earlier iterations. Notably, they all achieve significantly better accuracy than all other algorithms.

Stochastic Weight Averaging on CIFAR-10: Our final example is for the recently proposed Stochastic Weight Averaging (SWA) procedure of Izmailov et al. , for which good hyperparameter settings are not fully understood. The setting is 300 epochs of training on the VGG-16 architecture for CIFAR-10. We tune three SWA hyperparameters: learning rate, update frequency, and starting iteration using OKG. Izmailov et al. report the mean and standard deviation of the test accuracy over three runs to be $93.64$ and $0.18$, respectively, which corresponds to a 95% confidence interval of $93.64\pm 0.20$. We tune the problem to an average accuracy of $93.84\pm 0.03$.

Appendix D Additional Theoretical Results and Omitted Proofs

Then $\hat{\alpha}_{\!N}^{*}\xrightarrow{a.s.}\alpha^{*}$ and $\textnormal{dist}(\hat{\mathbf{x}}_{\!N}^{*},\mathcal{X}_{f}^{*})\xrightarrow{a.s.}0$.

The following proposition follows directly from Proposition 2.1, Theorem 2.3, and remarks on page 528 of .

D.2 Formal Statement of Theorem 1

$\hat{\alpha}_{\!N}^{*}\rightarrow\alpha^{*}$ a.s., and

$\textnormal{dist}(\hat{\mathbf{x}}_{\!N}^{*},\mathcal{X}^{*})\rightarrow 0$ a.s.

D.3 Randomized Quasi-Monte Carlo Sampling for Sample Average Approximation

In order to use randomized QMC methods with SAA for MC acquisition function, the base samples $E=\{\epsilon^{i}\}$ will need to be generated via RQMC. For the case of Normal base samples, this can be achieved in various ways, e.g. by using inverse CDF methods or a suitable Box-Muller transform of samples $\epsilon^{i}\in^{s}$ (both approaches are implemented in BoTorch). In the language of Section 4, such a transform will become part of the base sample transform $\epsilon\mapsto h(\mathbf{x},\epsilon)$ for any fixed $\mathbf{x}$.

For the purpose of this paper, we consider scrambled $(t,d)$-sequences as discussed by Owen , which are a particular class of RQMC method (BoTorch uses PyTorch’s implementation of scrambled Sobol sequences, which are $(t,d)$-nets in base 2). Using recent theoretical advances from Owen and Rudolf , it is possible to generalize the convergence results from Theorems 1 and 2 to the RQMC setting (to our knowledge, this is the first practical application of these theoretical results).

In the setting of Theorem 1, let $\{\epsilon^{i}\}$ be samples from a $(t,d)$-sequence in base $b$ with gain coefficients no larger than $\Gamma<\infty$, randomized using a nested uniform scramble as in . Then, the conclusions of Theorem 1 still hold. In particular,

$\hat{\alpha}_{\!N_{i}}^{*}\rightarrow\alpha^{*}$ a.s. as $i\rightarrow\infty$,

$\textnormal{dist}(\hat{\mathbf{x}}_{\!N_{i}}^{*},\mathcal{X}^{*})\rightarrow 0$ a.s. as $i\rightarrow\infty$,

In the setting of Theorem 2, let $\{\epsilon^{i}\}$ be samples from $(t,d)$-sequence in base $b$, with gain coefficients no larger than $\Gamma<\infty$, randomized using a nested uniform scramble as in . Then,

Theorem 2(q) as stated does not provide a rate on the convergence of the optimizer. We believe that such result is achievable, but leave it to future work.

Note that while the above results hold for any sequence $(N_{i})_{i}$ with $N_{i}\rightarrow\infty$, in practice the RQMC integration error can be minimized by using sample sizes that exploit intrinsic symmetry of the $(t,d)$-sequences. Specifically, for integers $b\geq 2$ and $M\geq 1$, let

In practice, we chose the MC sample size $N$ from the unique elements of $\mathcal{N}$.

D.4 Asymptotic Optimality of OKG

Suppose conditions (i) and (ii) of Theorem 1 and (iii) of Theorem 2 are satisfied. In addition, suppose that $\limsup_{n}N_{n}=\infty$. Then, $f(\chi_{n})\rightarrow f(x^{*})$ a.s. and in $L^{1}$.

D.5 Proofs

where $h_{\mathcal{D}}(\mathbf{x},\epsilon)=\mu_{\mathcal{D}}(\mathbf{x})+L_{\mathcal{D}}(\mathbf{x})\Phi^{-1}(\epsilon)$ with $\Phi^{-1}$ the inverse CDF of $\mathcal{N}(0,1)$, applied element-wise to the vector $\epsilon$ of qMC samples. Now choose $\epsilon_{0}=(0.5,\dotsc,0.5)$, then

If $f$ is a GP, then $f_{\mathcal{D}_{\mathbf{x}}}(x^{\prime})=h(x^{\prime},\mathbf{x},\epsilon,\epsilon_{I})$, where $\epsilon\sim\mathcal{N}(0,I_{q})$ and $\epsilon_{I}\sim\mathcal{N}(0,1)$ are independent and $h$ is linear in both $\epsilon$ and $\epsilon_{I}$.

This essentially follows from the property of a GP that the covariance conditioned on a new observation $(x,y)$ is independent of $y$.In some cases we may consider constructing a heteroskedastic noise model that results in the function $\sigma^{2}(\mathbf{x})$ changing depending on observations $y$, in which case this argument does not hold true anymore. We will not consider this case further here. We can write $f_{\mathcal{D}_{\mathbf{x}}}(x^{\prime})=\mu_{\mathcal{D}_{\mathbf{x}}}(x^{\prime})+L_{\mathcal{D}_{\mathbf{x}}}^{\sigma}(x^{\prime})\epsilon_{I}$ , where

$L_{\mathcal{D}}^{\sigma}(\mathbf{x})$ is the Cholesky decomposition of $K_{\mathcal{D}}^{\sigma}(\mathbf{x}):=K_{\mathcal{D}}(\mathbf{x},\mathbf{x})+\text{diag}(\sigma^{2}(\mathbf{x}_{1}),\dotsc,\sigma^{2}(\mathbf{x}_{q}))$, and $L_{\mathcal{D}_{\mathbf{x}}}^{\sigma}(x^{\prime})$ is the Cholesky decomposition of

Hence, we see that $f_{\mathcal{D}_{\mathbf{x}}}(x^{\prime})=h(x^{\prime},\mathbf{x},\epsilon,\epsilon_{I})$, with

Bect et al. provide a proof for the case $q=1$. Following their exposition, one finds that the only thing that needs to be verified in order to generalize their results to $q>1$ is that condition (c) in their Definition 3.18 holds also for the case $q>1$. What follows is the multi-point analogue of step (f) in the proof of their Theorem 4.8, which establishes this.

Under some abuse of notation we will use $\mu$ and $K$ also as the vector / matrix-valued mean / kernel function. Let $K^{\sigma}(\mathbf{x}):=K(\mathbf{x},\mathbf{x})+\text{diag}(\sigma(\mathbf{x}))$ and observe that

Moreover, by the analysis above, it holds that

The following Lemma will be used to prove Theorem 4:

For $t<0$, the result is clear because $\||f|\|\geq 0$. ∎

almost surely. Thus, we can use the true KG values as a “potential function” to quantify how the OKG policy performs asymptotically, even though we are never using the KG acquisition function for selecting points. We emphasize that the data that induce $\{\mu_{n}\}_{n\geq 0}$ and $\{\Sigma_{n}\}_{n\geq 0}$ are collected using the OKG policy.

Note that $H_{A}$, for all possible subsets $A$, partition the sample space. Consider some $A\neq\emptyset$. By Lemma A.7 of , if $\alpha_{\textnormal{KG}}(x;\mu_{\infty},\Sigma_{\infty})>0$, then $x$ is measured a finite number of times, meaning that there exists an almost surely finite random variable $M_{0}$ such that on iterations after $N_{0}$, OKG stops sampling from $A$. By the definition of $H_{A}$ in (15), there must exist another random iteration index $M_{1}\geq M_{0}$ such that when $n\geq N_{1}$,

implying that the exact KG policy must prefer points in $\mathcal{A}$ over all others after iteration $M_{1}$. This implies that

Appendix E Illustration of Sample Average Approximation

QMC methods have been used in other applications in machine learning, including variational inference and evolutionary strategies , but rarely in BO. Letham et al. use QMC in the context of a specific acquisition function. BoTorch’s abstractions make it straightforward (and mostly automatic) to use QMC integration with any acquisition function.

Using SAA, i.e., fixing the base samples $E=\{\epsilon^{i}\}$, introduces a consistent bias in the function approximation. While i.i.d. re-sampling in each evaluation ensures that $\hat{\alpha}_{\!N}(\mathbf{x},\Phi,\mathcal{D})$ and $\hat{\alpha}_{\!N}(\mathbf{y},\Phi,\mathcal{D})$ are conditionally independent given $(\mathbf{x},\mathbf{y})$, this no longer holds when fixing the base samples.

Figure 18 illustrates this behavior for EI (we consider the simple case of $q=1$ for which we have an analytic ground truth available). The top row shows the MC and QMC version, respectively, when re-drawing base samples for every evaluation. The solid lines correspond to a single realization, and the shaded region covers four standard deviations around the mean, estimated across 50 evaluations. It is evident that QMC sampling significantly reduces the variance of the estimate. The bottom row shows the same functions for 10 different realizations of fixed base samples. Each of these realizations is differentiable w.r.t. $x$ (and hence $\lambda$ in the slice parameterization). In expectation (over the base samples), this function coincides with the true function (the dashed black line). Conditional on the base sample draw, however, the estimate displays a consistent bias. The variance of this bias (across re-drawing the base samples) is much smaller for the QMC versions.

Figure 20 shows empirical mean and variance of the metrics from Figure 19 as a function of the number of MC samples $N$ on a log-log scale. The stochastic optimizer used is Adam with a learning rate of 0.025. Both for the SAA and the stochastic version we use the same number of random restart initial conditions generated from the same initialization heuristic.

Empirical asymptotic convergence rates can be obtained as the slopes of the OLS fit (dashed lines), and are given in Table 1. It is quite remarkable that in order to achieve the same error as the MC approximation with 4096 samples, the QMC approximation only requires 64 samples. This holds true for the bias and variance of the (relativized) optimal value as well as for the distance from the true optimizer. That said, as we are in a BO setting, we are not necessarily interested in the estimation error $\hat{\alpha}_{\!N}^{*}-\alpha^{*}$ of the optimum, but primarily in how far $x_{\!N}^{*}$ is from the true optimizer $x^{*}$.

A somewhat subtle point is that whether better optimization of the acquisition function results in improved closed-loop BO performance depends on the acquisition function as well as the underlying problem. More exploitative acquisition functions, such as EI, tend to show worse performance for problems with high noise levels. In these settings, not solving the EI maximization exactly adds randomness and thus induces additional exploration, which can improve closed-loop performance. While a general discussion of this point is outside the scope of this paper, BoTorch does provide a framework for optimizing acquisition functions well, so that these questions can be compartmentalized and acquisition function performance can be investigated independently from the quality of optimization.

Perhaps the most significant advantage of using deterministic optimization algorithms is that, unlike for algorithms such as SGD that require tuning the learning rate, the optimization procedure is essentially hyperparameter-free. Figure 9 shows the closed-loop optimization performance of qEI for both deterministic and stochastic optimization for different optimizers and learning rates. While some of the stochastic variants (e.g. ADAM with learning rate 0.01) achieve performance similar to the deterministic optimization, the type of optimizer and learning rate matters. In fact, the rank order of SGD and ADAM w.r.t. to the learning rate is reversed, illustrating that selecting the right hyperparameters for the optimizer is itself a non-trivial problem.

Appendix F Additional Implementation Details

The simplest approach is to use zeroth-order optimizers that do not require gradient information, such as DIRECT or CMA-ES . These approaches are feasible for lower-dimensional problems, but do not scale to higher dimensions. Note that performing parallel optimization over $q$ candidates in a $d$-dimensional feature space means solving a $qd$-dimensional optimization problem.

A more scalable approach incorporates gradient information into the optimization. As described in Section 4, BoTorch by default uses quasi-second order methods, such as L-BFGS-B. Because of the complex structure of the objective, the initial conditions for the algorithm are extremely important so as to avoid getting stuck in a potentially highly sub-optimal local optimum. To reduce this risk, one typically employs multi-start optimization (i.e. start the solver from multiple initial conditions and pick the best of the final solutions). To generate a good set of initial conditions, BoTorch heavily exploits the fast batch evaluation discussed in the previous section. Specifically, BoTorch by default uses $N_{\text{opt}}$ initialization candidates generated using the following heuristic:

F.2 Sequential Greedy Batch Optimization

The pending points approach discussed in Section 5 provides a natural way of generating parallel BO candidates using sequential greedy optimization, where candidates are chosen sequentially, while in each step conditioning on selected points and integrating over the uncertainty in their outcome (using MC integration). By using a full MC formulation, in which we jointly sample at new and pending points, we avoid constructing an individual “fantasy” model for each sampled outcome, a common (and costly) approach in the literature . In practice, the sequential greedy approach often performs well, and may even outperform the joint optimization approach, since it involves a sequence of small, simpler optimization problems, rather than a larger and complex one that is harder to solve.

provide a theoretical justification for why the sequential greedy approach works well with a class of acquisition functions that are submodular.

Appendix G Active Learning Example

Recall from Section 5 the negative integrated posterior variance (NIPV) of the model:

This acquisition function supports both parallel selection of points and asynchronous evaluation. Since MC integration requires evaluating the posterior variance at a large number of points, this acquisition function benefits significantly from the fast predictive variance computations in GPyTorch .

To illustrate how NIPV may be used in combination with scalable probabilistic modeling, we examine the problem of efficient allocation of surveys across a geographic region. Inspired by Cutajar et al. , we utilize publicly-available data from The Malaria Atlas Project (2019) dataset, which includes the yearly mean parasite rate (along with standard errors) of Plasmodium falciparum at a $4.5\text{km}^{2}$ grid spatial resolution across Africa. In particular, we consider the following active learning problem: given a spatio-temporal probabilistic model fit to data from 2011-2016, which geographic locations in and around Nigeria should one sample in 2017 in order to minimize the model’s error for 2017 across all of Nigeria?

We fit a heteroskedastic GP model to 2500 training points prior to 2017 (using a noise model that is itself a GP fit to the provided standard errors). We then select $q=10$ sample locations for 2017 using the NIPV acquisition function, and make predictions across the entirety of Nigeria using this new data. Compared to using no 2017 data, we find that our new dataset reduces MSE by 16.7% on average (SEM = 0.96%) across 60 subsampled datasets. By contrast, sampling the new 2017 points at a regularly spaced grid results only in a 12.4% reduction in MSE (SEM = 0.99%). The mean relative improvement in MSE reduction from NIPV optimization is 21.8% (SEM = 6.64%). Figure 21 shows the NIPV-selected locations on top of the base model’s estimated parasite rate and standard deviation.

Appendix H Additional Implementation Examples

Many of BoTorch’s benefits are qualitative, including the simplification and acceleration of implementing new acquisition functions. Quantifying this in a meaningful way is very challenging. Comparisons are often made in terms of Lines of Code (LoC) - while this metric is problematic when comparing across different design philosophies, non-congruent feature sets, or even programming languages, it does provides a general idea of the effort required for developing and implementing new methods.

MOE’s KG involves thousands of LoC in C++ and python spread across a large number of files,https://github.com/wujian16/Cornell-MOE while our more efficient implementation is <30 LoC. Astudillo and Frazier is a full paper in last year’s installment of this conference,Code available at https://github.com/RaulAstudillo06/BOCF whose composite function method we implement and significantly extend (e.g to support KG) in 7 LoC using BoTorch’s abstractions. The original NEI implementation is >250 LoC, while the one from Code Example LABEL:codex:Abstractions:ImplementationExamples:NEI:qNEI is 14 LoC.

H.1 Composite Objectives

We consider the Bayesian model calibration of a simulator with multiple outputs from Section 5.3 of Astudillo and Frazier . In this case, the simulator from Bliznyuk et al. models the concentrations of chemicals at 12 positions in a one-dimensional channel. Instead of modeling the overall loss function (which measures the deviation of the simulator outputs with a set of observations) directly, we follow Astudillo and Frazier and model the underlying concentrations while utilizing a composite objective approach. A powerful aspect of BoTorch’s modular design is the ability to easily combine different approaches into one. For the composite function problem in this section this means that we can easily extend the work by Astudillo and Frazier not only to use the Knowledge Gradient, but also to the “parallel BO” setting of jointly selecting $q>1$ points. Figures 23 and 23 show results for this with $q=1$ and $q=3$, repspectively. The plots show log regret evaluated at the maximizer of the posterior mean averaged over 250 trials. While the performance of EI-CF is similar for $q=1$ and $q=3$, KG-CF reaches lower regret significantly faster for $q=1$ compared to $q=3$, suggesting that “looking ahead“ is beneficial in this context.

H.2 Generalized UCB

Code Example LABEL:codex:appdx:Example:GenUCB:qUCB presents a generalized version of parallel UCB from Wilson et al. supporting pending candidates, generic objectives, and QMC sampling. If no sampler is specified, a default QMC sampler is used. Similarly, if no objective is specified, the identity objective is assumed.

H.3 Full Code Examples

In this section we provide full implementations for the code examples. Specifically, we include parallel Noisy EI (Code Example LABEL:codex:appdx:Example:NEI:qNEI), OKG (Code Example LABEL:codex:appdx:Example:OSKG:qKG), and (negative) Integrated Posterior Variance (Code Example LABEL:codex:appdx:Example:NIPV:qNIPV).