Predictive Entropy Search for Efficient Global Optimization of Black-box Functions

2014·Arxiv

Abstract

Abstract

We propose a novel information-theoretic approach for Bayesian optimization called Predictive Entropy Search (PES). At each iteration, PES selects the next evaluation point that maximizes the expected information gained with respect to the global maximum. PES codifies this intractable acquisition function in terms of the expected reduction in the differential entropy of the predictive distribution. This reformulation allows PES to obtain approximations that are both more accurate and efficient than other alternatives such as Entropy Search (ES). Furthermore, PES can easily perform a fully Bayesian treatment of the model hyperparameters while ES cannot. We evaluate PES in both synthetic and real-world applications, including optimization problems in machine learning, finance, biotechnology, and robotics. We show that the increased accuracy of PES leads to significant gains in optimization performance.

1 Introduction

Bayesian optimization techniques form a successful approach for optimizing black-box functions [5]. The goal of these methods is to find the global maximizer of a nonlinear and generally nonconvex function f whose derivatives are unavailable. Furthermore, the evaluations of f are usually corrupted by noise and the process that queries f can be computationally or economically very expensive. To address these challenges, Bayesian optimization devotes additional effort to modeling the unknown function f and its behavior. These additional computations aim to minimize the number of evaluations that are needed to find the global optima.

Optimization problems are widespread in science and engineering and as a result so are Bayesian approaches to this problem. Bayesian optimization has successfully been used in robotics to adjust the parameters of a robot’s controller to maximize gait speed and smoothness [16] as well as parameter tuning for computer graphics [6]. Another example application in drug discovery is to find the chemical derivative of a particular molecule that best treats a given disease [20]. Finally, Bayesian optimization can also be used to find optimal hyper-parameter values for statistical [28] and machine learning techniques [24].

As described above, we are interested in finding the global maximizer of a function f over some bounded domain, typically . We assume that f(x) can only be evaluated via queries to a black-box that provides noisy outputs of the form . We note, however, that our framework can be extended to other non-Gaussian likelihoods. In this setting, we describe a sequential search algorithm that, after n iterations, proposes to evaluate f at some location . To make this decision the algorithm conditions on all previous observations . After N iterations the algorithm makes a final recommendation for the global maximizer of the latent function f.

We take a Bayesian approach to the problem described above and use a probabilistic model for the latent function f to guide the search and to select . In this work we use a zero-mean Gaussian

process (GP) prior for f [22]. This prior is specified by a positive-definite kernel function . Given any finite collection of points , the values of f at these points are jointly zero-mean Gaussian with covariance matrix , where . For the Gaussian likelihood described above, the vector of concatenated observations is also jointly Gaussian with zero-mean. Therefore, at any location x, the latent function f(x) conditioned on past observations is then Gaussian with marginal mean and variance given by

where is a vector of cross-covariance terms between x and .

Bayesian optimization techniques use the above predictive distribution to guide the search for the global maximizer . In particular, is used during the computation of an acquisition function that is optimized at each iteration to determine the next evaluation location . This process is shown in Algorithm 1. Intuitively, the acquisition function should be high in areas where the maxima is most likely to lie given the current data. However, should also encourage exploration of the search space to guarantee that the recommendation is a global optimum of f, not just a global optimum of the posterior mean. Several acquisition functions have been proposed in the literature. Some examples are the probability of improvement [14], the expected improvement [19] or upper confidence bounds [26]. Alternatively, one can combine several of these acquisition functions [10].

The acquisition functions described above are based on optimistic estimates of the latent function f which implicitly trade off between exploiting the posterior mean and exploring based on the uncertainty. We instead follow the approach described in [9] and aim to maximize the expected gain of information on the posterior distribution of the global maximizer . In Section 2 we derive a rearrangement of this acquisition function and a corresponding approximation that we call Predictive Entropy Search (PES). PES is more accurate than the approximation used in [9]. In Section 3 we empirically evaluate this claim on both synthetic and real-world problems and show that this leads to real gains in performance.

