Batch Bayesian Optimization via Local Penalization

Introduction

Many problems, such as the configuration of machine learning algorithms (Snoek et al., 2012) or the experimental design of biological experiments (González et al., 2014) require the optimization of an unknown, possibly noisy, function $f$. Bayesian optimization (BO) has emerged in this scenario as an efficient heuristic to optimize $f$ if function evaluations are costly and the overall number of evaluations must be kept low (Jones et al., 1998).

The task is to solve the global optimization problem of finding

We assume that $f$ is a black-box from which only perturbed evaluations of the type $y_{i}=f(\mathbf{x}_{i})+\epsilon_{i}$, with $\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$, are available. We will assume that the objective of interest can be described well by a L-Lipschitz continuous function $f:{\mathcal{X}}\to{\rm I\!R}$ defined on a compact subset ${\mathcal{X}}\subseteq{\rm I\!R}^{d}$.

In sequential BO the goal is to make a series of evaluations $\mathbf{x}_{1},\dots,\mathbf{x}_{N}$ of $f$ such that the maximum of $f$ is evaluated as quickly as possible. After $n$ points are available, BO proposes a new location $\mathbf{x}_{n+1}$ using a probabilistic model for $f$, conditioned on all previous observations $\mathcal{D}_{n}=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{n}$. Typically the model is a Gaussian process (GP) $p(f)=\mathcal{GP}(\mu;k)$ with mean function $\mu$ and positive-definite covariance function (kernel) $k$ that in this work we assume is stationary. Under Gaussian likelihoods, the posterior distribution of $f$ (for a sample of size $n$) is also a GP, with posterior mean and variance given by

where $\textbf{K}_{n}$ is the matrix such that $(\textbf{K}_{n})_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$, $\textbf{k}_{n}(\mathbf{x}^{*})=[k(\mathbf{x}_{1},\mathbf{x}^{*}),\dots,k(\mathbf{x}_{n},\mathbf{x}^{*})]^{\top}$ (Rasmussen and Williams, 2005) and $\mathbf{x}^{*}$ is the point where the GP is evaluated.

This posterior is used to form the acquisition function $\alpha(\mathbf{x};\mathcal{I}_{n})$, where $\mathcal{I}_{n}$ represents the available data set $\mathcal{D}_{n}$ and the GP structure (kernel, likelihood and parameter values) when $n$ data points are available. The next evaluation is placed at the (numerically estimated) global maximum $\mathbf{x}_{n+1}$ of this acquisition function. A number of possible acquisition functions are now available, ranging from fast heuristics (Osborne, 2010; Jones et al., 1998) to non-local entropy-based approaches (Hennig and Schuler, 2012; Hernández-Lobato et al., 2014).

While the goal of Bayesian optimization is to keep the number of evaluations of $f$ as low as possible, in high-dimensional and or otherwise complex problems, the number of required evaluations can still be considerable. Parallel approaches arise as the natural solution to circumvent the computational bottleneck around these evaluations of $f$. We focus on cases in which the cost of evaluating $f$ in a batch of points of size $n_{b}$ is the same as evaluating $f$ in a single point. Such scenarios appear, for instance, in the optimization of computer models where several cores are available to run in parallel, or in wet-lab experiments when the cost of testing one experimental design is the same as testing a batch of them. In these settings, the set of available pairs $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{n}$ can be augmented with the evaluations of $f$ on batches of data points $\mathcal{B}_{t}^{n_{b}}=\{\mathbf{x}_{t,1},\dots,\mathbf{x}_{t,nb}\}$, for $t=1,\dots,m$, rather than on single observations. Our goal here is to define a design rule for such batches $\mathcal{B}_{1}^{n_{b}},\dots,\mathcal{B}_{m}^{n_{b}}$. The batch selection problem can be generalized further, e.g. by adapting the batch size (Azimi et al., 2012) or by collecting batches asynchronously (Ginsbourger et al., 2011; Janusevskis et al., 2012; Snoek et al., 2012). For simplicity of exposition these ideas will not feature further here.

