Max-value Entropy Search for Efficient Bayesian Optimization

Introduction

Bayesian optimization (BO) has become a popular and effective way for black-box optimization of nonconvex, expensive functions in robotics, machine learning, computer vision, and many other areas of science and engineering (Brochu et al., 2009; Calandra et al., 2014; Krause & Ong, 2011; Lizotte et al., 2007; Snoek et al., 2012; Thornton et al., 2013; Wang et al., 2017). In BO, a prior is posed on the (unknown) objective function, and the uncertainty given by the associated posterior is the basis for an acquisition function that guides the selection of the next point to query the function. The selection of queries and hence the acquisition function is critical for the success of the method.

Different BO techniques differ in this acquisition function. Among the most popular ones range the Gaussian process upper confidence bound (GP-UCB) (Auer, 2002; Srinivas et al., 2010), probability of improvement (PI) (Kushner, 1964), and expected improvement (EI) (Moc̆kus, 1974). Particularly successful recent additions are entropy search (ES) (Hennig & Schuler, 2012) and predictive entropy search (PES) (Hernández-Lobato et al., 2014), which aim to maximize the mutual information between the queried points and the location of the global optimum.

ES and PES are effective in the sense that they are query-efficient and identify a good point within competitively few iterations, but determining the next query point involves very expensive computations. As a result, these methods are most useful if the black-box function requires a lot of effort to evaluate, and are relatively slow otherwise. Moreover, they rely on estimating the entropy of the $\operatorname{arg\,max}$ of the function. In high dimensions, this estimation demands a large number of samples from the input space, which can quickly become inefficient.

We propose a twist to the viewpoint of ES and PES that retains the information-theoretic motivation and empirically successful query-efficiency of those methods, but at a much reduced computational cost. The key insight is to replace the uncertainty about the $\operatorname{arg\,max}$ with the uncertainty about the maximum function value. As a result, we refer to our new method as Max-value Entropy Search (MES). As opposed to the $\operatorname{arg\,max}$, the maximum function value lives in a one-dimensional space, which greatly facilitates the estimation of the mutual information via sampling. We explore two strategies to make the entropy estimation efficient: an approximation by a Gumbel distribution, and a Monte Carlo approach that uses random features.

Our contributions are as follows: (1) MES, a variant of the entropy search methods, which enjoys efficient computation and simple implementation; (2) an intuitive analysis which establishes the first connection between ES/PES and the previously proposed criteria GP-UCB, PI and EST (Wang et al., 2016), where the bridge is formed by MES; (3) a regret bound for a variant of MES, which, to our knowledge, is the first regret bound established for any variant of the entropy search methods; (4) an extension of MES to the high dimensional settings via additive Gaussian processes; and (5) empirical evaluations which demonstrate that MES identifies good points as quickly or better than ES/PES, but is much more efficient and robust in estimating the mutual information, and therefore much faster than its input-space counterparts.

After acceptance of this work, we learned that Hoffman & Ghahramani (2015) independently arrived at the acquisition function in Eq. (5). Yet, our approximation (Eq. (6)) is different, and hence the actual acquisition function we evaluate and analyze is different.

Background