2 Predictive entropy search

We propose to follow the information-theoretic method for active data collection described in [17]. We are interested in maximizing information about the location of the global maximum, whose posterior distribution is . Our current information about can be measured in terms of the negative differential entropy of . Therefore, our strategy is to select which maximizes the expected reduction in this quantity. The corresponding acquisition function is

where represents the differential entropy of its argument and the expectation above is taken with respect to the posterior predictive distribution of y given x. The exact evaluation of (2) is infeasible in practice. The main difficulties are i) must be computed for many different values of x and y during the optimization of (2) and ii) the entropy computations themselves are not analytical. In practice, a direct evaluation of (2) is only possible after performing many approximations [9]. To avoid this, we follow the approach described in [11] by noting that (2) can be equivalently written as the mutual information between and y

given . Since the mutual information is a symmetric function, can be rewritten as

where is the posterior predictive distribution for y given the observed data and the location of the global maximizer of f. Intuitively, conditioning on the location pushes the posterior predictions up in locations around and down in regions away from . Note that, unlike the previous formulation, this objective is based on the entropies of predictive distributions, which are analytic or can be easily approximated, rather than on the entropies of distributions on whose approximation is more challenging.

The first term in (3) can be computed analytically using the posterior marginals for f(x) in (1), that is, , where we add to because y is obtained by adding Gaussian noise with variance to f(x). The second term, on the other hand, must be approximated. We first approximate the expectation in (3) by averaging over samples drawn approximately from . For each of these samples, we then approximate the corresponding entropy function using expectation propagation [18]. The code for all these operations is publicly available at http://tobediscloseduponacceptance.org.

2.1 Sampling from the posterior over global maxima

In this section we show how to approximately sample from the conditional distribution of the global maximizer given the observed data , that is,

If the domain X is restricted to some finite set of m points, the latent function f takes the form of an m-dimensional vector f. The probability that the ith element of f is optimal can then be written as. This suggests the following generative process: i) draw a sample from the posterior distribution and ii) return the index of the maximum element in the sampled vector. This process is known as Thompson sampling or probability matching when used as an arm-selection strategy in multi-armed bandits [8]. This same approach could be used for sampling the maximizer over a continuous domain X. At first glance this would require constructing an infinite-dimensional object representing the function f. To avoid this, one could sequentially construct f while it is being optimized. However, evaluating such an f would ultimately have cost where m is the number of function evaluations necessary to find the optimum. Instead, we propose to sample and optimize an analytic approximation to f. We will briefly derive this approximation below, but more detail is given in Appendix A.

Given a shift-invariant kernel k, Bochner’s theorem [4] asserts the existence of its Fourier dual s(w), which is equal to the spectral density of k. Letting be the associated normalized density, we can write the kernel as the expectation

where . Let denote an m-dimensional feature mapping where W and b consist of m stacked samples from p(w, b). The kernel k can then be approximated by the inner product of these features, . This approach was used by [21] as an approximation method in the context of kernel methods. The feature mapping allows us to approximate the Gaussian process prior for f with a linear model where is a standard Gaussian. By conditioning on , the posterior for is also multivariate Gaussian, where and . Let and be a random set of features and the corresponding posterior weights sampled both according to the generative process given above. They can then be used to construct the function , which is an approximate posterior sample of f—albeit one with a finite parameterization. We can then maximize this function to obtain , which is approximately distributed according to . Note that for early iterations when n < m, we can efficiently sample with cost using the method described in Appendix B.2 of [23]. This allows us to use a large number of features in .

2.2 Approximating the predictive entropy

