Budgeted Batch Bayesian Optimization With Unknown Batch Sizes

Parameter settings profoundly impact the performance of machine learning algorithms and laboratory experiments. The classical grid search or trial-error methods are exponentially expensive in large parameter spaces, and Bayesian optimization (BO) offers an elegant alternative for global optimization of black box functions. In situations where the black box function can be evaluated at multiple points simultaneously, batch Bayesian optimization is used. Current batch BO approaches are restrictive in that they fix the number of evaluations per batch, and this can be wasteful when the number of specified evaluations is larger than the number of real maxima in the underlying acquisition function. We present the Budgeted Batch Bayesian Optimization (B3O) for hyper-parameter tuning and experimental design - we identify the appropriate batch size for each iteration in an elegant way. To set the batch size flexible, we use the infinite Gaussian mixture model (IGMM) for automatically identifying the number of peaks in the underlying acquisition functions. We solve the intractability of estimating the IGMM directly from the acquisition function by formulating the batch generalized slice sampling to efficiently draw samples from the acquisition function. We perform extensive experiments for both synthetic functions and two real world applications - machine learning hyper-parameter tuning and experimental design for alloy hardening. We show empirically that the proposed B3O outperforms the existing fixed batch BO approaches in finding the optimum whilst requiring a fewer number of evaluations, thus saving cost and time.


Introduction
Global optimization is fundamental to diverse real-world problems where parameter settings and design choices are pivotal -as an example, in algorithm performance (deep learning networks (Bengio, 2009)) or quality of the products (chemical processes or engineering design (Wang & Shan, 2007)). This requires us to find the global maximum of a non-concave objective function using sequential, and often, noisy observations. Critically, the objective functions are unknown and expensive to evaluate. Therefore, the challenge is to find the maximum of such expensive objective functions in few sequential queries, thus minimizing time and cost.
Bayesian optimization (BO) is an approach to find the global optimum of such expensive, blackbox objective functions f using limited evaluations (Snoek, Larochelle, & Adams, 2012;Shahriari, Swersky, Wang, Adams, & de Freitas, 2016;Wang, Hutter, Zoghi, Matheson, & de Feitas, 2016). Instead of optimizing the real function, BO utilizes a cheaper to evaluate surrogate function, called the acquisition function. This function is used to suggest the next point by balancing exploitation (knowledge of what has been observed) and exploration (where a function has not been investigated). Then, after evaluating the objective function at the suggested point, a Gaussian process (Rasc 2017 AI Access Foundation. All rights reserved. mussen, 2006) is updated and thus gradually a profile of the mean and uncertainty of the exploration space is built up. Bayesian optimization has been demonstrated to outperform other state-of-the-art black-box optimization techniques when function evaluations are expensive and the number of allowed function evaluations is low . Thus, BO has received increasing attention in machine learning community (Thornton, Hutter, Hoos, & Leyton-Brown, 2013;Khajah, Roads, Lindsey, Liu, & Mozer, 2016).
Most selection strategies in BO are sequential, wherein only one experiment is tested at a timethe experiment selection at each iteration is optimized by using the available observed information. However, such methods are inefficient when parallel evaluations are possible. For examples, many alloy samples can be placed in an oven simultaneously (for testing the alloy quality), or several machine learning algorithms can be run in parallel at the same time using multiple cores. This motivates batch Bayesian optimization algorithms that select multiple experiments, a batch of q experiments, at each iteration for evaluation.
Existing batch BO approaches are mostly greedy, sequentially visiting all the maxima of the acquisition function. After the first maximum is found, such methods modify the acquisition function by suppressing the current maximum and then move on to find the next best maximum. This repeats till a batch of maxima is collected. Different algorithms implement this philosophy - (Ginsbourger, Le Riche, & Carraro, 2008) modifies the found maximum point by replacing it with a "fake" or constant value, thus the algorithm is called "constant liar". Other approaches, including GP-BUCB (Desautels, Krause, & Burdick, 2014) and the GP-UCB-PE (Contal, Buffoni, Robicquet, & Vayatis, 2013), exploit the predictive variance of GPs that only depend on the features x, but not the outcome values y. Local Penalization (González, Dai, Hennig, & Lawrence, 2016), on the other hand, using the estimation of Lipschitz constant to penalize the peaks. All these methods update the posterior variance sequentially to modify the acquisition function and thereby derive a batch. Most importantly, these batch BO methods are essentially greedy, choosing individual points until the batch is filled. This is not optimal as after the "real" maxima in the function are found, noisy points will get added simply to complete the batch. An alternate and non-greedy approach is to select a batch of points that reduce the uncertainty of the global maximizer (Shah & Ghahramani, 2015) by extending the Predictive Entropy Search (PES) algorithm (Hernández-Lobato, Hoffman, & Ghahramani, 2014) for a batch setting. But again, a fixed batch size is used. Fixed batch size approaches can restrict the flexibility in finding the global solution and require unnecessary evaluations or miss evaluation of important points. If we over-specify the batch size, we waste resources in evaluating redundant points. In contrast, if we under-specify the batch size, we may miss important points that could potentially be the optimal solution of f .
The best solution is to fill each batch flexibly, but the challenge is to know exactly how many maxima are there in each batch. Our key contribution is to solve this problem, by estimating the number of maxima in the acquisition function, we let the algorithm adjust the batch size at each iteration. This flexible and non-greedy approach has the advantage of being efficient in that it suggests only as many maxima locations as needed. The solution comes from our intuition that the multiple peaks in the acquisition function can be approximated by a mixture of Gaussians -the means of the Gaussian correspond to the underlying peaks. However, since the number of underlying peaks is unknown, we use the infinite Gaussian mixture model (IGMM) (Rasmussen, 1999) so that the number of peaks in the acquisition function can be estimated automatically. Because fitting the IGMM directly to the acquisition function is intractable, we present an efficient batch generalized slice sampler approach to draw samples from acquisition function which are then used to learn the unknown number of peaks. Although the IGMM and the slice sampler have existed for more than a decade, the idea of making use of these techniques for batch Bayesian optimization is novel.
In the experiments, we first validate our methods in 8 benchmark functions and compare it against 9 baselines. We use several criteria for comparison: a) in best-found-value, our method achieves better results in 5 out of 8 cases; and b) our method matches the baseline performance but requires far fewer evaluation points. We further demonstrate the algorithm on machine learning hyper-parameter tuning for support vector regression, multi-label classification and deep learning. We show that our model outperforms the baselines in finding the best values (RMSE, F1 score and accuracy) whilst we require fewer experimental evaluations to reach targets. In addition, we perform a real world experimental design in which we formulate the heat-treatment process required for producing an Aluminum-scandium alloy. We show that our method produces the highest hardness using the fewest number of evaluations. Our result is significant as searching for the best heattreatment process is costly and making this process efficient results in cost and time savings.