Gaussian processes (GPs) are distributions over functions, and popular priors for Bayesian nonparametric regression. In a GP, any finite set of function values has a multivariate Gaussian distribution. A Gaussian process $GP(\mu,k)$ is fully specified by a mean function $\mu({\boldsymbol{x}})$ and covariance (kernel) function $k({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})$. Let $f$ be a function sampled from $GP(\mu,k)$. Given the observations $D_{t}=\{({\boldsymbol{x}}_{\tau},y_{\tau})\}_{\tau=1}^{t}$, we obtain the posterior mean $\mu_{t}({\boldsymbol{x}})=\boldsymbol{k}_{t}({\boldsymbol{x}})^{\textrm{T}}(\boldsymbol{K}_{t}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{t}$ and posterior covariance $k_{t}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})=k({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})-\boldsymbol{k}_{t}({\boldsymbol{x}})^{\textrm{T}}(\boldsymbol{K}_{t}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{k}_{t}({\boldsymbol{x}}^{\prime})$ of the function via the kernel matrix $\boldsymbol{K}_{t}=\left[k({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j})\right]_{{\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j}\in D_{t}}$ and $\boldsymbol{k}_{t}({\boldsymbol{x}})=[k({\boldsymbol{x}}_{i},{\boldsymbol{x}})]_{{\boldsymbol{x}}_{i}\in D_{t}}$ (Rasmussen & Williams, 2006). The posterior variance is $\sigma^{2}_{t}({\boldsymbol{x}})=k_{t}({\boldsymbol{x}},{\boldsymbol{x}})$.

2 Additive Gaussian Processes

Additive Gaussian processes (add-GP) were proposed in (Duvenaud et al., 2011), and analyzed in the BO setting in (Kandasamy et al., 2015). Following the latter, we assume that the function $f$ is a sum of independent functions sampled from Gaussian processes that are active on disjoint sets $A_{m}$ of input dimensions. Precisely, $f(x)=\sum_{m=1}^{M}f^{(m)}(x^{A_{m}})$, with $A_{i}\cap A_{j}=\emptyset$ for all $i\neq j$, $|\cup_{i=1}^{M}A_{i}|=d$, and $f^{(m)}\sim GP(\mu^{(m)},k^{(m)})$, for all $m\leq M$ ($M\leq d<\infty$). As a result of this decomposition, the function $f$ is distributed according to $GP(\sum_{m=1}^{M}\mu^{(m)},\sum_{m=1}^{M}k^{(m)})$. Given a set of noisy observations $D_{t}=\{({\boldsymbol{x}}_{\tau},y_{\tau})\}_{\tau=1}^{t}$ where $y_{\tau}\sim\mathcal{N}(f(x_{\tau}),\sigma^{2})$, the posterior mean and covariance of the function component $f^{(m)}$ can be inferred as $\mu_{t}^{(m)}({\boldsymbol{x}})=\boldsymbol{k}^{(m)}_{t}({\boldsymbol{x}})^{\textrm{T}}(\boldsymbol{K}_{t}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{t}$ and $k_{t}^{(m)}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})=k^{(m)}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})-\boldsymbol{k}^{(m)}_{t}({\boldsymbol{x}})^{\textrm{T}}(\boldsymbol{K}_{t}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{k}^{(m)}_{t}({\boldsymbol{x}}^{\prime})$, where $\boldsymbol{k}^{(m)}_{t}({\boldsymbol{x}})=[k^{(m)}({\boldsymbol{x}}_{i},{\boldsymbol{x}})]_{{\boldsymbol{x}}_{i}\in D_{t}}$ and $\boldsymbol{K}_{t}=\left[\sum_{m=1}^{M}k^{(m)}({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j})\right]_{{\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j}\in D_{t}}$. For simplicity, we use the shorthand $k^{(m)}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})=k^{(m)}({\boldsymbol{x}}^{A_{m}},{\boldsymbol{x}}^{\prime A_{m}})$.

3 Evaluation Criteria

Max-value Entropy Search

Entropy search methods use an information-theoretic perspective to select where to evaluate. They find a query point that maximizes the information about the location ${\boldsymbol{x}}_{*}=\operatorname{arg\,max}_{{\boldsymbol{x}}\in\mathfrak{X}}f(x)$ whose value $y_{*}=f({\boldsymbol{x}}_{*})$ achieves the global maximum of the function $f$. Using the negative differential entropy of $p({\boldsymbol{x}}_{*}|D_{t})$ to characterize the uncertainty about ${\boldsymbol{x}}_{*}$, ES and PES use the acquisition functions

ES uses formulation (2), in which the expectation is over $p(y|D_{t},{\boldsymbol{x}})$, while PES uses the equivalent, symmetric formulation (3), where the expectation is over $p({\boldsymbol{x}}_{*}|D_{t})$. Unfortunately, both $p({\boldsymbol{x}}_{*}|D_{t})$ and its entropy is analytically intractable and have to be approximated via expensive computations. Moreover, the optimum may not be unique, adding further complexity to this distribution.

We follow the same information-theoretic idea but propose a much cheaper and more robust objective to compute. Instead of measuring the information about the argmax ${\boldsymbol{x}}_{*}$, we use the information about the maximum value $y_{*}=f({\boldsymbol{x}}_{*})$. Our acquisition function is the gain in mutual information between the maximum $y_{*}$ and the next point we query, which can be approximated analytically by evaluating the entropy of the predictive distribution:

where $\psi$ is the probability density function and $\Psi$ the cumulative density function of a normal distribution, and $\gamma_{y_{*}}({\boldsymbol{x}})=\frac{y_{*}-\mu_{t}({\boldsymbol{x}})}{\sigma_{t}({\boldsymbol{x}})}$. The expectation in Eq. (5) is over $p(y_{*}|D_{n})$, which is approximated using Monte Carlo estimation by sampling a set of $K$ function maxima. Notice that the probability in the first term $p(y|D_{t},{\boldsymbol{x}})$ is a Gaussian distribution with mean $\mu_{t}({\boldsymbol{x}})$ and variance $k_{t}({\boldsymbol{x}},{\boldsymbol{x}})$. The probability in the second term $p(y|D_{n},{\boldsymbol{x}},y_{*})$ is a truncated Gaussian distribution: given $y_{*}$, the distribution of $y$ needs to satisfy $y<y_{*}$. Importantly, while ES and PES rely on the expensive, $d$-dimensional distribution $p({\boldsymbol{x}}_{*}|D_{t})$, here, we use the one-dimensional $p(y_{*}|D_{n})$, which is computationally much easier.

It may not be immediately intuitive that the value should bear sufficient information for a good search strategy. Yet, the empirical results in Section 5 will demonstrate that this strategy is typically at least as good as ES/PES. From a formal perspective, Wang et al. (2016) showed how an estimate of the maximum value implies a good search strategy (EST). Indeed, Lemma 3.1 will make the relation between EST and a simpler, degenerate version of MES explicit.