We now show how to approximate in (3). Note that we can write the argument to H in this expression as . Here is the posterior distribution on f(x) given and the location of the global maximizer of f. When the likelihood p(y|f(x)) is Gaussian, we have that is analytically tractable since it is the predictive distribution of a Gaussian process. However, by further conditioning on the location of the global maximizer we are introducing additional constraints, namely that for all . These constraints make intractable. To circumvent this difficulty, we instead use the following simplified constraints:

C1. is a local maximum. This is achieved by letting and ensuring that is negative definite. We further assume that the non-diagonal elements of , denoted by , are known, for example they could all be zero. This simplifies the negative-definite constraint. We denote by C1.1 the constraint given by and . We denote by C1.2 the constraint that forces the elements of to be negative.

C2. is larger than past observations. We also assume that for all . However, we only observe noisily via . To avoid making inference on these latent function values, we approximate the above hard constraints with the soft constraint , where and is the largest seen so far.

C3. f(x) is smaller than . This simplified constraint only conditions on the given x rather than requiring for all .

We incorporate these simplified constraints into to approximate . This is achieved by multiplying with specific factors that encode the above constraints. In what follows we briefly show how to construct these factors; more detail is given in Appendix B.

Consider the latent variable . To incorporate constraint C1.1 we can condition on the data and on the “observations” given by the constraints and . Since f is distributed according to a GP, the joint distribution between z and these observations is multivariate Gaussian. The covariance between the noisy observations and the extra noise-free derivative observations can easily be computed [25]. The resulting conditional distribution is also multivariate Gaussian with mean and covariance . These computations are similar to those performed in (1). Constraints C1.2 and C2 can then be incorporated by writing

where is the cdf of a zero-mean Gaussian distribution with variance . The first new factor in this expression guarantees that , where we have marginalized out, and the second set of factors guarantees that the entries in are negative.

Later integrals that make use of , however, will not admit a closed-form expression. As a result we compute a Gaussian approximation q(z) to this distribution using Expectation Propagation (EP) [18]. The resulting algorithm is similar to the implementation of EP for binary classification with Gaussian processes [22]. EP approximates each non-Gaussian factor in (6) with a Gaussian factor whose mean and variance are and , respectively. The EP approximation can then be written as . Note that these computations have so far not depended on x, so we can compute once and store them for later use, where and .

We will now describe how to compute the predictive variance of some latent function value f(x) given these constraints. Let be a vector given by the concatenation of the values of the latent function at x and . The joint distribution between f, z, the evaluations collected so far and the derivative “observations” and is multivariate Gaussian. Using q(z), we then obtain the following approximation:

Implicitly we are assuming above that f depends on our observations and constraint C1.1, but is independent of C1.2 and C2 given z. The computations necessary to obtain and are similar to those used above and in (1). The required quantities are similar to the ones used by EP to make predictions in the Gaussian process binary classifier [22]. We can then incorporate C3 by multiplying with a factor that guarantees . The predictive distribution for f(x) given and all the constraints can be approximated as

where Z is a normalization constant. The variance of the right hand size of (8) is given by

where , and and are the standard Gaussian density function and cdf, respectively. By further approximating (8) by a Gaussian distribution with the same mean and variance we can write the entropy as .

The computation of (9) can be numerically unstable when s is very close to zero. This occurs when is very similar to . To avoid these numerical problems, we multiply by the largest that guarantees that . This can be understood as slightly reducing the amount of dependence between f(x) and when x is very close to . Finally, fixing to be zero can also produce poor predictions when the actual f does not satisfy this constraint. To avoid this, we instead fix this quantity to , where is the ith sample function optimized in Section 2.1 to sample .

2.3 Hyperparameter learning and the PES acquisition function

We now show how the previous approximations are integrated to compute the acquisition function used by predictive entropy search (PES). This acquisition function performs a formal treatment of the hyperparameters. Let denote a vector of hyperparameters which includes any kernel parameters as well as the noise variance . Let denote the posterior distribution over these parameters where is a hyperprior and is the GP marginal likelihood. For a fully Bayesian treatment of we must marginalize the acquisition function (3) with respect to this posterior. The corresponding integral has no analytic expression and must be approximated using Monte Carlo. This approach is also taken in [24].