The goal of any batch criterion is to mimic the decisions that would be made under the equivalent (optimal) sequential policy: Consider the choice of selecting $\mathbf{x}_{t,k}$, the $k$-th element of the $t$-th batch. Under a sequential policy, in which the evaluations of $f$ at all locations prior to $\mathbf{x}_{t,k}$ are available, the decision is to take $\mathbf{x}_{t,k}$ as the maximizer of $\alpha(\mathbf{x};\mathcal{I}_{t,k-1})$. In the batch case, the decision about where to collect $\mathbf{x}_{t,k}$ has to incorporate the uncertainty about the locations $\mathbf{x}_{t,1},\dots,\mathbf{x}_{t,k-1}$, and the outcomes of the evaluation of $f$ there. Iteratively marginalizing these sources of uncertainty gives

is the predictive distribution of the GP at $\mathbf{x}_{t,j}$ when a total of $n$ points are available and

reflects the optimization step required to obtain $\mathbf{x}_{t,j}$ after the evaluations of $f$ at previous batch-elements have been marginalized.

The optimization in Eq. (1.1) is intractable even for small batch-sizes, due to the optimization-marginalization loop required to obtain $\mathbf{x}_{t,k}$. The literature in batch BO has tried to avoid this computational burden by means of different strategies, most of which involve the explicit use of the predictive distributions $p(y_{t,j}|\mathbf{x}_{t,j},\mathcal{I}_{t,j-1})$, for $j=1,\dots,n_{b}$. Exploratory approaches (Schonlau et al., 1998; Contal et al., 2013) search for a reduction in system uncertainty. This is using the property that the variance of $p(y_{t,j}|\mathbf{x}_{t,j},\mathcal{I}_{t,j-1})$ does not depend on the value of the objective there. Other methods use $p(y_{t,j}|\mathbf{x}_{t,j},\mathcal{I}_{t,j-1})$ to generate ‘fake’ observations of the model (Azimi et al., 2012, 2011; Bergstra et al., 2011) and avoid the marginalization step. In statistics, the suitability of the expected improvement utility has been studied for the design of batches (Chevalier and Ginsbourger, 2013; Frazier, 2012). In contrast to the previous mentioned works, these methods use the joint distribution of $y_{t_{1}},\dots y_{t,nb}$ to simultaneously optimize elements on the batch (Azimi et al., 2010). These non-greedy strategies are very well founded from a theoretical perspective in practice but tend to scale poorly with the dimension of the problem and the sizes of the batches. Other theoretical properties of batch BO have been studied in the context of Bayesian networks (Očenášek and Schwarz, 2000), multi-armed bandits (Desautels et al., 2012), and the optimal balance between exploration and exploitation in batch designs (Jalali et al., 2013).

2 Goal and Contributions of this work

Using $p(y_{t,j}|\mathbf{x}_{t,j},\mathcal{I}_{t,j-1})$ to model the interaction between batch elements has a computational overhead of $\mathcal{O}(n^{3})$, since the GP needs to be updated after each batch location is selected to jointly optimize all the elements in the batch. The motivation of this work is to develop a heuristic approximation of Eq. (1.1) at lower computational cost, while incorporating information about global properties of $f$ from the GP model into the batch design.

Our approach rests on the hypothesis that $f$ is a Lipschitz continuous function, which is a common assumption in global optimization (Floudas and Pardalos, 2009). For easy reference: a real-valued function $f:{\mathcal{X}}\to{\rm I\!R}$ on a compact subset ${\mathcal{X}}\subseteq{\rm I\!R}^{d}$ of the $d$-dimensional real vector space is said to be $L$-Lipschitz if it satisfies

In the context of parallelizing Bayesian optimization, a beneficial aspect of the Lipschitzian assumption is that it naturally allows us to place bounds on how far the optimum of $f$ is from a certain location. See Figure 1 for details. As explained below, this information can be used to define policies to collect a batch of points multiple steps ahead without evaluating $f$, by mimicking the hypothesized behavior of a sequential policy. The main challenge is that, in practice, the constant $L$ is unknown. In the literature, this problem has been addressed from different angles (Floudas and Pardalos, 2009). We explore a new alternative: inferring the Lipschitz constant directly from the Gaussian process model for $f$.

