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.
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
.
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.
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.
[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.
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).
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 V
K
.
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
.