Batched High-dimensional Bayesian Optimization via Structural Kernel Learning

Introduction

Optimization is one of the fundamental pillars of modern machine learning. Considering that most modern machine learning methods involve the solution of some optimization problem, it is not surprising that many recent breakthroughs in this area have been on the back of more effective techniques for optimization. A case in point is deep learning, whose rise has been mirrored by the development of numerous techniques like batch normalization.

While modern algorithms have been shown to be very effective for convex optimization problems defined over continuous domains, the same cannot be stated for non-convex optimization, which has generally been dominated by stochastic techniques. During the last decade, Bayesian optimization has emerged as a popular approach for optimizing black-box functions. However, its applicability is limited to low-dimensional problems because of computational and statistical challenges that arise from optimization in high-dimensional settings.

In the past, these two problems have been addressed by assuming a simpler underlying structure of the black-box function. For instance, Djolonga et al. (2013) assume that the function being optimized has a low-dimensional effective subspace, and learn this subspace via low-rank matrix recovery. Similarly, Kandasamy et al. (2015) assume additive structure of the function where different constituent functions operate on disjoint low-dimensional subspaces. The subspace decomposition can be partially optimized by searching possible decompositions and choosing the one with the highest GP marginal likelihood (treating the decomposition as a hyper-parameter of the GP). Fully optimizing the decomposition is, however, intractable. Li et al. (2016) extended (Kandasamy et al., 2015) to functions with a projected-additive structure, and approximate the projective matrix via projection pursuit with the assumption that the projected subspaces have the same and known dimensions. The aforementioned approaches share the computational challenge of learning the groups of decomposed subspaces without assuming the dimensions of the subspaces are known. Both (Kandasamy et al., 2015) and subsequently (Li et al., 2016) adapt the decomposition by maximizing the GP marginal likelihood every certain number of iterations. However, such maximization is computationally intractable due to the combinatorial nature of the partitions of the feature space, which forces prior work to adopt randomized search heuristics.

In this paper, we develop a new formulation of Bayesian optimization specialized for high dimensions. One of the key contributions of this work is a new formulation that interprets prior work on high-dimensional Bayesian optimization (HDBO) through the lens of structured kernels, and places a prior on the kernel structure. Thereby, our formulation enables simultaneous learning of the decomposition of the function domain.

Prior work on latent decomposition of the feature space considers the setting where exploration/evaluation is performed once at a time. This approach makes Bayesian optimization time-consuming for problems where a large number of function evaluations need to be made, which is the case for high dimensional problems. To overcome this restriction, we extend our approach to a batched version that allows multiple function evaluations to be performed in parallel (Desautels et al., 2014; González et al., 2016; Kathuria et al., 2016). Our second contribution is an approach to select the batch of evaluations for structured kernel learning-based HDBO.

In the past half century, a series of different acquisition functions was developed for sequential BO in relatively low dimensions (Kushner, 1964; Moc̆kus, 1974; Srinivas et al., 2012; Hennig & Schuler, 2012; Hernández-Lobato et al., 2014; Kawaguchi et al., 2015; Wang et al., 2016a; Kawaguchi et al., 2016; Wang & Jegelka, 2017). More recent developments address high dimensional BO by making assumptions on the latent structure of the function to be optimized, such as low-dimensional structure (Wang et al., 2016b; Djolonga et al., 2013) or additive structure of the function (Li et al., 2016; Kandasamy et al., 2015). Duvenaud et al. (2013) explicitly search over kernel structures.

While the aforementioned methods are sequential in nature, the growth of computing power has motivated settings where at once a batch of points is selected for observation (Contal et al., 2013; Desautels et al., 2014; González et al., 2016; Snoek et al., 2012; Wang et al., 2017). For example, the UCB-PE algorithm (Contal et al., 2013) exploits that the posterior variance of a Gaussian Process is independent of the function mean. It greedily selects points with the highest posterior variance, and is able to update the variances without observations in between selections. Similarly, B-UCB (Desautels et al., 2014) greedily chooses points with the highest UCB score computed via the out-dated function mean but up-to-date function variances. However, these methods may be too greedy in their selection, resulting in points that lie far from an optimum. More recently, Kathuria et al. (2016) tries to resolve this issue by sampling the batch via a diversity-promoting distribution for better randomized exploration, while Wang et al. (2017) quantifies the goodness of the batch with a submodular surrogate function that trades off quality and diversity.