We draw M samples from using slice sampling [27]. Let denote a sampled global maximizer drawn from as described in Section 2.1. Furthermore, let and denote the predictive variances computed as described in Section 2.2 when the model hyperparameters are fixed to . We then write the marginalized acquisition function as

Note that PES is effectively marginalizing the original acquisition function (2) over . This is a significant advantage with respect to other methods that optimize the same information-theoretic acquisition function but do not marginalize over the hyper-parameters. For example, the approach of [9] approximates (2) only for fixed . The resulting approximation is computationally very expensive and recomputing it to average over multiple samples from is infeasible in practice.

Algorithm 2 shows pseudo-code for computing the PES acquisition function. Note that most of the computations necessary for evaluating (10) can be done independently of the input x, as noted in the pseudo-code. This initial cost is dominated by a matrix inversion necessary to pre-compute V for each hyperparameter sample. The resulting complexity is . This cost can be reduced to by ignoring the derivative observations imposed on by constraint C1.1. Nevertheless, in the problems that we consider d is very small (less than 20). After these precomputations are done, the evaluation of (10) is .

3 Experiments

In our experiments, we use Gaussian process priors for f with squared-exponential kernels . The corresponding spectral density is zero-mean Gaus- sian with covariance given by and normalizing constant . The model hyperpa- rameters are . We use broad, uninformative Gamma hyperpriors.

Figure 1: Comparison of different estimates of the objective function given by (2). Left, ground truth obtained by the rejection sampling method RS. Middle, approximation produced by the ES method. Right, approximation produced by the proposed PES method. These plots show that the PES objective is much more similar to the RS ground truth than the ES objective.

First, we analyze the accuracy of PES in the task of approximating the differential entropy (2). We compare the PES approximation (10), with the approximation used by the entropy search (ES) method [9]. We also compare with the ground truth for (2) obtained using a rejection sampling (RS) algorithm based on (3). For this experiment we generate the data using an objective function f sampled from the Gaussian process prior. The domain X of f is fixed to be and data are generated using , and . To compute (10) we avoid sampling the hyperparameters and use the known values directly. We further fix M = 200 and m = 1000.

The ground truth rejection sampling scheme works as follows. First, X is discretized using a uniform grid. The expectation with respect to in (3) is then approximated using sampling. For this, we sample by evaluating a random sample from on each grid cell and then selecting the cell with highest value. Given , we then approximate by rejection sampling. We draw samples from and reject those whose corresponding grid cell with highest value is not . Finally, we approximate by first, adding zero-mean Gaussian noise with variance to the the evaluations at x of the functions not rejected during the previous step and second, we estimate the differential entropy of the resulting samples using kernels [1].

Figure 1 shows the objective functions produced by RS, ES and PES for a particular with 10 measurements whose locations are selected uniformly at random in . The locations of the collected measurements are displayed with an “x” in the plots. The particular objective function used to generate the measurements in is displayed in the left part of Figure 2. The plots in Figure 1 show that the PES approximation to (2) is more similar to the ground truth given by RS than the approximation produced by ES. In this figure we also see a discrepancy between RS and PES at locations near x = (0.572, 0.687). This difference is an artifact of the discretization used in RS. By zooming in and drawing many more samples we would see the same behavior in both plots.

We now evaluate the performance of PES in the task of finding the optimum of synthetic black-box objective functions. For this, we reproduce the within-model comparison experiment described in [9]. In this experiment we optimize objective functions defined in the 2-dimensional unit domain . Each objective function is generated by first sampling 1024 function values from the GP prior assumed by PES, using the same and as in the previous experiment. The objective function is then given by the resulting GP posterior mean. We generated a total of 1000 objective functions by following this procedure. The left plot in Figure 2 shows an example function.