Preliminary
We first review the Gaussian process (GP), Bayesian optimization (BO) and the acquisition functions. Then, we summarize the batch Bayesian optimization setting and existing approaches.

Gaussian process
Gaussian processes (GP) (Rasmussen, 2006) extends a multivariate Gaussian distribution to infinite dimensionality. Formally, Gaussian process generates data located throughout some domains such that any finite subset of the range follows a multivariate Gaussian distribution. Given N observations Y = {y 1 , y 2 , ...y N } which can always be imagined as a single point sampled from some multivariate Gaussian distributions.
The mean of GP is assumed to be zero everywhere. What relates one observation to another in such cases is just the covariance function, k (x, x ′ ). From the assumption of GP, we have y ∼ N (0, K) where the covariance matrix is defined as follows: A popular choice for the covariance function k is the squared exponential function: The length parameter l will control this separation when x is not closed to x ′ .
For prediction on a new data point y * , we can update the covariance matrix with k * * = k (x * , x * ) and k * = k (x * , x 1 ) k (x * , x 2 ) · · · k (x * , x N ) . Hence, we can write The conditional probability is followed Gaussian distribution as p (y * | y) ∼ N µ (x * ) , σ 2 (x * ) (Rasmussen, 2006;Ebden, 2008) where its mean and variance are given by GPs provide a full probabilistic model of the data, and allow us to compute not only the model's prediction at input points but also to quantify the uncertainty in the predictions. Therefore, GP is flexible as a nonparametric prior for Bayesian optimization.

Bayesian optimization
We assume that f is a black-box function, that is, its form is unknown and further it is expensive to evaluate. Perturbed evaluations of the type y i = f (x i ) + ε i , where ε i ∼ N 0, σ 2 , are available. Example of the Branin function is in Fig. 1. Bayesian optimization makes a series of evaluations x 1 , ..., x T of f such that the maximum of f is found in the fewest iterations (Snoek et al., 2012;Shahriari et al., 2016;Dai Nguyen, Gupta, Rana, Nguyen, Venkatesh, Deane, & Sanders, 2016;Nguyen, Rana, Gupta, Li, & Venkatesh, 2016c). Formally, let f : X → R be a well behaved function defined on a compact subset X ⊆ R D . Our goal is solving the global optimization problem Bayesian optimization reasons about f by building a Gaussian process through evaluations (Rasmussen, 2006). This flexible distribution allows us to associate a normally distributed random variable at every point in the continuous input space.