Background

Following (Kandasamy et al., 2015), we assume a latent decomposition of the feature dimensions $[D]=\{1,\ldots,D\}$ into disjoint subspaces, namely, $\bigcup_{m=1}^{M}A_{m}=[D]$ and $A_{i}\cap A_{j}=\emptyset$ for all $i\neq j$, $i,j\in[D]$. Further, $f$ can be decomposed into the following additive form:

To make the problem tractable, we assume that each $f_{m}$ is drawn independently from $\mathcal{G}\mathcal{P}(0,k^{(m)})$ for all $m\in[M]$. The resulting $f$ will also be a sample from a GP: $f\sim\mathcal{G}\mathcal{P}(\mu,k)$, where the priors are $\mu(x)=\sum_{m\in[M]}\mu_{m}(x^{A_{m}})$ and $k(x,x^{\prime})=\sum_{m\in[M]}k^{(m)}(x^{A_{m}},{x^{\prime}}^{A_{m}})$. Let $\mathcal{D}_{n}=\{(x_{t},y_{t})\}_{t=1}^{n}$ be the data we observed from $f$, where $y_{t}\sim\mathcal{N}(f(x_{t}),\sigma)$. The log data likelihood for $\mathcal{D}_{n}$ is

where $\bm{K}_{n}=\left[\sum_{m=1}^{M}k^{(m)}(x_{i}^{A_{m}},x_{j}^{A_{m}})\right]_{i\leq n,j\leq n}$ is the gram matrix associated with $\mathcal{D}_{n}$, and $\bm{y}=[y_{t}]_{t\leq n}$ are the concatenated observed function values. Conditioned on the observations $\mathcal{D}_{n}$, we can infer the posterior mean and covariance function of the function component $f^{(m)}$ to be

where $\bm{k}^{(m)}_{n}(x^{A_{m}})=[k^{(m)}(x_{t}^{A_{m}},x^{A_{m}})]_{t\leq n}$.

Learning Additive Kernel Structure

We take a Bayesian view on the task of learning the latent structure of the GP kernel. The decomposition of the input space $\mathcal{X}$ will be learned simultaneously with optimization as more and more data is observed. Our generative model draws mixing proportions $\theta\sim\textsc{Dir}(\alpha)$. Each dimension $j$ is assigned to one out of $M$ groups via the decomposition assignment variable $z_{j}\sim\textsc{Multi}(\theta)$. The objective function is then $f(x)=\sum_{m=1}^{M}f_{m}(x^{A_{m}})$, where $A_{m}=\{{j:z_{j}=m}\}$ is the set of support dimensions for function $f_{m}$, and each $f_{m}$ is drawn from a Gaussian Process. Finally, given an input $x$, we observe $y\sim\mathcal{N}(f(x),\sigma)$. Figure 1 illustrates the corresponding graphical model.

Given the observed data $\mathcal{D}_{n}=\{(x_{t},y_{t})\}_{t=1}^{n}$, we obtain a posterior distribution over possible decompositions $z$ (and mixing proportions $\theta$) that we will include later in the BO process:

Marginalizing over $\theta$ yields the posterior distribution of the decomposition assignment

where $p(\mathcal{D}_{n}|z)$ is the data likelihood (2.1) for the additive GP given a fixed structure defined by $z$. We learn the posterior distribution for $z$ via Gibbs sampling, choose the decomposition among the samples that achieves the highest data likelihood, and then proceed with BO. The Gibbs sampler repeatedly draws coordinate assignments $z_{j}$ according to

and $\bm{K}_{n}^{(z_{j}=m)}$ is the gram matrix associated with the observations $\mathcal{D}_{n}$ by setting $z_{j}=m$. We can use the Gumbel trick to efficiently sample from this categorical distribution. Namely, we sample a vector of i.i.d standard Gumbel variables $\omega_{i}$ of length $M$, and then choose the sampled decomposition assignment $z_{j}=\operatorname{arg\,max}_{i\leq M}\phi_{i}+\omega_{i}$.

With a Dirichlet process, we could make the model nonparametric and the number $M$ of possible groups in the decomposition infinite. Given that we have a fixed number of input dimension $D$, we set $M=D$ in practice.