In these experiments we compared the performance of PES with that of ES [9] and expected improvement (EI) [13], a widely used acquisition function in the Bayesian optimization literature. We again assume that the optimal hyper-parameter values are known to all methods. Predictive performance is then measured in terms of the immediate regret (IR) , where is the known location of the global maximum and is the recommendation of each algorithm had we stopped at step n—for all methods this is given by the maximizer of the posterior mean. The right plot in Figure 2 shows the decimal logarithm of the median of the IR obtained by each method across the 1000 different objective functions. Confidence bands equal to one standard deviation are obtained using the bootstrap method. Note that while averaging these results is also interesting, corresponding to the expected performance averaged over the prior, here we report the median IR

Figure 2: Left, example of objective functions f. Right, median of the immediate regret (IR) for the methods PES, ES and EI in the experiments with synthetic objective functions.

Figure 3: Median of the immediate regret (IR) for the methods EI, ES, PES and PES-NB in the experiments with well-known synthetic benchmark functions.

because the empirical distribution of IR values is very heavy-tailed. In this case, the median is more representative of the exact location of the bulk of the data. These results indicate that the best method in this setting is PES, which significantly outperforms ES and EI. The plot also shows that in this case ES is significantly better than EI.

We perform another series of experiments in which we optimize well-known synthetic benchmark functions including a mixture of cosines [2] and Branin-Hoo (both functions defined in ) as well as the Hartmann-6 (defined in ) [15]. In all instances, we fix the measurement noise to . For both PES and EI we marginalize the hyperparameters using the approach described in Section 2.3. ES, by contrast, cannot average its approximation of (2) over the posterior on . Instead, ES works by fixing to an estimate of its posterior mean (obtained using slice sampling) [27]. To evaluate the gains produced by the fully Bayesian treatment of in PES, we also compare with a version of PES (PES-NB) which performs the same non-Bayesian (NB) treatment of as ES. In PES-NB we use a single fixed hyperparameter as in previous sections with value given by the posterior mean of . All the methods are initialized with three random measurements collected using latin hypercube sampling [5].

The plots in Figure 3 show the median IR obtained by each method on each function across 250 random initializations. Overall, PES is better than PES-NB and ES. Furthermore, PES-NB is also significantly better than ES in most of the cases. These results show that the fully Bayesian treatment of in PES is advantageous and that PES can produce better approximations than ES. Note that PES performs better than EI in the Branin and cosines functions, while EI is significantly better on the Hartmann problem. This appears to be due to the fact that entropy-based strategies explore more aggressively which in higher-dimensional spaces takes more iterations. The Hartmann problem, however, is a relatively simple problem and as a result the comparatively more greedy behavior of EI does not result in significant adverse consequences. Note that the synthetic functions optimized in the previous experiment were much more multimodal that the ones considered here.

3.1 Experiments with real-world functions

We finally optimize different real-world cost functions. The first one (NNet) returns the predictive accuracy of a neural network on a random train/test partition of the Boston Housing dataset [3].

Figure 4: Median of the immediate regret (IR) for the methods PES, PES-NB, ES and EI in the experiments with non-analytic real-world cost functions.

The variables to optimize are the weight-decay parameter and the number of training iterations for the neural network. The second function (Hydrogen) returns the amount of hydrogen production of a particular bacteria in terms of the PH and Nitrogen levels of the growth medium [7]. The third one (Portfolio) returns the ratio of the mean and the standard deviation (the Sharpe ratio) of the 1-year ahead returns generated by simulations from a multivariate time-series model that is adjusted to the daily returns of stocks AXP, BA and HD. The time-series model is formed by univariate GARCH models connected with a Student’s t copula [12]. These three functions (NNet, Hydrogen and Portfolio) have as domain . Furthermore, in these examples, the ground truth function that we want to optimize is unknown and is only available through noisy measurements. To obtain a ground truth, we approximate each cost function as the predictive distribution of a GP that is adjusted to data sampled from the original function (1000 uniform samples for NNet and Portfolio and all the available data for Hydrogen [7]). Finally, we also consider another real-world function that returns the walking speed of a bipedal robot [29]. This function is defined in and its inputs are the parameters of the robot’s controller. In this case the ground truth function is noiseless and can be exactly evaluated through expensive numerical simulation. We consider two versions of this problem (Walker A) with zero-mean, additive noise of and (Walker B) with .