Our contributions are: (i) A new batch BO heuristic, BBO-LP, that selects batches of points by an iterative maximization-penalization loop around the the acquisition function. This leads to efficient parallelization of BO and can be used with any acquisition function. (ii) A probabilistic framework to approximately infer the Lipschitz constant of $f$, termed GP-LCA, that uses the properties of the gradients of the GP. The inferred value of $L$ is used to improve batch selection. (iii) A python implementation of several batch BO methods is published in conjunction with this work.http://sheffieldml.github.io/GPyOpt/ (iv) Confirmation of the effectiveness of the approach is demonstrated through several simulated experiments, an algorithm configuration problem, and a real wet-lab experimental design. In particular, the local penalization approach performs equal or better than current batch BO methods in terms of the convergence to the maximum, but shows better performance in terms of gained information per second.

Maximization-Penalization Strategy for Batch Design

The intuition behind our approach is that for most GP priors in practical use for BO, the dominant effect of a function evaluation on the acquisition function is a local exclusion around the new evaluation. This shape of the acquisition function will be modeled through the Lipschitz properties of $f$, to distribute the elements in each batch. This should be understood as a heuristic to the shape of $\alpha(\mathbf{x};\mathcal{I}_{t,k-1})$ if all previous observations were available, mimicking the effect a sequential policy. This is especially useful in cases in which the acquisition function shows multi-modal shape, a common situation in the first iterations of BO algorithms. The following definition is helpful for the formalization of the algorithm:

A function $\varphi(\mathbf{x};\mathbf{x}_{j})$, $\mathbf{x}\in\mathcal{X}$, is a local penalizer of a generic acquisition function $\alpha(\mathbf{x})$ at $\mathbf{x}_{j}$ if $\varphi(\mathbf{x};\mathbf{x}_{j})$ is differentiable, $0\leq\varphi(\mathbf{x};\mathbf{x}_{j})\leq 1$ and $\varphi(\mathbf{x};\mathbf{x}_{j})$ is an non-decreasing function in $\|\mathbf{x}-\mathbf{x}_{j}\|$.

We propose to replace the maximization-marginalization loop in Eq. 1.1 by a maximization-penalization strategy: while the optimization is carried out in a similar fashion, the marginalization step is replaced by the direct penalization of $\alpha(\mathbf{x};\mathcal{I}_{t,k-1})$ around its most recent maximum, i.e, the previous batch element. Figure 2 gives a graphical illustration. The maximization-penalization strategy selects $\mathbf{x}_{t,k}$ as

where $\varphi(\mathbf{x};\mathbf{x}_{t,j})$ are local local penalizers centered at $\mathbf{x}_{t,j}$ and $g:{\rm I\!R}\rightarrow{\rm I\!R}^{+}$ is a differentiable transformation of $\alpha(\mathbf{x})$ that keeps it strictly positive without changing the location of its extrema. We will use $g(z)=z$ if $\alpha(\mathbf{x})$ is already positive and the soft-plus transformation $g(z)=\ln(1+e^{z})$ elsewhere. This does not require re-estimation of the GP model after each location is selected, just a new optimization of the penalized utility.

The effect of a local penalizer is to smoothly reduce the value of the acquisition function in a neighborhood of $\mathbf{x}_{j}$. A ‘good’ local penalizer centered at $\mathbf{x}_{j}$ should reflect the belief about the distance from $\mathbf{x}_{j}$ to $\mathbf{x}_{M}$: If we suspect that $\mathbf{x}_{M}$ is far from $\mathbf{x}_{j}$, a broad $\varphi(\mathbf{x};\mathbf{x}_{j})$ will discard a large portion of $\mathcal{X}$ in which we don’t need to collect any sample. On the other hand, if we believe that $\mathbf{x}_{M}$ and $\mathbf{x}_{j}$ are close, ideally we want to minimize the penalization of $\alpha(\mathbf{x})$ and keep collecting samples is a close neighborhood. This local penalization mimics the acquisition function’s dynamics under a sequential policy in the following sense: the modes of the acquisition functions correspond to regions in which either $\mu_{n}(\mathbf{x})$ or $\sigma^{2}_{n}(\mathbf{x})$ (or both) are large. Evaluating, for instance, where $\sigma_{n}(\mathbf{x})$ is large will reduce uncertainty in that region, decreasing $\alpha(\mathbf{x})$ in a neighborhood. The functions $\varphi(\mathbf{x};\mathbf{x}_{j})$ are surrogates for this neighborhood.