Diverse Batch Sampling

In real-world applications where function evaluations translate into time-intensive experiments, the typical sequential exploration strategy – observe one function value, update the model, then select the next observation – is undesirable. Batched Bayesian Optimization (BBO) (Azimi et al., 2010; Contal et al., 2013; Kathuria et al., 2016) instead selects a batch of $B$ observations to be made in parallel, then the model is updated with all simultaneously.

Extending this scenario to high dimensions, two questions arise: (1) the acquisition function is expensive to optimize and (2), by itself, does not sufficiently account for exploration. The additive kernel structure improves efficiency for (1). For batch selection (2), we need an efficient strategy that enourages observations that are both informative and non-redundant. Recent work (Contal et al., 2013; Kathuria et al., 2016) selects a point that maximizes the acquisition function, and adds additional batch points via a diversity criterion. In high dimensions, this diverse selection becomes expensive. For example, if each dimension has a finite number of possible valuesWhile we use this discrete categorical domain to illustrate the batch setting, our proposed method is general and is applicable to continuous box-constrained domains., the cost of sampling batch points via a Determinantal Point Process (DPP), as proposed in (Kathuria et al., 2016), grows exponentially with the number of dimensions. The same obstacle arises with the approach by Contal et al. (2013), where points are selected greedily. Thus, naïve adoptions of these approaches in our setting would result in intractable algorithms. Instead, we propose a general approach that explicitly takes advantage of the structured kernel to enable relevant, non-redundant high-dimensional batch selection.

The determinant $\det(\mathbf{K}^{(m)}_{n}(S))$ measures diversity, and hence the DPP assigns higher probability to diverse subsets $S$. An alternative to sampling is to directly maximize the determinant. While this is NP-hard, a greedy strategy gives an approximate solution, and is used in (Kathuria et al., 2016), and in (Contal et al., 2013) as Pure Exploration (PE). We too test this strategy in the experiments. In the beginning, if the GP is not approximating the function well, then greedy may perform no better than a stochastic combination of coordinates, as we observe in Fig. 6.

Sample Combination

Besides this random combination, we can also combine samples greedily. We define a quality function $\psi_{t}^{(m)}$ for each group $m\in[M]$ at time $t$, and combine samples to maximize this quality function. Concretely, for the first point, we combine the maximizers $x_{*}^{(m)}=\arg\max_{x^{(m)}\in\mathcal{X}_{m}}\psi_{t}^{(m)}(x^{(m)})$ from each group. We remove those used parts, $\mathcal{X}_{m}\leftarrow\mathcal{X}_{m}\backslash\{x_{*}^{(m)}\}$, and repeat the procedure until we have $(B-1)$ samples. In each iteration, the sample achieving the highest quality score gets selected, while diversity is retained.

Both selection strategies can be combined with a wide range of existing quality and acquisition functions.

Add-UCB-DPP-BBO

We illustrate the above framework with GP-UCB (Srinivas et al., 2012) as both the acquisition and quality functions. The Upper Confidence Bound $(f_{t}^{(m)})^{+}$ and Lower Confidence Bound $(f_{t}^{(m)})^{-}$ with parameter $\beta_{t}$ for group $m$ at time $t$ are

and combine the expected value $\mu_{t-1}^{(m)}(x)$ of $f_{t}^{(m)}$ with the uncertainty $\beta_{t}^{1/2}\sigma_{t}^{(m)}(x)$. We set both the acquisition function and quality function $\psi_{t}^{(m)}$ to be $(f_{t}^{(m)})^{+}$ for group $m$ at time $t$.

To ensure that we select points with high acquisition function values, we follow (Contal et al., 2013; Kathuria et al., 2016) and define a relevance region $\mathcal{R}_{t}^{(m)}$ for each group $m$ as

where $(y_{t}^{(m)})^{\bullet}=\max_{x^{(m)}\in\mathcal{X}_{m}}(f_{t}^{(m)})^{-}(x^{(m)})$. We then use $\mathcal{R}_{t}^{(m)}$ as the ground set to sample with PE/DPP. The full algorithm is shown in the appendix.

Empirical Results