Hence, it remains to determine how to sample $y_{*}$. We propose two strategies: (1) sampling from an approximation via a Gumbel distribution; and (2) sampling functions from the posterior Gaussian distribution and maximizing the functions to obtain samples of $y_{*}$. We present the MES algorithm in Alg. 1.

The marginal distribution of $f(x)$ for any $x$ is a one-dimensional Gaussian, and hence the distribution of $y^{*}$ may be viewed as the maximum of an infinite collection of dependent Gaussian random variables. Since this distribution is difficult to compute, we make two simplifications. First, we replace the continuous set $\mathfrak{X}$ by a discrete (finite), dense subset $\hat{\mathfrak{X}}$ of representative points. If we select $\hat{\mathfrak{X}}$ to be an $\epsilon$-cover of $\mathfrak{X}$ and the function $f$ is Lipschitz continuous with constant $L$, then we obtain a valid upper bound on $f(\mathfrak{X})$ by adding $\epsilon L$ to any upper bound on $f(\hat{\mathfrak{X}})$.

Second, we use a “mean field” approximation and treat the function values at the points in $\hat{\mathfrak{X}}$ as independent. This approximation tends to over-estimate the maximum; this follows from Slepian’s lemma if $k(x,x^{\prime})\geq 0$. Such upper bounds still lead to optimization strategies with vanishing regret, whereas lower bounds may not (Wang et al., 2016).