We now construct penalizing functions $\varphi(\mathbf{x};\mathbf{x}_{j})$ that incorporate into $\alpha(\mathbf{x})$ the current belief about the distance from the batch locations to $\mathbf{x}_{M}$. Take $M=\max_{\mathbf{x}\in{\mathcal{X}}}f(\mathbf{x})$, and a valid Lipschitz constant $L$. Consider the ball

To simplify the notation we write $r_{j}=r(\mathbf{x}_{j})$ for the radius of the ball around $\mathbf{x}_{j}$. If $f$ in (5) is the true optimization objective, then $\mathbf{x}_{M}\notin B_{r_{j}}(\mathbf{x})$—otherwise the Lipschitz condition would be violated. The size of $B_{r_{j}}(\mathbf{x}_{j})$ depends on $L$, $M$ and the value of $f$ at $\mathbf{x}_{j}$. Both large variability in $f$ (large $L$) and proximity of $f(\mathbf{x}_{j})$ to the optimum $M$ shrink $B_{r_{j}}(\mathbf{x}_{j})$.

In the BO context, under the assumption $f(\mathbf{x})\sim\mathcal{GP}(\mu(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime}))$, we choose $\varphi(\mathbf{x};\mathbf{x}_{j})$ as the probability that $\mathbf{x}$, any point in $\mathcal{X}$ that is a potential candidate to be a maximum, does not belong to $B_{r_{j}}(\mathbf{x}_{j})$:

The following proposition (proof in Supp. Materials A) shows that this local penalizer can be computed in closed form.

Let $f(\mathbf{x})$ be a $\mathcal{GP}$ with posterior mean $\mu_{n}(\mathbf{x})$ and posterior variance $\sigma^{2}_{n}(\mathbf{x})$. The function $\varphi(\mathbf{x};\mathbf{x}_{j})$ in Eq. (6) is a valid local penalizer of $\alpha(\mathbf{x})$ at $\mathbf{x}_{j}$ such that:

where $z=\frac{1}{\sqrt{2\sigma_{n}^{2}(\mathbf{x}_{j})}}\left(L\|\mathbf{x}_{j}-\mathbf{x}\|-M+\mu_{n}(\mathbf{x}_{j})\right),$ for erfc the complementary error function, $M=\max_{\mathbf{x}\in{\mathcal{X}}}f(\mathbf{x})$ and $L$ a valid Lipschitz constant.

The functions $\varphi(\mathbf{x};\mathbf{x}_{j})$ thus create exclusion zones whose size is governed by $L$. If $\mu_{n}(\mathbf{x}_{j})$ is close to $M$, then $\varphi(\mathbf{x};\mathbf{x}_{j})$ will have a smaller and more localized effect on $\alpha(\mathbf{x})$ (a smaller exclusion area). On the other hand, if $\mu_{n}(\mathbf{x}_{j})$ is far from $M$, $\varphi(\mathbf{x};\mathbf{x}_{j})$ will produce a wider yet less intense correction on $\alpha(\mathbf{x})$. The value of $L$ also affects the size of the effect of $\varphi(\mathbf{x};\mathbf{x}_{j})$ on $\alpha(\mathbf{x})$, decreasing it as $L$ increases.

2 Selecting the parameters L𝐿L and M𝑀M

The values of $M$ and $L$ are unknown in general. To approximate $M$, one can take $\hat{M}=\max_{\mathcal{X}}\mu_{n}(\mathbf{x})$ or, to avoid solving this maximization problem, use the even rougher approximation $\hat{M}=\max_{i}\{y_{i}\}$.