We empirically evaluate our approach in two parts: First, we verify the effectiveness of using our Gibbs sampling algorithm to learn the additive structure of the unknown function, and then we test our batch BO for high dimensional problems with the Gibbs sampler. Our code is available at https://github.com/zi-w/Structural-Kernel-Learning-for-HDBBO.

We first probe the effectiveness of using the Gibbs sampling method described in Section 3 to learn the decomposition of the input space. More details of the experiments including sensitivity analysis for $\alpha$ can be found in the appendix.

First, we sample test functions from a known additive Gaussian Process prior with zero-mean and isotropic Gaussian kernel with bandwidth $=0.1$ and scale $=5$ for each function component. For $D=2,5,10,20,50,100$ input dimensions, we randomly sample decomposition settings that have at least two groups in the decomposition and at most 3 dimensions in each group.

We set the burn-in period to be 50 iterations, and the total number of iterations for Gibbs sampling to be 100. In Tables 4 and 5, we show two quantities that are closely related to the learned empirical posterior of the decompositions with different numbers of randomly sampled observed data points ($N$). Table 4 shows the probability of two dimensions being correctly grouped together by Gibbs sampling in each iteration of Gibbs sampling after the burn-in period, namely, $(\sum_{i<j\leq D}\mathds{1}_{z^{g}_{i}\equiv z^{g}_{j}\wedge z_{i}\equiv z_{j}})/(\sum_{i<j\leq D}\mathds{1}_{z_{i}\equiv z_{j}})$. Table 5 reports the probability of two dimensions being correctly separated in each iteration of Gibbs sampling after the burn-in period, namely, $(\sum_{i<j\leq D}\mathds{1}_{z^{g}_{i}\neq z^{g}_{j}\wedge z_{i}\neq z_{j}})/(\sum_{i<j\leq D}\mathds{1}_{z_{i}\neq z_{j}})$. The results show that the more data we observe, the more accurate the learned decompositions are. They also suggest that the Gibbs sampling procedure can converge to the ground truth decomposition with enough data for relatively small numbers of dimensions. The higher the dimension, the more data we need to recover the true decomposition.

Effectiveness of Learning Decompositions for Bayesian Optimization

To verify the effectiveness of the learned decomposition for Bayesian optimization, we tested on 2, 10, 20 and 50 dimensional functions sampled from a zero-mean Add-GP with randomly sampled decomposition settings (at least two groups, at most 3 dimensions in each group) and isotropic Gaussian kernel with bandwidth $=0.1$ and scale $=5$. Each experiment was repeated 50 times. An example of a 2-dimensional function component is shown in the appendix. For Add-GP-UCB, we used $\beta^{(m)}_{t}=|A_{m}|\log 2t$ for lower dimensions ($D=2,5,10$), and $\beta^{(m)}_{t}=|A_{m}|\log 2t/5$ for higher dimensions ($D=20,30,50$). We show parts of the results on averaged cumulative regret and simple regret in Fig. 10, and the rest in the appendix. We compare Add-GP-UCB with known additive structure (Known), no partitions (NP), fully partitioned with one dimension for each group (FP) and the following methods of learning the decomposition: Gibbs sampling (Gibbs), randomly sampling the same number of decompositions sampled by Gibbs and select the one with the highest data likelihood (PL-1), randomly sampling 5 decompositions and selecting the one with the highest data likelihood (PL-2). For the latter two learning methods are referred to as “partial learning” in (Kandasamy et al., 2015). The learning of the decomposition is done every 50 iterations. Fig. 3 shows the improvement of learning decompositions with Gibbs over optimizing without partitions (NP).

Overall, the results show that Gibbs outperforms both of the partial learning methods, and for higher dimensions, Gibbs is sometimes even better than Known. Interestingly, similar results can be found in Fig. 3 (c) of (Kandasamy et al., 2015), where different decompositions than the ground truth may give better simple regret. We conjecture that this is because Gibbs is able to explore more than Known, for two reasons:

Empirically, Gibbs changes the decompositions across iterations, especially in the beginning. With fluctuating partitions, even exploitation leads to moving around, because the supposedly “good” points are influenced by the partition. The result is an implicit “exploration” effect that is absent with a fixed partition.

Gibbs sometimes merges “true” parts into larger parts. The parameter $\beta_{t}$ in UCB depends on the size of the part, $|A_{m}|(\log 2t)/5$ (as in (Kandasamy et al., 2015)). Larger parts hence lead to larger $\beta_{t}$ and hence more exploration.