Figure 4 shows the median IR values obtained by each method on each function across 250 random initializations, except in Hydrogen where we used 500 due to its higher level of noise. Overall, PES, ES and PES-NB perform similarly in NNet, Hydrogen and Portfolio. EI performs rather poorly in these first three functions. This method seems to make excessively greedy decisions and fails to explore the search space enough. This strategy seems to be advantageous in Walker A, where EI obtains the best results. By contrast, PES, ES and PES-NB tend to explore more in this latter dataset. This leads to worse results than those of EI. Nevertheless, PES is significantly better than PES-NB and ES in both Walker datasets and better than EI in the noisier Walker B. In this case, the fully Bayesian treatment of hyper-parameters performed by PES produces improvements in performance.

4 Conclusions

We have proposed a novel information-theoretic approach for Bayesian optimization. Our method, predictive entropy search (PES), greedily maximizes the amount of one-step information on the location of the global maximum using its posterior differential entropy. Since this objective function is intractable, PES approximates the original objective using a reparameterization that measures entropy in the posterior predictive distribution of the function evaluations. PES produces more accurate approximations than Entropy Search (ES), a method based on the original, non-transformed acquisition function. Furthermore, PES can easily marginalize its approximation with respect to the posterior distribution of its hyper-parameters, while ES cannot. Experiments with synthetic and real-world functions show that PES often outperforms ES in terms of immediate regret. In these experiments, we also observe that PES often produces better results than expected improvement (EI), a popular heuristic for Bayesian optimization. EI often seems to make excessively greedy decisions, while PES tends to explore more. As a result, EI seems to perform better for simple objective functions while often getting stuck with noisier objectives or for functions with many modes.

References

[1] I. Ahmad and P.-E. Lin. A nonparametric estimation of the entropy for absolutely continuous distributions. IEEE Transactions on Information Theory, 22(3):372–375, 1976.

[2] B. S. Anderson, A. W. Moore, and D. Cohn. A nonparametric approach to noisy and costly optimization. In ICML, pages 17–24, 2000.

[3] K. Bache and M. Lichman. UCI machine learning repository, 2013.

[4] S. Bochner. Lectures on Fourier integrals. Princeton University Press, 1959.

[5] E. Brochu, V. M. Cora, and N. de Freitas. A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. Technical Report UBC TR-2009-23 and arXiv:1012.2599v1, Dept. of Computer Science, University of British Columbia, 2009.

[6] E. Brochu, N. de Freitas, and A. Ghosh. Active preference learning with discrete choice data. In NIPS, pages 409–416, 2007.

[7] E. H. Burrows, W.-K. Wong, X. Fern, F. W. R. Chaplen, and R. L. Ely. Optimization of ph and nitrogen for enhanced hydrogen production by synechocystis sp. pcc 6803 via statistical and machine learning methods. Biotechnology Progress, 25(4):1009–1017, 2009.

[8] O. Chapelle and L. Li. An empirical evaluation of Thompson sampling. In NIPS, pages 2249–2257, 2011.

[9] P. Hennig and C. J. Schuler. Entropy search for information-efficient global optimization. Journal of Machine Learning Research, 13, 2012.

[10] M. W. Hoffman, E. Brochu, and N. de Freitas. Portfolio allocation for Bayesian optimization. In UAI, pages 327–336, 2011.

[11] N. Houlsby, J. M. Hern´andez-lobato, F. Huszar, and Z. Ghahramani. Collaborative gaussian processes for preference learning. In NIPS, pages 2096–2104. 2012.