Regarding the parameter $L$ note that the definition of Lipschitz continuity in Eq. (3) does not uniquely identify $L$. In the BO penalization context, small but feasible values of $L$ are preferred, because they produce large exclusion zones and thus more efficient search. Given access to the true objective $f$, one can show that $L_{\nabla}=\max_{\mathbf{x}\in\mathcal{X}}\|{\nabla f(\mathbf{x})}\|$ is a valid Lipschitz constant (see Supp. Material C for further details). Note that $L_{\nabla}$ is the smallest value of $L$ that satisfies the Lipschitz condition Eq. 3 in the limit $\mathbf{x}_{1}\to\mathbf{x}_{2}$ in (3).

We now construct an approximation for $L_{\nabla}$. Assuming that $f$ is a draw from a GP with a (at least) twice differentiable kernel $k$, the gradient of $f$ at $\mathbf{x}^{*}$ is distributed as a multivariate Gaussian $\nabla f(\mathbf{x}^{*})|\textbf{X},\textbf{y},\mathbf{x}^{*}\sim\mathcal{N}(\mu_{\nabla}(\mathbf{x}^{*}),\Sigma_{\nabla}^{2}(\mathbf{x}^{*}))$ with mean vector

and call this the Gaussian Process Lipschitz Constant Approximation criterion (GP-LCA). Note that this definition of $\hat{L}_{GP-LCA}$ ignores the variance of the gradient, which could be used to identify candidate points to improve the approximation of $L_{\nabla}$ in a Bayesian optimization fashion. The supplement contains further experiments supporting the quality of this approximation. See Algorithm 1 for a description of all the steps described in this section.

3 Heteroscedastic scenarios

The use of an unique value of $L$ assumes that the function to optimize is Lipschitz homocedastic. Although this is a typical hypothesis for most BO methods, recent works have pointed out that some real problems do not satisfy this condition (Assael et al., 2014). It is not the goal of this work to analyze this case further but, interestingly, the method proposed here can be extended to non-Lipschitz cases by replacing $L$ in the penalizers $\varphi(\mathbf{x};\mathbf{x}_{j},\hat{L})$ by a local values of $L$. For instance, a possible approach would be to replace the local penalizers by $\varphi(\mathbf{x};\mathbf{x}_{j},\hat{L}_{j})$ where $\hat{L}_{j}=\|\mu_{\nabla}(\mathbf{x}_{j})\|$.

4 Optimizing the penalized acquisition function

where $\nabla\alpha(\textbf{x};\mathcal{I}_{t,0})$ is the (assumed known) gradient of the original acquisition function and $\nabla\varphi(\textbf{x};\textbf{x}_{t,j})$ are the gradients of the local penalizers

Experimental Section

This section compares the performance of Algorithm 1 with the state-of-the-art methods for batch BO. We label the different methods by means of the batch design type followed by the acquisition used: Rand is used when the first element in the batch is collected maximizing the acquisition and the remaining ones randomly, B and PE denote the exploratory approaches in (Schonlau et al., 1998) and (Contal et al., 2013), Pred is used in cases when the model is used to generate ‘fake’ batch observations as in (Azimi et al., 2012), SM identifies the simulating and matching method (Azimi et al., 2010) and LP stands for our local penalization method. The multi-point expected improvement (Chevalier and Ginsbourger, 2013) is denoted by qEI. Two acquisition functions are used: the expected improvement (EI) defined as $\alpha_{EI}(\mathbf{x};\mathcal{I}_{n})=\left(u\Phi(u)+\phi(u)\right)\sigma_{n}(\mathbf{x}),$ where $u=(\mu_{n}(\mathbf{x})-y_{\min})/\sigma_{n}(\mathbf{x})$ and $\Phi(\cdot)$, $\phi(\cdot)$ are the standard Gaussian distribution and density functions respectively and $y_{\min}$ is the best current location and the Upper Confidence Bound (UCB) defined as $\alpha_{UCB}(\mathbf{x};\mathcal{I}_{n})=\mu_{n}(\mathbf{x})+\kappa\sigma_{n}(\mathbf{x})$, with $\kappa\geq 0$. The batch methods that can be used with an arbitrary acquisition function are tested using both, with the exception of the SM whose implementation is only available with the UCB. When used in a sequential setting (for baseline reference) the EI and UCB are referred by their acronyms. In total, we use 2 sequential and 10 batch methods. To run the B, PE, SM, methods, we use the available Matlab code.http://econtal.perso.math.cnrs.fr/software/. Note that an alternative implementation of the GP-B-UCB code is available at http://www.its.caltech.edu/ tadesaut/GPBUCBCode/ but we used the former one for consistency in the comparisons. The implementation of these methods optimize $f$ by searching its optimum in a fine grid, which is an advantage computationally but a drawback in terms of precision. The qEI was taken from the R-package DiceOptim.http://cran.r-project.org/web/packages/DiceOptim Unless specified otherwise, the default implemented settings of all the previous methods are used.