Of course, more exploration is not always better, but Gibbs was able to find a good balance between exploration and exploitation, which leads to better performance. Our preliminary experiments indicate that one solution to ensure that the ground truth decomposition produces the best result is to tune $\beta_{t}$. Hyperparameter selection (such as choosing $\beta_{t}$) for BO is, however, very challenging and an active topic of research (e.g. (Wang et al., 2016a)).

Next, we test the decomposition learning algorithm on a real-world function, which returns the distance between a designated goal location and two objects being pushed by two robot hands, whose trajectory is determined by 14 parameters specifying the location, rotation, velocity, moving direction etc. This function is implemented with a physics engine, the Box2D simulator (Catto, 2011). We use add-GP-UCB with different ways of setting the additive structure to tune the parameters for the robot hand so as to push the object closer to the goal. The regrets are shown in Fig. 4. We observe that the performance of learning the decomposition with Gibbs dominates all existing alternatives including partial learning. Since the function we tested here is composed of the distance to two objects, there could be some underlying additive structure for this function in certain regions of the input space, e.g. when the two robots hands are relatively distant from each other so that one of the hands only impacts one of the objects. Hence, it is possible for Gibbs to learn a good underlying additive structure and perform effective BO with the structures it learned.

2 Diverse Batch Sampling

Next, we probe the effectiveness of batch BO in high dimensions. In particular, we compare variants of the Add-UCB-DPP-BBO approach outlined in Section 4, and a baseline:

Rand: All batch points are chosen uniformly at random from $\mathcal{X}$.

Batch-UCB-*: *$\in\{\texttt{PE},\texttt{DPP}\}$. All acquisition functions are UCB (Eq. 4.1). Exploration is done via PE or DPP with posterior covariance kernels for each group. Combination is via sampling without replacement.

*-Fnc: *$\in\{\texttt{Batch-UCB-PE},\texttt{Batch-UCB-DPP}\}$. All quality functions are also UCB’s, and combination is done by maximizing the quality functions.

A direct application of existing batch selection methods is very inefficient in the high-dimensional settings where they differ more, algorithmically, from our approach that exploits decompositions. Hence, we only compare to uniform sampling as a baseline.

We tested on $2$, $10$, $20$ and $50$-dimensional functions sampled the same way as in Section 5.1; we assume the ground-truth decomposition of the feature space is known. Since Rand performs the worst, we show relative averaged cumulative regret and simple regret of all methods compared to Rand in Fig. 5. Results for absolute values of regrets are shown in the appendix. Each experiment was repeated for $20$ times. For all experiments, we set $\beta^{m}_{t}=|A_{m}|\log 2t$ and $B=10$. All diverse batch sampling methods perform comparably well and far better than Rand, although there exist slight differences. While in lower dimensions ($D\in\{2,10\}$), Batch-UCB-PE-Fnc performs among the best, in higher dimensions ($D\in\{20,50\}$), Batch-UCB-DPP-Fnc performs better than (or comparable to) all other variants. We will see a larger performance gap in later real-world experiments, showing that biasing the combination towards higher quality functions while retaining diversity across the batch of samples provides a better exploration-exploitation trade-off.

For a real-data experiment, we tested the diverse batch sampling algorithms for BBO on the Walker function which returns the walking speed of a three-link planar bipedal walker implemented in Matlab (Westervelt et al., 2007). We tune 25 parameters that may influence the walking speed, including 3 sets of 8 parameters for the ODE solver and 1 parameter specifying the initial velocity of the stance leg. We discretize each dimension into $40$ points, resulting in a function domain of $|\mathcal{X}|=40^{25}$. This size is very inefficient for existing batch sampling techniques. We learn the additive structure via Gibbs sampling and sample batches of size $B=5$. To further improve efficiency, we limit the maximum size of each group to $2$. The regrets for all methods are shown in Fig. 6. Again, all diverse batch sampling methods outperform Rand by a large gap. Moreover, Batch-UCB-DPP-Fnc is a bit better than other variants, suggesting that a selection by quality functions is useful.

Batch Sizes