We sample from the approximation $\hat{p}(y^{*}|D_{n})$ via its cumulative distribution function (CDF) $\widehat{\Pr}[y_{*}<z]=\prod_{{\boldsymbol{x}}\in\hat{\mathfrak{X}}}\Psi(\gamma_{z}({\boldsymbol{x}}))$. That means we sample $r$ uniformly from $$and find$z$such that$\Pr[y_{*}$. A binary search for$z$to accuracy$\delta$requires$O(\log\tfrac{1}{\delta})$queries to the CDF, and each query takes$O(|\hat{\mathfrak{X}}|)\approx O(n^{d})$time, so we obtain an overall time of$O(M|\hat{\mathfrak{X}}|\log\frac{1}{\delta})$for drawing$M$ samples.

To sample more efficiently, we propose a $O(M+|\hat{\mathfrak{X}}|\log\frac{1}{\delta})$-time strategy, by approximating the CDF by a Gumbel distribution: $\widehat{\Pr}[y_{*}<z]\approx\mathcal{G}(a,b)=e^{-e^{-\frac{z-a}{b}}}$. This choice is motivated by the Fisher-Tippett-Gnedenko theorem (Fisher, 1930), which states that the maximum of a set of i.i.d. Gaussian variables is asymptotically described by a Gumbel distribution (see the appendix for further details). This does not in general extend to non-i.i.d. Gaussian variables, but we nevertheless observe that in practice, this approach yields a good and fast approximation.

We sample from the Gumbel distribution via the Gumbel quantile function: we sample $r$ uniformly from $$, and let the sample be$y=\mathcal{G}^{-1}(a,b)=a-b\log(-\log r)$. We set the appropriate Gumbel distribution parameters$a$and$b$by percentile matching and solve the two-variable linear equations$a-b\log(-\log r_{1})=y_{1}$and$a-b\log(-\log r_{2})=y_{2}$, where$\Pr[{y}_{*}$and$\Pr[{y}_{*}$. In practice, we use$r_{1}=0.25$and$r_{2}=0.75$so that the scale of the approximated Gumbel distribution is proportional to the interquartile range of the CDF$\hat{\Pr}[y_{*}

3 Relation to Other BO Methods

As a side effect, our new acquisition function draws connections between ES/PES and other popular BO methods. The connection between MES and ES/PES follows from the information-theoretic viewpoint; the following lemma makes the connections to other methods explicit.

MES, where we only use a single sample $y_{*}$ for $\alpha_{t}(x)$;

GP-UCB with $\beta^{\frac{1}{2}}=\min_{{\boldsymbol{x}}\in\mathfrak{X}}\frac{y_{*}-\mu_{t}({\boldsymbol{x}})}{\sigma_{t}({\boldsymbol{x}})}$;

This equivalence no longer holds if we use $M>1$ samples of $y_{*}$ in MES.

The equivalence among 2,3,4 is stated in Lemma 2.1 in (Wang et al., 2016). What remains to be shown is the equivalence between 1 and 2. When using a single $y_{*}$ in MES, the next point to evaluate is chosen by maximizing $\alpha_{t}({\boldsymbol{x}})=\gamma_{y_{*}}({\boldsymbol{x}})\frac{\psi(\gamma_{y_{*}}({\boldsymbol{x}}))}{2\Psi(\gamma_{y_{*}}({\boldsymbol{x}}))}-\log(\Psi(\gamma_{y_{*}}({\boldsymbol{x}})))$ and $\gamma_{y_{*}}=\frac{y_{*}-\mu_{t}({\boldsymbol{x}})}{\sigma_{t}({\boldsymbol{x}})}$. For EST with $m=y_{*}$, the next point to evaluate is chosen by minimizing $\gamma_{y_{*}}({\boldsymbol{x}})$. Let us define a function $g(u)=u\frac{\psi(u)}{2\Psi(u)}-\log(\Psi(u))$. Clearly, $\alpha_{t}({\boldsymbol{x}})=g(\gamma_{y_{*}}({\boldsymbol{x}}))$. Because $g(u)$ is a monotonically decreasing function, maximizing $g(\gamma_{y_{*}}({\boldsymbol{x}}))$ is equivalent to minimizing $\gamma_{y_{*}}({\boldsymbol{x}})$. Hence 1 and 2 are equivalent. ∎

4 Regret Bound

The connection with EST directly leads to a bound on the simple regret of MES, when using only one sample of $y_{*}$. We prove Theorem 3.2 in the appendix.

Let $F$ be the cumulative probability distribution for the maximum of any function $f$ sampled from $GP(\mu,k)$ over the compact search space $\mathfrak{X}\subset R^{d}$, where $k({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})\leq 1,\forall{\boldsymbol{x}},{\boldsymbol{x}}^{\prime}\in\mathfrak{X}$. Let $f_{*}=\max_{{\boldsymbol{x}}\in\mathfrak{X}}f({\boldsymbol{x}})$ and $w=F(f_{*})\in(0,1)$, and assume the observation noise is iid $\mathcal{N}(0,\sigma)$. If in each iteration $t$, the query point is chosen as ${\boldsymbol{x}}_{t}=\operatorname{arg\,max}_{{\boldsymbol{x}}\in\mathfrak{X}}\gamma_{y^{t}_{*}}({\boldsymbol{x}})\frac{\psi(\gamma_{y^{t}_{*}}({\boldsymbol{x}}))}{2\Psi(\gamma_{y^{t}_{*}}({\boldsymbol{x}}))}-\log(\Psi(\gamma_{y^{t}_{*}}({\boldsymbol{x}})))$, where $\gamma_{y^{t}_{*}}({\boldsymbol{x}})=\frac{y^{t}_{*}-\mu_{t}({\boldsymbol{x}})}{\sigma_{t}({\boldsymbol{x}})}$ and $y^{t}_{*}$ is drawn from $F$, then with probability at least $1-\delta$, in $T^{\prime}=\sum_{i=1}^{T}\log_{w}\frac{\delta}{2\pi_{i}}$ number of iterations, the simple regret satisfies

where $C=2/\log(1+\sigma^{-2})$ and $\zeta_{T}=(2\log(\frac{\pi_{T}}{\delta}))^{\frac{1}{2}}$; $\pi$ satisfies $\sum_{i=1}^{T}\pi_{i}^{-1}\leq 1$ and $\pi_{t}>0$, and $t^{*}=\operatorname{arg\,max}_{t}\nu_{t}$ with $\nu_{t}\triangleq\min_{{\boldsymbol{x}}\in\mathfrak{X},y^{t}_{*}>f_{*}}\gamma_{y^{t}_{*}}({\boldsymbol{x}})$, and $\rho_{T}$ is the maximum information gain of at most $T$ selected points.

5 Model Adaptation

In practice we do not know the hyper-parameters of the GP, so we must adapt our GP model as we observe more data. A standard way to learn the GP hyper-parameters is to optimize the marginal data likelihood with respect to the hyper-parameters. As a full Bayesian treatment, we can also draw samples of the hyper-parameters using slice sampling (Vanhatalo et al., 2013), and then marginalize out the hyper-parameters in our acquisition function in Eq. (6). Namely, if we use $E$ to denote the set of sampled settings for the GP hyper-parameters, our acquisition function becomes

where $\gamma^{\eta}_{y_{*}}({\boldsymbol{x}})=\frac{y_{*}-\mu^{\eta}_{t}({\boldsymbol{x}})}{\sigma^{\eta}_{t}({\boldsymbol{x}})}$ and the posterior inference on the mean function $\mu_{t}^{\eta}$ and $\sigma_{t}^{\eta}$ depends on the GP hyper-parameter setting $\eta$. Similar approaches have been used in (Hernández-Lobato et al., 2014; Snoek et al., 2012).

High Dimensional MES with Add-GP

The high-dimensional input setting has been a challenge for many BO methods. We extend MES to this setting via additive Gaussian processes (Add-GP). In the past, Add-GP has been used and analyzed for GP-UCB (Kandasamy et al., 2015), which assumed the high dimensional black-box function is a summation of several disjoint lower dimensional functions. Utilizing this special additive structure, we overcome the statistical problem of having insufficient data to recover a complex function, and the difficulty of optimizing acquisition functions in high dimensions.

Since the function components $f^{(m)}$ are independent, we can maximize the mutual information between the input in the active dimensions $A_{m}$ and maximum of $f^{(m)}$ for each component separately. Hence, we have a separate acquisition function for each component, where $y^{(m)}$ is the evaluation of $f^{(m)}$:

where $\gamma^{(m)}_{y_{*}}({\boldsymbol{x}})=\frac{y^{(m)}_{*}-\mu^{(m)}_{t}({\boldsymbol{x}})}{\sigma^{(m)}_{t}({\boldsymbol{x}})}$. Analogously to the non-additive case, we sample $y^{(m)}_{*}$, separately for each function component. We select the final $x_{t}$ by choosing a sub-vector $x_{t}^{(m)}\in\arg\max_{{\boldsymbol{x}}^{(m)}\in A_{m}}\alpha^{(m)}_{t}({\boldsymbol{x}}^{(m)})$ and concatenating the components.

The Gumbel sampling from Section 3.1 directly extends to sampling $y^{(m)}_{*}$, approximately. We simply need to sample from the component-wise CDF $\widehat{\Pr}[{y}^{(m)}_{*}<z]=\prod_{{\boldsymbol{x}}\in\hat{\mathfrak{X}}}\Psi(\gamma^{(m)}_{y}({\boldsymbol{x}})))$, and use the same Gumbel approximation.

The algorithm for the additive max-value entropy search method (add-MES) is shown in Algorithm 2. The function Approx-MI does the pre-computation for approximating the mutual information in a similar way as in Algorithm 1, except that it only acts on the active dimensions in the $m$-th group.

Experiments

In this section, we probe the empirical performance of MES and add-MES on a variety of tasks. Here, MES-G denotes MES with $y_{*}$ sampled from the approximate Gumbel distribution, and MES-R denotes MES with $y_{*}$ computed by maximizing a sampled function represented by random features. Following (Hennig & Schuler, 2012; Hernández-Lobato et al., 2014), we adopt the zero mean function and non-isotropic squared exponential kernel as the prior for the GP. We compare to methods from the entropy search family, i.e., ES and PES, and to other popular Bayesian optimization methods including GP-UCB (denoted by UCB), PI, EI and EST. The parameter for GP-UCB was set according to Theorem 2 in (Srinivas et al., 2010); the parameter for PI was set to be the observation noise $\sigma$. For the functions with unknown GP hyper-parameters, every 10 iterations, we learn the GP hyper-parameters using the same approach as was used by PES (Hernández-Lobato et al., 2014). For the high dimensional tasks, we follow (Kandasamy et al., 2015) and sample the additive structure/GP parameters with the highest data likelihood when they are unknown. We evaluate performance according to the simple regret and inference regret as defined in Section 2.3. We used the open source Matlab implementation of PES, ES and EST (Hennig & Schuler, 2012; Hernández-Lobato et al., 2014; Wang et al., 2016). Our Matlab code and test functions are available at https://github.com/zi-w/Max-value-Entropy-Search/.

We begin with a comparison on synthetic functions sampled from a 3-dimensional GP, to probe our conjecture that MES is much more robust to the number of $y_{*}$ sampled to estimate the acquisition function than PES is to the number of $x_{*}$ samples. For PES, we sample 100 (PES 100), 10 (PES 10) and 1 (PES 1) argmaxes for the acquisition function. Similarly, we sample 100, 10, 1 $y_{*}$ values for MES-R and MES-G. We average the results on 100 functions sampled from the same Gaussian kernel with scale parameter $5.0$ and bandwidth parameter $0.0625$, and observation noise $\mathcal{N}(0,0.01^{2})$.

Figure 1 shows the simple and inference regrets. For both regret measures, PES is very sensitive to the the number of $x_{*}$ sampled for the acquisition function: 100 samples lead to much better results than 10 or 1. In contrast, both MES-G and MES-R perform competitively even with 1 or 10 samples. Overall, MES-G is slightly better than MES-R, and both MES methods performed better than other ES methods. MES methods performed better than all other methods with respect to simple regret. For inference regret, MES methods performed similarly to EST, and much better than all other methods including PES and ES.

In Table 1, we show the runtime of selecting the next input per iterationAll the timing experiments were run exclusively on an Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz. The function evaluation time is excluded. using GP-UCB, PI, EI, EST, ES, PES, MES-R and MES-G on the synthetic data with fixed GP hyper-parameters. For PES and MES-R, every $x_{*}$ or $y_{*}$ requires running an optimization sub-procedure, so their running time grows noticeably with the number of samples. MES-G avoids this optimization, and competes with the fastest methods EI and UCB.

In the following experiments, we set the number of $x_{*}$ sampled for PES to be 200, and the number of $y_{*}$ sampled for MES-R and MES-G to be 100 unless otherwise mentioned.

2 Optimization Test Functions

We test on three challenging optimization test functions: the 2-dimensional eggholder function, the 10-dimensional Shekel function and the 10-dimensional Michalewicz function. All of these functions have many local optima. We randomly sample 1000 points to learn a good GP hyper-parameter setting, and then run the BO methods with the same hyper-parameters. The first observation is the same for all methods. We repeat the experiments 10 times. The averaged simple regret is shown in the appendix, and the inference regret is shown in Table 2. On the 2-d eggholder function, PES was able to achieve better function values faster than all other methods, which verified the good performance of PES when sufficiently many $x_{*}$ are sampled. However, for higher-dimensional test functions, the 10-d Shekel and 10-d Michalewicz function, MES methods performed much better than PES and ES, and MES-G performed better than all other methods.

3 Tuning Hyper-parameters for Neural Networks

Next, we experiment with Levenberg-Marquardt optimization for training a 1-hidden-layer neural network. The 4 parameters we tune with BO are the number of neurons, the damping factor $\mu$, the $\mu$-decrease factor, and the $\mu$-increase factor. We test regression on the Boston housing dataset and classification on the breast cancer dataset (Bache & Lichman, 2013). The experiments are repeated 20 times, and the neural network’s weight initialization and all other parameters are set to be the same to ensure a fair comparison. Both of the datasets were randomly split into train/validation/test sets. We initialize the observation set to have 10 random function evaluations which were set to be the same across all the methods. The averaged simple regret for the regression L2-loss on the validation set of the Boston housing dataset is shown in Fig. 2(a), and the classification accuracy on the validation set of the breast cancer dataset is shown in Fig. 2(b). For the classification problem on the breast cancer dataset, MES-G, PES and UCB achieved a similar simple regret. On the Boston housing dataset, MES methods achieved a lower simple regret. We also show the inference regrets for both datasets in Table 3.

4 Active Learning for Robot Pushing

We use BO to do active learning for the pre-image learning problem for pushing (Kaelbling & Lozano-Pérez, 2017). The function we optimize takes as input the pushing action of the robot, and outputs the distance of the pushed object to the goal location. We use BO to minimize the function in order to find a good pre-image for pushing the object to the designated goal location. The first function we tested has a 3-dimensional input: robot location $(r_{x},r_{y})$ and pushing duration $t_{r}$. We initialize the observation size to be one, the same across all methods. The second function has a 4-dimensional input: robot location and angle $(r_{x},r_{y},r_{\theta})$, and pushing duration $t_{r}$. We initialize the observation to be 50 random points and set them the same for all the methods. We select 20 random goal locations for each function to test if BO can learn where to push for these locations. We show the simple regret in Fig. 4 and the inference regret in Table 4. MES methods performed on a par with or better than their competitors.

5 High Dimensional BO with Add-MES

In this section, we test our add-MES algorithm on high dimensional black-box function optimization problems. First we compare add-MES and add-GP-UCB (Kandasamy et al., 2015) on a set of synthetic additive functions with known additive structure and GP hyper-parameters. Each function component of the synthetic additive function is active on at most three input dimensions, and is sampled from a GP with zero mean and Gaussian kernel (bandwidth = $0.1$ and scale = $5$). For the parameter of add-GP-UCB, we follow (Kandasamy et al., 2015) and set $\beta^{(m)}_{t}=|A_{m}|\log 2t/5$. We set the number of $y^{(m)}_{*}$ sampled for each function component in add-MES-R and add-MES-G to be 1. We repeat each experiment for 50 times for each dimension setting. The results for simple regret are shown in Fig. 3. Add-MES methods perform much better than add-GP-UCB in terms of simple regret. Interestingly, add-MES-G works better in lower dimensional cases where $d=10,20,30$, while add-MES-R outperforms both add-MES-G and add-GP-UCB for higher dimensions where $d=50,100$. In general, MES-G tends to overestimate the maximum of the function because of the independence assumption, and MES-R tends to underestimate the maximum of the function because of the imperfect global optimization of the posterior function samples. We conjecture that MES-R is better for settings where exploitation is preferred over exploration (e.g., not too many local optima), and MES-G works better if exploration is preferred.

To further verify the performance of add-MES in high dimensional problems, we test on two real-world high dimensional experiments. One is a function that returns the distance between a goal location and two objects being pushed by a robot which has 14 parametersWe implemented the function in (Catto, 2011).. The other function returns the walking speed of a planar bipedal robot, with 25 parameters to tune (Westervelt et al., 2007). In Fig. 5, we show the simple regrets achieved by add-GP-UCB and add-MES. Add-MES methods performed competitively compared to add-GP-UCB on both tasks.

Conclusion

We proposed a new information-theoretic approach, max-value entropy search (MES), for optimizing expensive black-box functions. MES is competitive with or better than previous entropy search methods, but at a much lower computational cost. Via additive GPs, MES is adaptable to high-dimensional settings. We theoretically connected MES to other popular Bayesian optimization methods including entropy search, GP-UCB, PI, and EST, and showed a bound on the simple regret for a variant of MES. Empirically, MES performs well on a variety of tasks.

Acknowledgements

We thank Prof. Leslie Pack Kaelbling and Prof. Tomás Lozano-Pérez for discussions on active learning and Dr. William Huber for his solution to “Extreme Value Theory - Show: Normal to Gumbel” at stats.stackexchange.com, which leads to our Gumbel approximation in Section 3.1. 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 Related work

Our work is largely inspired by the entropy search (ES) methods (Hennig & Schuler, 2012; Hernández-Lobato et al., 2014), which established the information-theoretic view of Bayesian optimization by evaluating the inputs that are most informative to the $\operatorname{arg\,max}$ of the function we are optimizing.

Our work is also closely related to probability of improvement (PI) (Kushner, 1964), expected improvement (EI) (Moc̆kus, 1974), and the BO algorithms using upper confidence bound to direct the search (Auer, 2002; Kawaguchi et al., 2015, 2016), such as GP-UCB (Srinivas et al., 2010). In (Wang et al., 2016), it was pointed out that GP-UCB and PI are closely related by exchanging the parameters. Indeed, all these algorithms build in the heuristic that the next evaluation point needs to be likely to achieve the maximum function value or have high probability of improving the current evaluations, which in turn, may also give more information on the function optima like how ES methods queries. These connections become clear as stated in Section 3.1 of our paper.

Finding these points that may have good values in high dimensional space is, however, very challenging. In the past, high dimensional BO algorithms were developed under various assumptions such as the existence of a lower dimensional function structure (Djolonga et al., 2013; Wang et al., 2013), or an additive function structure where each component is only active on a lower manifold of the space (Li et al., 2016; Kandasamy et al., 2015). In this work, we show that our method also works well in high dimensions with the additive assumption made in (Kandasamy et al., 2015).

To sample the function maximum $y_{*}$, our first approach is to approximate the distribution for $y^{*}$ and then sample from that distribution. We use independent Gaussians to approximate the correlated $f({\boldsymbol{x}}),\forall{\boldsymbol{x}}\in\hat{\mathfrak{X}}$ where $\hat{\mathfrak{X}}$ is a discretization of the input search space $\mathfrak{X}$ (unless $\mathfrak{X}$ is discrete, in which case $\hat{\mathfrak{X}}=\mathfrak{X}$). A similar approach was adopted in (Wang et al., 2016). We can show that by assuming $\{f({\boldsymbol{x}})\}_{{\boldsymbol{x}}\in\hat{\mathfrak{X}}}$, our approximated distribution gives a distribution for an upperbound on $f({\boldsymbol{x}})$.

By the Slepian’s lemma, if the covariance $k_{t}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})\geq 0,\forall{\boldsymbol{x}},{\boldsymbol{x}}^{\prime}\in\hat{\mathfrak{X}}$, using the independent assumption with give us a distribution on the upperbound $\hat{y}_{*}$ of $f({\boldsymbol{x}})$, $\Pr[\hat{y}_{*}<y]=\prod_{{\boldsymbol{x}}\in\hat{\mathfrak{X}}}\Psi(\gamma_{y}({\boldsymbol{x}})))$.