ACQUISITION FUNCTIONS
As the original function is expensive to evaluate, the acquisition function acts as a surrogate that determines which point should be selected next. Therefore, instead of maximizing the original function f , we maximize the acquisition function α to select the next point to evaluate In this auxiliary maximization problem, the objective is known and easy to evaluate and can be easily carried out with standard numerical techniques such as multi-start or DIRECT (Jones, Perttunen, & Stuckman, 1993). We consider the case of computing the acquisition function α (x) from the posterior distribution of GP. These acquisition functions are carefully designed to trade off exploration of the search space and exploitation of current promising regions. Among many existing acquisition functions in literature (Hennig & Schuler, 2012;Hernández-Lobato et al., 2014;Srinivas, Krause, Kakade, & Seeger, 2010;Mockus, Tiesis, & Zilinskas, 1978;Jones, 2001;Freitas, Zoghi, & Smola, 2012;Nguyen, Gupta, Rana, Li, & Venkatesh, 2016b), we briefly describe three common acquisition functions including probability of improvement, expected improvement and upper confidence bound. The early work of (Kushner, 1964) suggested maximizing the probability of improvement (PI) over the incumbent α PI ( , where Φ is the standard normal cumulative distribution function (cdf) and the incumbent y + = max x i ∈D t f (x i ). However, the PI exploits quite aggressively. Thus, one could instead measure the expected improvement (EI) (Mockus et al., 1978;Jones, Schonlau, & Welch, 1998). The expected improvement chooses the next point with highest expected improvement over the current best outcome, i.e. the maxi- where φ is the standard normal probability distribution function (pdf), τ is the current best value, u (x) = µ(x)−τ σ (x) , for σ > 0 and zero otherwise. The Gaussian process upper confidence bound (GP-UCB) acquisition function is given by α GP-UCB (x) = µ(x) + β σ (x), where β is a domain-specific positive parameter which trades exploration and exploitation. There are theoretically motivated guidelines for setting and scheduling the hyper-parameter β to achieve sublinear regret (Srinivas et al., 2010).

Batch Bayesian optimization
Bayesian optimization is conventionally posed as a sequential problem where each experiment is completed before taking a new one. In practice it may be advantageous to run multiple function evaluations in parallel. Therefore, we consider evaluating f using a batch of points. An example of batch Bayesian optimization versus sequential Bayesian optimization is illustrated in Fig. 2. As discussed earlier, such scenarios appear, for instance, in the optimization of computer models where several machines (or cores) are available to run experiments in parallel, or in wet-lab experiments when the time of testing one experimental design is the same as testing a batch. Formally, we maximize the acquisition function by finding a collection of points as where n t is the batch size at iteration t. Existing approaches often fix the batch size n t to a constant value for all iterations.

Existing batch Bayesian optimization methods
Since the optimization in Eq. 2 is intractable, batch BO techniques avoid this computational burden by resorting to different strategies. The first BO parallelization was used in the context of the learning Bayesian networks (Očenášek & Schwarz, 2000), to the best of our knowledge. Then batch BO simulation matching (Azimi, Fern, & Fern, 2010) aims to select a batch of size q, which includes points 'close to' the best point. Developing from the expected improvement acquisition function (EI), the multi-point expected improvement (q-EI) (Ginsbourger, Le Riche, & Carraro, 2007) treats the acquisition function as the conditional expectation of the improvement obtained by q points. Constant Liar (CL) is a heuristic q-EI algorithm (Ginsbourger et al., 2008), which uses a greedy approach to iteratively construct a batch of q points. At each iteration, the heuristic uses the sequential algorithm to find a point that maximizes the expected improvement as follows: first, the maximum of the acquisition is found and to move to the next maximum by suppressing this point. This is done by inserting the outcome at this point as a constant value. The GP posterior is updated through the "fake" outcome which then is used to compute the acquisition function. This process is repeated until the batch is filled. Further parallel extension for q-EI are proposed in (Frazier & Clark, 2012;Wang, Clark, Liu, & Frazier, 2015).
Another direction (Contal et al., 2013;Desautels et al., 2014) in batch Bayesian optimization exploits an interesting fact about GPs: the predictive variance of GPs depends only on the feature x, but not the outcome values y. The GP-BUCB algorithm (Desautels et al., 2014) and GP-UCB-PE (Contal et al., 2013) extend the sequential UCB to a batch setting by first selecting the next point, updating the predictive variance which in turn alters the acquisition function, and then selecting the next point. This is repeated till the batch is filled. In particular, the GP-UCB-PE (Contal et al., 2013) chooses the first point of the batch via the UCB score and then defines a "relevance region" and selects the remaining points from this region greedily to maximize the information gain, to focus on pure exploration (PE).
Batch BO can also be developed using information-based policies (Hernández-Lobato et al., 2014). PES aim to select the point which maximizes the information gain of the model. This is solved by expressing the expected reduction in the differential entropy of the predicted distribution. Parallel Predictive Entropy Search (PPES) (Shah & Ghahramani, 2015) extends the Predictive Entropy Search (PES) algorithm of (Hernández-Lobato et al., 2014) to a batch setting.
More recently, Local Penalization (LP) (González et al., 2016) presents a heuristic approach for batch BO by iteratively penalizing the current peak in the acquisition function to find the next peak. LP depends on the estimation of Lipschitz constant to flexibly penalize the peaks. However, in general scenarios, the Lipschitz constant L is unknown. For ease of computation, (González et al., 2016) estimates the Lipschitz constant of the GP predictive mean instead of the actual acquisition function. Furthermore, the use of a unique value of L assumes that the function is Lipschitz homoscedastic. Gonzalez et al (González et al., 2016) has demonstrated that LP outperforms a wide range of baselines in batch Bayesian optimization. We thus consider LP as the most competitive baseline in batch BO.
Although the batch algorithms can speedup the selection of experiments, there are two main drawbacks. First, most of the proposed batch BO (Ginsbourger et al., 2008;Azimi et al., 2010;Azimi, Jalali, & Zhang-fern, 2012;Contal et al., 2013;Desautels et al., 2014;González et al., 2016) involve a greedy algorithm, which chooses individual points until the batch is filled. This is often detrimental as the requirement to fill a batch enforces the selection of not only the true maxima but also noisy peaks. Second, most approaches use a predefined and fixed batch size q for all iterations. Fixed batch size may be inefficient because the number of peaks in the acquisition function is unknown and constantly changing when the GP posterior is updated. For example, if a function has two real maxima, a batch size of more than two will force the selection of noisy points. Therefore, time and resources are wasted for evaluating these noisy points that contravening the goal of BO is to save the number of evaluations. In contrast, if we under-specify the number of peaks, we could miss important points that may be the optimal solution of f . Hence, it is necessary to identify the appropriate, if not exact, number of peaks for evaluation while preserving the ability of finding the optimal value of f .

Budgeted Batch Bayesian Optimization
We propose the novel batch Bayesian optimization method that learns the suitable batch size for each iteration. We term our approach budgeted batch Bayesian optimization (B3O). The proposed method is economic in terms of number of optimal evaluations whilst preserving the performance.
We first describe our approach to approximate the acquisition function using the infinite Gaussian mixture model (IGMM). Then, we present the batch generalized slice sampler to efficiently draw samples under the acquisition function. Next, we briefly describe the existing variational inference technique for IGMM. Finally, we summarize our algorithm.

Batch Bayesian optimization
We consider the batch Bayesian optimization. As a Bayesian optimization task, our ultimate goal is solving the global optimization problem of finding x * = argmax f (x) x∈X by making a series of batch such that the maximum of f is found (Snoek et al., 2012;Shahriari et al., 2016) where n t is the batch size.

Approximating the acquisition function with infinite Gaussian mixture models
The acquisition function α is often multi-modal with the unknown number of peaks. It is intuitive to assume that the acquisition function can be approximated by a mixture of Gaussians (cf. Fig. 3) wherein the peaks in the acquisition function are equivalent to the mean locations of the underlying Gaussians. Because the number of peaks are unknown, we borrow the elegance of Bayesian nonparametrics (Hjort, Holmes, Müller, & Walker, 2010) to identify the unknown number of Gaussian components. We use the infinite Gaussian mixture model (IGMM) (Rasmussen, 1999) which induces the Dirichlet Process prior over possibly infinite number of Gaussian components. In IGMM, each Gaussian component representing a peak k is parameterized by the mean µ k and the covariance Σ k . It is interesting that the estimated mean µ k is the location of the peaks while the estimated covariance Σ k will capture the shape of the peaks in the acquisition function. In contrast, previous work (González et al., 2016) relies on an unique Lipschitz constant to represent the shape of the peak that is not reasonable for the heteroscedastic setting when the peaks have different shapes. Formally, the probability density function of the IGMM is defined as follows ∑ ∞ k=1 π k × N (s | µ k , Σ k ) where s is an observation, π k is the mixing proportion in IGMM and N (s | µ k , Σ k ) is the Gaussian density function of the component k.
Our primary goal in utilizing IGMM is to approximately find the mean µ k (s) of Gaussian distributions as the unknown peaks from the acquisition function. However, directly estimating IGMM from acquisition function is intractable. Therefore, we use the intermediate step to draw samples under the acquisition function and then fit the IGMM to learn the mean locations as the batch of points to evaluate.
Input: #MaxIter, acquisition function α, x min , x max 1: s 0 ∼ uniform (x min , x max ) 2: α min = min x∈X α (x) 3: for i = 1 to #MaxIter do 4: while (notAccept) do 6: if α(s i ) > u i #Accept then We now present the sampling technique to draw samples under the acquisition function for fitting to the IGMM. We note that the sampling process, by its nature, generates more samples from high probability regions. This implies that most of the samples come from the peaks. Thus, even with small number of samples, the sampling process approximates the peaks and is then efficient. There are many existing sampling methods to draw samples from the probability distributions, for example, Metropolis-Hastings (Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953). However, to sample under the D-dimensional acquisition function which is not strictly probability distribution, we utilize the accept-reject sampling (Casella, Robert, & Wells, 2004) where we keep samples in the region under the density function and ignore samples if they are outside the curve.
In particular, we select to use the slice sampling (Neal, 2003) because it is easily implemented for univariate distributions, and can be used to sample from a multivariate distribution by updating each variable in turn. Suppose we wish to sample from a distribution for a variable, x, taking values in some subset of R D , whose density is proportional to some functions g(x). We can do this by sampling uniformly from the (D + 1)-dimensional region that lies under the plot of g(x). This idea can be formalized by introducing an auxiliary real variable, u, and defining a joint distribution over x and u that is uniform over the region R = {(x, u) : 0 < u < g(x)} below the curve or surface defined by g(x). The joint density for (x, u) is Z . To sample for x, we can sample jointly for (x, u), and then ignore u.
We extend the standard slice sampler to the generalized slice sampler that draws samples from D-dimensional acquisition function α which is not a proper distribution and can be negative. To draw u uniformly over the region R below curve of the acquisition function, we get the α min = min x∈X α (x) obtained by one of the non-convex optimization toolbox (e.g., DIRECT (Jones et al., . Accept Figure 4: The illustration of the proposed generalized slice sampling to draw samples under the acquisition functions. We consider an example as follows. The first point s 1 is uniformly drawn from Uniform (x min , x max ), then u 1 ∼ Uniform (α min , α(s 1 )). Next, the second point is sampled s 2 ∼ Uniform (x min , x max ). Due to α(s 2 ) < u 1 , we reject s 2 . We continue to draw s 3 and accept it as α(s 3 ) ≥ u 1 . We repeat the algorithm by again sampling u 3 ∼ Uniform (α min , α(s 3 )).

1993)
). Therefore, we can define the joint density To sample for x, we can sample jointly for (x, u), and then ignore u following the Algorithm 1 and Fig. 4 where uniform (x min , x max ) denotes for drawing randomly from the uniform distribution within the bound of [x min , x max ], which is predefined.
Slice sampling approach (Neal, 2003) also brings a challenge for high-dimensional data that it is hard to find accepted area R such that g (R) > u. There are two main reasons. First, the set of accepted samples could be very small with respect to the hyper-rectangle R. Second, R usually lies in a high-dimensional space, implying the curse of dimensionality. Therefore, it is a big drawback of slice sampling to deal with high dimensions (Tibbits, Haran, & Liechty, 2011;Pietrabissa & Rusconi, 2014).
To overcome the problem in high-dimensional functions, we utilize the generalized slice sampler in batch setting where a bunch of samples are drawn i.i.d at different places. Since our goal is to draw a collection of i.i.d. samples x 1,2,...N from the acquisition function α (to later estimate the IGMM). We define a collection of joint density distribution p (x 1 , u 1 ) , ...p (x M , u) and perform the generalized slice sampling for each joint distribution independently. In the experiment, we set M as 200. In Sec 4.5, we further study the computational efficiency of the batch generalized slice sampling (BGSS) for drawing samples from high dimensional acquisition functions (e.g., D = 50).

Variational inference for infinite Gaussian mixture model
IGMM is the nonparametric mixture model where the prior distribution over the mixing proportion is a Dirichlet process (Ferguson, 1973). There are few existing approaches to learn a IGMM, such as collapsed Gibbs sampler (Kamper, 2013) and variational inference (Blei, Jordan, et al., 2006). In this paper, we follow (Blei et al., 2006) to derive the variational inference for IGMM (Rasmussen, 1999) since the variational approach is generally faster than the Gibbs sampler. After fitting the IGMM using variational inference, we obtain the mean locations µ 1....K as the selected points for evaluations in the batch BO setting. We note that K is unknown and identified automatically in this Bayesian nonparametric setting.
Using the collection of samples (s i ) N i=1 drawn from our batch generalized slice sampler, we consider the IGMM generated as follows. We first draw a mixing proportion from the stick-breaking process given the concentration parameter γ as π where v k ∼ Beta (1, γ). Next, for each component k in the model, we sample the mean µ k ∼ N (µ 0 , I) and the covariance matrix Σ k ∼ Wish (τ k , I). Finally, we generate the assignment z i ∼ Stick (π) and the observation For posterior inference using variational inference, we first write the bound on the likelihood of the data as follows: To define the variational distribution q in the Eq.
(3), we need to construct a distribution on an infinite set of random variables v k ,µ k and Σ k . For this approach be tractable, we truncate the variational distribution at some value K by setting q(v K = 1) = 1 and we can ignore µ k , Σ k for k > K (Blei et al., 2006). Then, we factorize the variational distribution using mean field assumption as the following Particularly, the variational distributions for these variables are defined as q After factorizing the variational distribution using mean field, we maximize the lower bound in Eq. (3). We first take the partial derivative w.r.t. each variational variables. We then equate it to zero and solve the optimization for each variables. Due to the space restriction, we shall refer the detailed inference of the infinite Gaussian mixture model to (Rasmussen, 1999;Blei et al., 2006).

Algorithm
Although the technique of variational inference for IGMM and the slice sampling are existed for years, the idea of connecting these things at the right place for batch BO is novel. We summarize the steps for B3O in Algorithm 2. At an iteration t, we will find a batch including n t points where the number of point n t varies at each iteration.
We note that steps 2,3,7 and 8 are standard in Bayesian optimization techniques. Our proposed method is highlighted in steps 4,5 and 6 to find a batch of points X t from the acquisition function. Particularly, Fig. 3 illustrates and summarizes our steps 4,5 and 6 in 1D and 2D, respectively. We summarize the computation time of steps 4, 5, and 6 in Sec 4.4 and further analyze step 4 by simulation in Sec 4.5 which requires more computation than other steps.

Algorithm 2 Algorithm for Budgeted Batch Bayesian Optimization (B3O).
Input: D 0 = {x i , y i } n 0 i=1 , #iter T , acquisition function α 1: for t = 1 to T do 2: Fit a GP from the data D t .

5:
Fit the Infinite GMM using s.

Experiments
We evaluate the B3O algorithm using both synthetic functions and real-world experimental settings. For the synthetic function evaluation, we utilize 8 functions across dimensions (1-10). In addition, we perform hyper-parameter tuning for three machine learning algorithms -support vector regression , Bayesian nonparametric multi-label classification (Nguyen, Gupta, Rana, Li, & Venkatesh, 2016a) and multi layer perceptron (Ruck, Rogers, Kabrisky, Oxley, & Suter, 1990). We further consider real-world experimental design for Aluminum-Scandium hardening. This heat treatment design involves searching for the best combination of temperatures and times in two stages to achieve the requisite alloy properties -in our case, hardness. We compare our results with baselines on best-found-value and the total number of evaluations, given a fixed number of iterations. We show that we outperform baselines in the best-found-value whilst requiring fewer evaluations. We compare the performance of our proposed approach with the state-of-the-art methods for batch BO.
• Gaussian process -Batch Upper Confidence Bound (GP-BUCB) (Desautels et al., 2014): this is a batch BO utilizing the variance of GP for finding the next point until the batch is filled.
• Rand (EI and UCB): this optimizes the acquisition function to find the first element in the batch (this step is identical to the sequential BO), then chooses a random sample within the bounds until the batch is filled.

Function, Dimension Objective function f (x)
Optimum f (x * ) Table 1: Benchmark functions and dimensions used ranging from 1 to 10.
• Local Penalization (LP) (González et al., 2016): This is currently the state-of-the-art method for batch BO which has been demonstrated to outperform most of other baselines (González et al., 2016). The source code is available in GPyOpt toolbox 1 .

EXPERIMENTAL SETTINGS
We denote the number of iterations by T , and the number of evaluations by N. For the sequential setting of BO, the number of iterations is identical to the number of evaluations, that is N = T while N > T for batch setting. The number of initial points n 0 for all methods is set as 3× D where D is the dimension. For fixed-batch approaches, the batch size at each iteration is set as n t = 3 for functions with D < 5 or n t = D for functions with D ≥ 5. In contrast, our method automatically determines the batch size n t which could differ at each iteration. Thus, we do not need to set the batch size n t for B3O in advance. The performance of the algorithms is compared for a fixed number of iterations T , given as 10 × D, and the total number of evaluated points is N = ∑ T t=0 n t . In all the experiments, the squared exponential (SE) kernel given as k(x, x ′ ) = exp −γ||x − x ′ || 2 is used in the underlying GP, where γ = 0.1 × D and D is the dimension. For methods using the UCB, β is fixed to 2 (following the setting used in (González et al., 2016)), which allows us to compare the different batch designs using the same acquisition function. The results are taken over 20 replicates with different initial values. We always maximize the objective function, maximizing − f for cases in which the goal is to find the minimum. All implementations are in Python. All simulations are done on Windows machine Core i7 Ram 24GB.
Although B3O is designed to work with any kind of acquisition function, in this paper we use B3O with UCB since we empirically notice that the UCB generally works better than EI for B3O. All the source codes and data are available for reproducibility at the link 2 .

TEST FUNCTIONS AND EVALUATION CRITERIA
Test functions are important to validate the performance of optimization algorithms. There have been many benchmark functions reported in the literature (Jamil & Yang, 2013). We select 5 popular 1. https://github.com/SheffieldML/GPyOpt 2. https://github.com/ntienvu/ICDM2016 B3O  benchmark functions, see Table 1. We consider two major criteria for evaluation including the bestfound-value f (x best ) and the total number of evaluations ∑ T t=0 n t .

Comparison -best found value
We examine the performance of B3O in finding the optimum of the chosen benchmark functions. We demonstrate that our method can find better optimal values (minimum) for a fixed amount of iterations T . We report the best-found-value w.r.t. iterations in Fig. 5. Our method achieves significantly better values than the baselines in 5 over 8 functions. In particular, B3O outperforms all baselines in Forrester 1D, Alpine2 5D, gSobol 5D, Alpine2 10D, gSobol 10D. In other three cases, B3O still provides relatively good values. The fixed batch approaches (e.g., LP) may suffer the under-specification (the number of specified batch size n t is smaller than the number of real peaks) at some iterations and thus negatively affect the final performance. As expected, the batch approaches are always better than the sequential approaches using EI and UCB because the number of evaluated points for batch methods is much higher than the sequential. This fact also highlights the advantage of batch BO over the standard (sequential) BO.

Comparison -number of evaluations
The next criteria for comparison is the total number of evaluations, N = ∑ T t=0 n t . The requirement of Bayesian optimization is to keep the number of evaluations to find the optimal value as low as possible, as these evaluations can be costly. For example, it can very expensive to perform a single experiment in material science (to cast an alloy testing), or it takes a few days to train a deep network on a large-scale dataset.
It is not natural to fix the number of evaluations per batch n t , since the number of peaks is unknown, and importantly, these peaks will be changed after new evaluations are done. Therefore, we may waste time and resources if we over-evaluate the number of points than we need. In particular, after all the actual peaks are detected, if we set the batch size n t large, we will get the While the existing approaches fix the batch size n t to a constant for all iterations, our B3O is more efficient by flexibly creating a suitable batch size per iterations. One the one hand, at the iteration 2 in gSobol function, the fixed-batch-size approaches suggest 5 points while we have 6 peaks. Because of underspecifying the peaks, the baselines may miss the optimum. On the other hand, at the iterations 8, 9 and 10, the baselines suggest 5 points while we have only a single peak. This overspecifying the peaks will waste time and resources to evaluate at redundant points.
noisy points which are possibly close to the already detected ones due to the effect of penalizing the peaks (González et al., 2016). Hence, these noisy points are not useful for evaluation. In contrast, under-evaluating the number of necessary points also brings detrimental effects of losing the optimal points.
B3O automatically identifies the suitable number of peaks from the acquisition function in each iteration, and thus does not have to resort to a fixed batch sizes. This is efficient without suffering any performance loss. Given a fixed number of iterations T , we summarize the total number of evaluations in Fig. 6   compared to the fixed-batch baselines, especially for high dimensional settings, e.g., Alpine2 10D and gSobol 10D.

Analysis of B3O -number of points per iteration
We demonstrate that B3O is flexible in identifying the required number of evaluations per batch. We study the number of estimated points in B3O at each iteration. In particular, we record the number of points n t at each iteration, see Fig. 7 for two functions -Dropwave 2D and gSobol 5D. The estimated points per batch in our approach is flexible and determined automatically in this Bayesian nonparametric setting. For some iterations, B3O only recommends one point per batch (see gSobol 5D in Fig. 7). This fact is essential for the homogeneous acquisition function which contains a single peak, whereas the fixed-batch approaches may be wasted to evaluate at some unnecessary points. Moreover, given the fixed batch size for all iterations, other approaches may suffer the effect of over-specification at some iterations and under-specification at other iterations.

Analysis of computational time
Next we compare the computational time at each iteration for the different batch approaches. Specifically, we are interested in the CPU time for finding the batch between different approaches. The sequential methods and the random batch approaches are the cheapest to compute than the other batch approaches. Thus, we do not compare with the sequential and the random batch approaches. Constant Liar (Ginsbourger et al., 2008) and BUCB (Desautels et al., 2014) take time for reestimation of the GP when the fake observations are added and noticeably slower when the number of data points N is large. The Local Penalization approach (González et al., 2016) consumes a considerable amount of time to estimate the Lipschitz constant. In addition, LP computes and maintains the penalized cost around the visited points. The cost for penalizing and optimizing will grow for high dimensions and/or large number of observed points.
B3O is generally competitive to LP (González et al., 2016) in terms of computation. We run faster for low dimensional functions (e.g., less than 6 dimensions). However, LP tends to run faster than B3O for 10D functions. We note that in high dimensional functions, the bottle neck in LP is estimating the Lipschitz constant whilst for B3O the batch generalized slice sampling is the bottle neck due to the nature of the sampling algorithm.

Analysis of batch generalized slice sampling (BGSS)
We investigate the efficiency of the BGSS presented in Section 3.3 for drawing samples under the acquisition function. We vary the observation dimension D -5, 10, 20, 40, 50 -and build the acquisition function α UCB . Then, we design two BGSS settings of size 100 and 200 and compare with the generalized slice sampling which is presented in Fig. 4. We record two important factors including computational time and the number of accepted samples under the curve. In Fig. 8, we present the simulation results that BGSS significantly outperforms the standard slice sampling in terms of computation (running faster) and efficiency (higher number of accepted points). For example, BGSS of the size 200 takes 32 seconds to get 2485 accepted data points at D = 10 dimension. These accepted points are lying under the curve and generally surrounding the peak. Therefore, these samples are beneficial to fit the IGMM since we are particularly interested in finding the peaks.

Tuning hyper-parameters for machine learning algorithms
Hyper-parameter settings greatly impact prediction accuracy of machine learning algorithms. The traditional way of performing hyper-parameter tuning has been to perform a grid search. In practice, however, Bayesian optimization has been shown (Bergstra, Bardenet, Bengio, & Kégl, 2011;Snoek et al., 2012;Thornton et al., 2013) to obtain better results in fewer experiments than required by a full grid search. We employ the batch Bayesian optimization to find the optimal hyper-parameters for support vector regression (SVR), Bayesian nonparametric multi-label classification (BNMC), and multi layer perceptron (MLP).   : Performance comparison on real-world experiments. We consider tuning hyperparameters for three machine learning algorithm: Support vector regression (SVR), Bayesian nonparametric multi-label classification (BNMC), and multi layer perceptron (MLP). We also conduct the Alluminium-scandium heat-treatment design.
Support Vector Regression. The support vector machines (Cortes & Vapnik, 1995) are wellknown for classification problem. It is also applicable for regression task as the support vector regression (SVR) model (Drucker, Burges, Kaufman, Smola, Vapnik, et al., 1997). Using Abalone dataset, we consider tuning the hyper-parameters for SVR with three parameters: C (regularizer parameter), epsilon (ε-insensitive loss) for regression and γ (RBF kernel function). We utilize the Python code in sklearn package.
Bayesian Nonparametric Multi-label Classification. We select to tune the hyper-parameters for the multi-label classification machine learning algorithm of BNMC (Nguyen et al., 2016a) on Scene dataset using the released source code, already constructed as a black-box function 3 . BNMC uses stochastic variational inference (SVI) and stochastic gradient descent (SGD) for learning the hidden correlation of label and feature for multi-label prediction. In particular, we optimize 6 hyperparameters: Dirichlet symmetric for feature and for label, learning rate for SVI and for SGD, truncation threshold and stick-breaking parameter. We aim to maximize the F1 score.
Multi Layer Perceptron. For the deep learning model, we use three layers in the model and evaluate it on the MNIST dataset. Our MLP model has 7 parameters including 4 parameters that include the number of hidden units and dropout rates for the first two layers and 3 parameters of learning rate, decay, and momentum for SGD which is used to learn the MLP. We utilize the Python code in sklearn package.
We compare B3O with baselines in Table 9a (columns 2-3). Our model performs better than all baselines in terms of RMSE for SVR, F1 score for BNMC and has the highest accuracy of 98.47% for MLP. We run 10 iterations for all methods in which the number of initial points is 9 and the batch size n t is 3 for the fixed batch BO. The number of evaluations in B3O is the lowest compared to other methods in this setting. Our model identifies the unknown number of peaks in the acquisition function, then evaluate at these points. Other fixed-batch approaches may evaluate points which are unnecessary. Thus, B3O requires less number of evaluation than the baselines (cf. Fig. 9b). Although we can run these machine learning algorithms in parallel, each evaluation takes a significant amount of time.

Bayesian optimization for experimental design
We consider the alloy hardening process of Aluminum-scandium (Wagner, Kampmann, & Voorhees, 1991) consisting of two stages: nucleation and precipitation coarsening. We aim to maximize the hardness for Aluminum-scandium alloys by designing the appropriate times and temperatures for the two stages in the process. The traditional approach is to use iterative trial and error, in which we choose some material design (temperature and time) based on intuition, past experience, or theoretical knowledge; and test the material in physical experiments; and then use what we learn from these experiments in choosing the material design to try next. This iterative process is repeated until some combination of success or exhaustion is reached! Since the parameter space grows exponentially, trial and error approach is ill-affordable in terms of both time and cost. Therefor, Bayesian optimization is a necessary choice for this experimental design.
We use B3O for the experimental design problem above and show that our method achieves the highest hardness using the smallest number of experiments (see Table 9a). The heat-treatment profile is shown in Fig. 10 using the best estimated parameters (times and temperatures at two steps) 3. https://github.com/ntienvu/ACML2016 BNMC Hour Celsius degree 0 Nucleation Precipitate coarsening Hardness = 88.3 Figure 10: Aluminum-scandium hardening process. We optimize the hardness y w.r.t four features x as temperatures T 1 , T 2 and times t 1 ,t 2 . and the best hardness achieved is 88.3. Given the times and temperatures setting, we evaluate the hardness using the standard kinematic KWN model used by metallurgists (Kampmann & Wagner, 1983) .

Conclusion
We have introduced a novel approach for batch Bayesian optimization. The proposed approach of B3O can identify the suitable batch size at each iteration that most existing approaches in batch BO are unable to do. We have presented the batch generalized slice sampling for drawing samples under the acquisition function. We perform extensive experiments on finding the optimal solution for 8 synthetics functions and evaluate the performance further on 4 real-world tasks. The experimental results highlight the ability of B3O in finding the optimum whilst requiring fewer evaluations than the baselines.

Acknowledgment
We thank Dr Paul G Sanders and Kyle J Deane for the real-world case study of Aluminum-scandium hardening process.