Finally, we show how the batch size $B$ affects the performance of the proposed methods. We test the algorithms on the $14$-dimensional Robot dataset with $B\in\{5,10\}$. The regrets are shown in Fig. 4. With larger batches, the differences between the batch selection approaches become more pronounced. In both settings, Batch-UCB-DPP-Fnc performs a bit better than other variants, in particular with larger batch sizes.

Conclusion

In this paper, we propose two novel solutions for high dimensional BO: inferring latent structure, and combining it with batch Bayesian Optimization. The experimental results demonstrate that the proposed techniques are effective at optimizing high-dimensional black-box functions. Moreover, their gain over existing methods increases as the dimensionality of the input grows. We believe that these results have the potential to enable the increased use of Bayesian optimization for challenging black-box optimization problems in machine learning that typically involve a large number of parameters.

Acknowledgements

We gratefully acknowledge support from NSF CAREER award 1553284, NSF grants 1420927 and 1523767, from ONR grant N00014-14-1-0486, and from ARO grant W911NF1410433. We thank MIT Supercloud and the Lincoln Laboratory Supercomputing Center for providing computational resources. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of our sponsors.

References

Appendix A Add-UCB-DPP-BBO Algorithm

We present four variants of Add-UCB-DPP-BBO in Algorithm 1. The algorithm framework is general in that, one can plug in other acquisition and quality functions other than UCB to get different algorithms.

Appendix B Additional experiments

In this section, we provide more details in our experiments.

We decompose the acquisition function into $M$ subacquisition functions, one for each part, and optimize those separately. We randomly sample 10000 points in the low dimensional space, and then choose the one with the best value to start gradient descent in the search space (i.e. the range of the box on $R^{|A_{m}|}$). In practice, we observe this approach optimizes low-dimensional ($<5$ dimensions) functions very well. As the number of dimensions grows, the known difficulties of high dimensional BO (and global nonconvex optimization) arise.

B.2 Effectiveness of Decomposition Learning

Sensitivity Analysis for α𝛼\alpha

Empirically, we found that the quality of the learned decompositions is not very sensitive to the scale of $\alpha$ (see Table 6), because the log data likelihood plays a much more important role than $\log(|A_{m}|+\alpha)$ when $\alpha$ is less than the total number of dimensions. The reported results correspond to alpha = 1 for all the partitions.

BO for Synthetic Functions

We show an example of a 2 dimensional function component in the additive synthetic function in Fig. 8. Because of the numerous local maxima, it is very challenging to achieve the global optimum even for 2 dimensions, let alone maximizing an additive sum of them, only by observing their sum. The full results of the simple and cumulative regrets for the synthetic functions comparing Add-GP-UCB with known additive structure (Known), no partitions (NP), fully partitioned with one dimension for each group (FP) and the following methods of learning partition: Gibbs sampling (Gibbs), random sampling the same number of partitions sampled by Gibbs and select the one with the highest data likelihood (PL-1), random sampling 5 partitions and select the one with the highest data likelihood (PL-2) are shown in Fig. 10. The learning was done every 50 iterations, starting from the first iteration. For $D=20,30$, it is quite obvious that when a new partition is learned from the newly observed data (e.g. at iteration 100 and 150), the simple regret gets a boost.

BO for Real-world functions

In addition to be 14 parameter robot pushing task, we tested on the walker function which returns the walking speed of a three-link planar bipedal walker implemented in Matlab (Westervelt et al., 2007). We tune 25 parameters that may influence the walking speed, including 3 sets of 8 parameters for the ODE solver and 1 parameter specifying the initial velocity of the stance leg. To our knowledge, this function does not have an additive structure. The regrets of each decomposition learning methods are shown in Fig. 9. In addition to Gibbs, we test learning decomposition via constrained Gibbs sampling (Gibbs-L), where the maximum size of each group of dimensions does not exceed 2. Because the function does not have additive structure, Gibbs performed poorly since it groups together many dimensions of the input. As a result, its performance is similar to that of no partition (NP). However, Gibbs-L appears to learn a good decomposition with the group size limit, and manages to achieve a slightly lower regret than other methods. Gibbs, PL-1, PL-2 and FP all performed relatively well in for this function, indicating that using the additive structure may benefit the BO procedure even if the function itself is not additive.

B.3 Diverse Batch Sampling

In Fig. 11, we show the full results of the simple and the cumulative regrets on the synthetic functions described in Section 5.2 of the paper.