We then use the Gumbel distribution to approximate the distribution for the maximum of the function values for $\hat{\mathfrak{X}}$, $\Pr[\hat{y}_{*}<y]=\prod_{{\boldsymbol{x}}\in\hat{\mathfrak{X}}}\Psi(\gamma_{y}({\boldsymbol{x}})))$. If for all ${\boldsymbol{x}}\in\hat{\mathfrak{X}}$, $f({\boldsymbol{x}})$ have the same mean and variance, the Gumbel approximation is in fact asymptotically correct by the Fisher-Tippett-Gnedenko theorem (Fisher, 1930).

In particular, for i.i.d. Gaussians, the limit distribution of the maximum of them belongs to the Gumbel distribution (Von Mises, 1936). Though the Fisher-Tippett-Gnedenko theorem does not hold for independent and differently distributed Gaussians, in practice we still find it useful in approximating $\Pr[\hat{y}_{*}<y]$. In Figure 6, we show an example of the result of the approximation for the distribution of the maximum of $f({\boldsymbol{x}})\sim GP(\mu_{t},k_{t})\forall{\boldsymbol{x}}\in\hat{\mathfrak{X}}$ given 50 observed data points randomly selected from a function sample from a GP with 0 mean and Gaussian kernel.

Appendix C Regret bounds