We perform: (i) a simulation in which the performance of the algorithms is compared for a fixed time budget across different problem dimensions, batch sizes and acquisition functions and (ii) a comparison of the gained information per second rate in three objective functions with different evaluation costs. We always minimize the objective, minimizing $-f$ in examples in which the goal is to find maximum of $f$. In all the experiments the exponentiated quadratic (EQ) covariance $k(\mathbf{x},\mathbf{x}^{\prime})=\theta\exp(-\gamma\|\mathbf{x}-\mathbf{x}^{\prime}\|^{2})$, $\theta,\gamma>0$ is used in the GP, whose parameters are optimized by maximizing the marginal likelihood from the best of 10 random initializations. The results are taken over 20 replicates with different initial values. All the simulations were done on Amazon EC2 servers with Intel Xeon E5-2666 processors and 2 virtual CPUs except the SVR tuning with 16 virtual CPUs.

We consider the gSobol function (see Supp. Materials D) to compare the above mentioned methods across dimensions $d=2,5,10$ and batch sizes, $n_{b}=5,10,20$. For methods using the UCB, $\kappa$ was fixed to 2, which allows us to compare the different batch designs using the same acquisition function. For dimension 2, 5 and 10, we use a time budget of 1, 5 and 10 mins. respectively. Table 1 shows the averaged best value found by each algorithm for all the iterations completed within the time limit. In general, the batch methods using the UCB show a better performance that methods using the EI, especially in dimensions 5 and 10. The overall best technique is the LP-UCB, that achieves the best results in 5 of the 9 cases. It is also notable that it exhibits fairly small standard deviations compared with the rest of the methods and it is coherent accumulating information about the optimum of $f$ in terms of the batch size: as $n_{b}$ increases the results are consistently better. In dimension 2 and batch size 5, the LP-EI is the best method. There are three cases in which the LP batch designs are not the most competitive (although still providing good results). Exploratory approaches works well in low dimensional cases, being the B-UCB the best method in two scenarios.

2 Comparisons in terms of the cost to evaluate the objective

We choose three scenarios to compare the algorithms in terms of the running time. The examples correspond to three functions that are cheap, moderate and expensive to evaluate. More specifically, the first experiment uses a function (Cosines) that is inexpensive to evaluate but quite multi-modal. The second experiment is motivated by a wet-lab experimental design. We work with a surface that emulates the performance of mammalian cells in protein production given different gene designs. The function has dimension 71 and is is moderately expensive to evaluate since it corresponds to the predictive mean of a GP trained over 1,500 data instances. The qEI was not used in this experiment due to the huge computational effort required to jointly optimize the batches in dimension 71. The third experiment involves the tuning of the three parameters of a support vector regression (SVR) (Drucker et al., 1997) in a example with 45730 instances and 9 continuous attributes (Bache and Lichman, 2013). The objective function is the cross-validation error of the model, which is expensive to evaluate due to the amount of data used. See Supp. Materials D for further details. We take a batch size of $n_{b}=5$ for the Cosines function and $n_{b}=10$ for the wet-lab and SVR experiments. We compare the averaged best found results in terms of the number of collected batches and the wall-clock time. In the last experiment we use the SVR implementation available in scikit–learnhttp://scikit-learn.org/stable/index.html. and only the methods implemented in python are used (EI, UCB, Rand-EI, Rand-UCB, Pred-EI, Pred-UCB, LP-EI and LP-UCB).