[12] E. Jondeau and M. Rockinger. The copula-garch model of conditional dependencies: An international stock market application. Journal of international money and finance, 25(5):827–853, 2006.

[13] D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of expensive black-box func- tions. Journal of Global optimization, 13(4):455–492, 1998.

[14] H. Kushner. A new method of locating the maximum of an arbitrary multipeak curve in the presence of noise. Journal of Basic Engineering, 86, 1964.

[15] D. Lizotte. Practical Bayesian Optimization. PhD thesis, University of Alberta, Canada, 2008.

[16] D. Lizotte, T. Wang, M. Bowling, and D. Schuurmans. Automatic gait optimization with Gaussian process regression. In IJCAI, pages 944–949, 2007.

[17] D. J. MacKay. Information-based objective functions for active data selection. Neural Computation, 4(4):590–604, 1992.

[18] T. P. Minka. A family of algorithms for approximate Bayesian inference. PhD thesis, Massachusetts Institute of Technology, 2001.

[19] J. Moˇckus, V. Tiesis, and A. ˇZilinskas. The application of bayesian methods for seeking the extremum. In L. Dixon and G. Szego, editors, Toward Global Optimization, volume 2. Elsevier, 1978.

[20] D. M. Negoescu, P. I. Frazier, and W. B. Powell. The knowledge-gradient algorithm for sequencing experiments in drug discovery. INFORMS Journal on Computing, 23(3):346–363, 2011.

[21] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In NIPS, pages 1177–1184, 2007.

[22] C. E. Rasmussen and C. K. Williams. Gaussian processes for machine learning. The MIT Press, 2006.

[23] M. W. Seeger. Bayesian inference and optimal design for the sparse linear model. Journal of Machine Learning Research, 9:759–813, 2008.

[24] J. Snoek, H. Larochelle, and R. P. Adams. Practical Bayesian optimization of machine learning algo- rithms. In NIPS, pages 2960–2968, 2012.

[25] E. Solak, R. Murray-smith, W. E. Leithead, D. J. Leith, and C. E. Rasmussen. Derivative observations in gaussian process models of dynamic systems. In NIPS, pages 1057–1064. 2003.

[26] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In ICML, pages 1015–1022, 2010.

[27] J. Vanhatalo, J. Riihim¨aki, J. Hartikainen, P. Jyl¨anki, V. Tolvanen, and A. Vehtari. Bayesian modeling with gaussian processes using the matlab toolbox gpstuff (v3.3). CoRR, abs/1206.5754, 2012.

[28] Z. Wang, S. Mohamed, and N. de Freitas. Adaptive Hamiltonian and Riemann Monte Carlo samplers. In ICML, 2013.

[29] E. Westervelt and J. Grizzle. Feedback Control of Dynamic Bipedal Robot Locomotion. Control and Automation Series. CRC PressINC, 2007.

A Details on approximating GP sample paths

In this section we give further details about the approach used in Section 2.1 to approximate a GP using random features. These random features can be used to approximate sample paths from the GP posterior. By optimizing these sample paths we obtain posterior samples over the global maxima . We derive in more detail the kernel approximation from (5). Formally, the theorem of [4] states Theorem 1 (Bochner’s theorem). A continuous, shift-invariant kernel is positive definite if and only if it is the Fourier transform of a non-negative, finite measure.

As a result given some kernel there must exist an associated density s(w), known as its spectral density, which is the Fourier dual of k. This can be written as

Further, we can treat this measure as a probability density where is the normalizing constant. Consequently, the kernel can be written as

and due to the symmetry of p(w) [see 22] we can write the expectation as

We can then note thatfor any constant offset . As a result, for b uniformly distributed between 0 and we can write

The last equality can be derived from the sum of angles formula, which leads to the identity: . Finally, we can average over m weights and phases