Based on the connection of MES to EST, we show the bound on the learning regret for MES with a point estimate for $\alpha(x)$. See 3.2

Before we continue to the proof, notice that if the function upper bound $\hat{y}_{*}$ is sampled using the approach described in Section 3.1 and $k_{t}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})\geq 0,\forall{\boldsymbol{x}},{\boldsymbol{x}}^{\prime}\in\hat{\mathfrak{X}}$, we may still get the regret guarantee by setting $y_{*}=\hat{y}_{*}$ (or $y_{*}=\hat{y}_{*}+\epsilon L$ if $\mathfrak{X}$ is continuous) since $\Pr[\max_{\hat{\mathfrak{X}}}\leq y]\geq\Pr[\hat{y}_{*}<y]$. Moreover, Theorem 3.2 assumes $y_{*}$ is sampled from a universal maximum distribution of functions from $GP(\mu,k)$, but it is not hard to see that if we have a distribution of maximums adapted from $GP(\mu_{t},k_{t})$, we can still get the same regret bound by setting $T^{\prime}=\sum_{i=1}^{T}\log_{w_{i}}\frac{\delta}{2\pi_{i}}$, where $w_{i}=F_{i}(f_{*})$ and $F_{i}$ corresponds to the maximum distribution at an iteration where $y_{*}>f_{*}$. Next we introduce a few lemmas and then prove Theorem 3.2.