In the Cosines experiment both the sequential EI and UCB policies achieve the best results during the first 10 iterations of the algorithms (2 full batches). As the algorithms progress, however, a significant improvement is observed by the LP-EI and LP-UCB methods in terms of the number iterations and in terms of the wall-clock-time. When many points are collected, the update of the GP is more expensive and a good batch design is able to explore regions that the sequential method cannot. The rest of the batch methods, however, are not able to do this exploration efficiently, which leads to poorer results. Similar results are obtained for the wet-lab experiment. The LP-EI and LP-UCB are again the most competitive techniques improving the rest of the batch methods and the sequential policies. The differences are even more significant in this scenario. Since $f$ is now more expensive to evaluate, the parallelization of the evaluations makes the search much more efficient, specially for the LP-UCB method. Regarding the last experiment, the cost of evaluating the function dominates the cost of designing the batch. In this case the performance of the different batch methods is comparable but significantly better than the sequential policies due to the parallel evaluations of $f$. The results for the three functions are coherent with those observed in Section 3.1 showing that the BBP-LP methods is overall the most efficient method for batches collection in BO.

Discussion

We have investigated a new heuristic for batch BO, BBO-LP, that significantly reduces the computational burden of non-parallelizable tasks. The resulting method can be used with any acquisition function and it is able to make fast and appropriate decisions about the locations where $f$ should be evaluated. When the batch evaluations of $f$ are parallelizable this is an important advantage, meaning that they don’t lead to considerable additional computational overhead. We have found other interesting results. In simple scenarios, batch policies based on random exploration work reasonably well in terms of the information gained per second. When the complexity of the problem increases, however, methods that make use of some information about $f$ improve the random policy. In particular, the approach here proposed makes use of the Lipschitz continuity of $f$ to model the interaction between the elements in the batch. In spirit, this is similar to use the GP to predict the evaluations of $f$ but, in practice, is much more efficient because it avoids the re-computation of the GP after every point is selected. The limitations of this approach are, however, determined by the ability to learn correctly a small enough, and valid, Lipschitz constant for $f$.

One could also wonder whether it is necessary to require that sample paths from the GP measure on $f$ should be Lipschitz-continuous themselves. This would severely restrict the applicability of this notion, because the relationship between regularity of the kernel and the sample paths is complicated. Even if the kernel is Lipschitz-continuous, sample paths may not be Lipschitz (Adler, 1981). However, our approach only tries to model the effect of evaluations on the BO objective, not the GP probability measure itself. Many BO objectives, in particular the EI and UCB, are smooth functions of only the sufficient statistics (mean and covariance function) of the GP posterior. Both the posterior mean and covariance function are members of the Reproducing kernel Hilbert space induced by the kernel (i.e. they are weighted sums of kernel functions). Thus, if the kernel is Lipschitz, so is the acquisition function, even if the GP measure itself has non-Lipschitz sample paths. Finally, note that our local repulsion criterion naturally suggests a Latin square design for the case when no functional values have been been acquired. The latin square design is widely suggested for this domain (Jones et al., 1998).

References

Appendix A Proof of Proposition 1

We compute the explicit form of the penalization functions $\varphi(\textbf{x};\textbf{x}_{j})$. The distribution of $r_{j}$ is Gaussian with mean $(M-\mu_{n}(\textbf{x}_{j}))/L$ and variance $\sigma_{n}^{2}(\textbf{x}_{j})/L^{2}$ by the properties of $f(\textbf{x}_{j})$. Then we obtain that

Appendix B Optimization of the penalized acquisition function

Under the proposed local penalization method, to select the $k$-th element of the $t$-th batch requires the optimization of the function

which can be done by any gradient descend method as follows. We fist map the problem into the natural log space by observing that

Applying the properties of the logarithms we transform the problem into the maximization of