where and are stacked versions of the original random variables. The resulting quantity has the same expectation but results in a lower variance estimator. If we let denote a random m-dimensional feature generated by this model we can also write the kernel as .

We now briefly show the equivalence between a Bayesian linear model using random features and a GP with kernel k. Consider now a linear model where has a standard Gaussian distribution and observations of the form . The posterior of given is also be normal N(m, V) where

and where and consist of the stacked features and observations respectively. We can also easily write the predictive distribution over f evaluated at a test point x, which is Gaussian distributed with mean and variance given by

By a simple application of the matrix-inversion lemma these quantities can be rewritten in terms which only make use of the inner products between features,

the expectations of which are equivalent to the kernel k and we obtain the same expressions as that in (1).

B Details on approximating the predictive variance

We now provide further details on approximating the predictive variance of inputs x given the position of the global optimizer . In particular we include all steps omitted in the presentation of Section 2.2.

B.1 Incorporating the analytic latent constraints (C1.1)

We first turn to the random variables

Here c contains the random variables that we will condition on in order to enforce constraint C1.1. Given the input locations x and we can construct a kernel matrix K containing the covariance evaluated on the stacked vector [z; c]. We again refer to [25] in constructing this matrix which includes derivative observations, the computations of which are tedious but not overly complicated. Note also that the portions of K which correspond to will have an additional due to the observation noise. Next let , and denote the corresponding diagonal and off-diagonal blocks of the kernel matrix. We can now condition on the observed values of c to write

where mc and VK.

B.2 Incorporating the non-analytic latent constraints (C1.2 and C2)

The additional constraints C1.2 and C2 can be introduced explicitly as in (6), which takes the form of a single Gaussian factor and d + 1 non-Gaussian factors

We approximate this distribution using a single multivariate Gaussian q(z) where each non-Gaussian factor is replaced by a Gaussian approximation such that

where this approximation is parameterized by and . The parameters of the approximate factors are combined to form the vector and the diagonal matrix .

To compute the approximate factors we use expectation propagation (EP). EP is a procedure that starts from some initial values for the approximate factors and iteratively refines these quantities; here we initialize and which corresponds to and . At each iteration, for every factor i, we remove the contribution of the ith approximate factor to form the cavity distribution . Given the independent factors we consider here we can focus on each individual component separately with mean and variance

Let denote the tilted distribution where the ith approximate factor has been replaced by the corresponding real factor. EP proceeds by finding the approximation that minimizes the KL-divergence where is restricted to be Gaussian. This amounts to matching the first two moments. Finally, by removing the influence of the cavity distribution and setting we can update the approximate factors. This can be performed using the same procedure which forms the cavity distribution.

For both sets of constraints used in this work the moments can easily be obtained by computing the normalizing constant and using the following identities:

For the factors corresponding to constraints on the diagonal Hessian, i.e. where 0], the distribution is also simply a truncated Gaussian. Given these moments we can remove the contribution of the cavity distribution as above and write

For the final soft-maximum constraint, , the moments can be calculated in a similar fashion. Using the same procedure as above we arrive at very similar updates:

B.3 Incorporating the prediction constraint (C3)

Given some test input x we now turn to the problem of making predictions about f(x). We again note that both the “prior” terms and the EP factors, and , are independent of x and can be precomputed once for later use at prediction time.

Let be a vector given by the concatenation of the latent function at x and . The distribution for f given the first two constraints can be written as

By writing p(f|z, c) above we are assuming that f is independent of C1.2 and C2 given z and as a result the above is simply an integral over the product of two Gaussians. Let be the cross-covariance matrix evaluated between f and [z; c] and the covariance matrix associated with f. The posterior above will then be Gaussian with mean and variance

where is a block-diagonal matrix where the first block is zero and the second is (note this matrix is also diagonal since is diagonal). Finally, these values can be plugged into (8–9) in order to arrive at .

designed for accessibility and to further open science