Pick $\delta\in(0,1)$ and set $\zeta_{t}=(2\log(\frac{\pi_{t}}{2\delta}))^{\frac{1}{2}}$, where $\sum_{t=1}^{T}\pi_{t}^{-1}\leq 1$, $\pi_{t}>0$. Then, it holds that $\Pr[\mu_{t-1}({\boldsymbol{x}}_{t})-f({\boldsymbol{x}}_{t})\leq\zeta_{t}\sigma_{t-1}({\boldsymbol{x}}_{t}),\forall t\in[1,T]]\geq 1-\delta$.

(Theorem 3.2) By lemma 3.1 in our paper, we know that the theoretical results from EST (Wang et al., 2016) can be adapted to MES if $y_{*}\geq f_{*}$. The key question is when a sampled $y_{*}$ that can satisfy this condition. Because the cumulative density $w=F(f_{*})\in(0,1)$ and $y^{t}_{*}$ are independent samples from $F$, there exists at least one $y^{t}_{*}$ that satisfies $y^{t}_{*}>f_{*}$ with probability at least $1-w^{k_{i}}$ in $k_{i}$ iterations.

Let $T^{\prime}=\sum_{i=1}^{T}k_{i}$ be the total number of iterations. We split these iterations to $T$ parts where each part have $k_{i}$ iterations, $i=1,\cdots,T$. By union bound, with probability at least $1-\sum_{i=1}^{T}w^{k_{i}}$, in all the $T$ parts of iterations, we have at least one iteration $t_{i}$ which samples $y^{t_{i}}_{*}$ satisfying $y^{t_{i}}_{*}>f_{*},\forall i=1,\cdots,T$.