The gradient with respect to x is now easy to calculate since the problem is in additive form. First, note that the gradients of the local penalizers $\nabla\varphi(\textbf{x};\textbf{x}_{t,j})$ are

When $\alpha(\textbf{x};\mathcal{I}_{t,0})$ is not necessarily positive one can take $g(z)=\exp(z)$ and the gradient simplifies to

Appendix C Lipschitz constant approximation

In this section we elaborate in the approximation of the Lipschitz constant. First, we include the following proposition that allows to uniquely identify a valid value of $L$.

Let $f:{\mathcal{X}}\to{\rm I\!R}$ be a L-Lipschitz continuous function defined on a compact subset ${\mathcal{X}}\subseteq{\rm I\!R}^{d}$. Take

where $\nabla f(\textbf{x})=\left(\frac{\partial f}{\partial\textbf{x}_{1}},\cdots,\frac{\partial f}{\partial\textbf{x}_{p}}\right)^{\top}$. Then, $L_{p}$ is a valid Lipschitz constant such that the Lipschitz condition

where $\frac{1}{s}+\frac{1}{l}=1$, holds.

Using the mean value theorem for every $\textbf{x}_{1},\textbf{x}_{2}\in\mathcal{X}$ there exist a $\textbf{w}=\textbf{x}_{1}+\beta\textbf{x}_{2}$, with $\beta\in(0,1)$ such that,

Since $\textbf{w}\in\mathcal{X}$ by definition, we have that

for $L_{p}=\max_{\textbf{x}\in\mathcal{X}}\|{\nabla f(\textbf{x})}\|_{p}.$

In order to test the empirical approximation of the Lipschitz constant detailed in Section 2.2 we use the Cosines function described in the experimental section of this work. The true $L_{\nabla}$ for this function is $8.808636$, that was calculated by maximizing the norm of gradient of $f$ in a very fine grid. We check the quality of our approximation to $L_{\nabla}$ for increasing sample size up to 50 observations, where the locations of the points are randomly selected along the domain of $f$ using a bivariate uniform distribution. The evaluations of $f$ at the selected locations were perturbed with Gaussian noise with standard deviations $\sigma=0,0.1,0.25$. In Figure 4 we show the results for 30 replicates of the experiment. The average approximation of L converges to the true $L_{\nabla}$, being this convergence slower when the evaluation errors increase.

Appendix D Detailed description of the experiments

Table 2 contains the the details of the functions used in the experiments of this work.

D.2 Gene design experiment

There is an increasing interest in the pharmacological industry in the the design of synthetic genes capable of transforming cells into ‘factories’ able to produce drugs of interest. In this experiment we emulate a gene design process.

The function to maximize is the production of cell proteins, that it is known depends on certain features of the gene sequences. We built a GP to link gene features and protein production efficiency based the model described in González et al. . A total of 71 gene features are considered, which correspond to the dimension of the final design space. We validated the model with the remaining 2908 genes of the dataset and we used its posterior mean as the function to optimize. We can understand this model as a mathematical surrogate of the cell behavior in which the mean evaluations play the role of physical wet-lab gene design tests, many of which can be run in parallel for the same price of one.

D.3 SVR parameter tuning experiment

Support Vector Machines (SVR) for regression Drucker et al. with an EQ kernel, depend on three parameters: the kernel lengthscale ($\gamma$), the soft margin parameter ($C$) and the band size ($\epsilon$). A proper choice of the parameters is crucial to guarantee a good performance of the SVR, which is typically done by minimizing the mean square error (RMSE) in a test dataset. This task can be expensive, specially for large datasets. We use BO to optimize the parameters of the SVR using the ‘Physiochemical’ properties of protein tertiary structure’ dataset available in the UCI Machine Learning repository Bache and Lichman . This dataset has 45,730 instances and 9 continuous attributes that are used to predict the coordinate root mean square distance (RMSD), a measure that describes the distance per residue between to optimally aligned protein sequences. We trained the SVR using a randomly selected subset of 22,000 proteins and we tested the results of using the rest. Every iteration takes around 300 seconds.