Let $\sum_{i=1}^{T}w^{k_{i}}=\frac{\delta}{2}$, we can set $k_{i}=\log_{w}\frac{\delta}{2\pi_{i}}$ for any $\sum_{i=1}^{T}(\pi_{i})^{-1}=1$. A convenient choice for $\pi_{i}$ is $\pi_{i}=\frac{\pi^{2}i^{2}}{6}$. Hence with probability at least $1-\frac{\delta}{2}$, there exist a sampled $y^{t_{i}}_{*}$ satisfying $y^{t_{i}}_{*}>f_{*},\forall i=1,\cdots,T$.

Now let $\zeta_{t_{i}}=(2\log\frac{\pi_{t_{i}}}{\delta})^{\frac{1}{2}}$. By Lemma C.1 and Lemma C.2, the immediate regret $r_{t_{i}}=f_{*}-f({\boldsymbol{x}}_{t_{i}})$ can be bounded as

As a result, our learning regret is bounded as

where $T^{\prime}=\sum_{i=1}^{T}k_{i}=\sum_{i=1}^{T}\log_{w}\frac{\delta}{2\pi_{i}}$ is the total number of iterations. ∎

At first sight, it might seem like MES with a point estimate does not have a converging rate as good as $EST$ or $GP-UCB$. However, notice that $\min_{{\boldsymbol{x}}\in\mathfrak{X}}\gamma_{y_{1}}({\boldsymbol{x}})<\min{{\boldsymbol{x}}\in\mathfrak{X}}\gamma_{y_{2}}({\boldsymbol{x}})$ if $y_{1}<y_{2}$, which decides the rate of convergence in Eq. 7. So if we use $y_{*}$ that is too large, the regret bound could be worse. If we use $y_{*}$ that is smaller than $f_{*}$, however, its value won’t count towards the learning regret in our proof, so it is also bad for the regret upper bound. With no principled way of setting $y_{*}$ since $f_{*}$ is unknown. Our regret bound in Theorem 3.2 is a randomized trade-off between sampling large and small $y_{*}$.

For the regret bound in add-GP-MES, it should follow add-GP-UCB. However, because of some technical problems in the proofs of the regret bound for add-GP-UCB, we haven’t been able to show a regret bound for add-GP-MES either. Nevertheless, from the experiments on high dimensional functions, the methods worked well in practice.

Appendix D Experiments

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

In Fig. 7, we show the simple regret comparing BO methods on the three challenging optimization test functions: the 2-D eggholder function, the 10-D Shekel function, and the 10-D Michalewicz function.

Choosing the additive decomposition

We follow the approach in (Kandasamy et al., 2015), and sample 10000 random decompositions (at most 2 dimensions in each group) and pick the one with the best data likelihood based on 500 data points uniformly randomly sampled from the search space. The decomposition setting was fixed for all the 500 iterations of BO for a fair comparison.