# practical implementation detail of Bayesian Optimization

I'm giving Bayesian Optimization a go, following Snoek, Larochelle, and Adams [http://arxiv.org/pdf/1206.2944.pdf], using GPML [http://www.gaussianprocess.org/gpml/code/matlab/doc/]. I've implemented the Expected Improvement acquisition function described on page 3, and I'm assuming I'm correct that to decide where to next query my objective I should take the $\bf{x}$ that maximizes:

$a_{EI}(\bf{x}; (\bf{x}_n,y_n,\theta))$

But I can't seem to find guidance on what set of candidates $\bf{x}$'s to consider. Theoretically, I'd like to find the best $\bf{x}$ over the entire domain, and the paper is written in such a way that seems to suggest this is possible ("[EI] also has closed form under the Gaussian process"). But as a practical matter, I need to calculate the posterior predictive means and variances at any $\bf{x}^*$ I might consider before I can calculate $a_{EI}(\bf{x}^*)$ and while these posteriors have a closed form, I still need to calculate them using matrix algebra, so I can't see a way to get around picking a bunch of $\bf{x}$'s.

The question: what is a practical method for choosing the large (medium? small?) set of candidate $\bf{x}$'s over which I maximize EI (or any other acquisition function)? (Is this in the paper somewhere and I just missed it?)

At the moment, I'm just taking my current set ${x_i}$, sampling it with replacement 2000 times, and then adding some Gaussian noise to each point. Seems okay, I guess.

• I haven't read this paper yet, but have you considered sampling from the domain of interest using a Latin Hypercube Design? en.wikipedia.org/wiki/Latin_hypercube_sampling – RustyStatistician Feb 5 '16 at 17:14
• An alternative would be to grid the domain (if its possible) and then evaluate every point on that gridded domain. – RustyStatistician Feb 5 '16 at 17:15
• These are both sensible suggestions, thanks! Don't know much about Latin hypercubes, but regular gridded domains would mean Toeplitz and Kronecker algebra could be used, which would make things nice and efficient even with a very large grid. – stackoverflax Feb 6 '16 at 17:36

The norm is to use any global optimizer you like. The problem is that the EI surface is highly multi-modal and disconnected; optimizing this acquisition function is a nontrivial problem in itself.

A common choice that I have seen in various papers is the DIRECT algorithm; sometimes I've seen CMA-ES which is a state-of-the-art method in nonlinear optimization. In my experience for other forms of optimization, MCS (Multi-Level Coordinate Search) tends to work relatively well. You can find a review of derivative-free global optimizers here:

• Rios and Sahinidis, "Derivative-free optimization: a review of algorithms and comparison of software implementations", Journal of Global Optimization (2013).

By the way, the EI is analytical so if you want you can also compute its gradient to guide the optimization, but this is not necessary. An effective technique is to run a global optimizer first to find promising solutions and then run a local optimizer to refine it (e.g., a quasi-Newton method such as BFGS, that is fminunc in MATLAB; or fmincon if you have constraints).

Finally, if speed of the optimization of the acquisition function is a factor (which is not the "traditional" BO scenario), I have found decent results by starting with a Latin Hypercube design or a quasi-random Sobol sequence design, then refined with a few steps of a local optimizer from the best point(s); see also @user777 comment. Since this is not the standard BO scenario, I don't have any specific reference that actually uses this method.

Examples of papers that refer to DIRECT or CMA-ES:

• Calandra, R., Seyfarth, A., Peters, J., & Deisenroth, M. P. (2015). Bayesian optimization for learning gaits under uncertainty. Annals of Mathematics and Artificial Intelligence, 1-19 (link).
• Mahendran, N., Wang, Z., Hamze, F., & Freitas, N. D. (2012). Adaptive MCMC with Bayesian optimization. In International Conference on Artificial Intelligence and Statistics (pp. 751-760) (link).
• Gunter, T., Osborne, M. A., Garnett, R., Hennig, P., & Roberts, S. J. (2014). Sampling for inference in probabilistic models with fast Bayesian quadrature. In Advances in neural information processing systems (pp. 2789-2797) (link).

You can just google "Bayesian optimization" + the desired global optimization algorithm, and you'll find a bunch of papers. Plus, in pretty much every other paper about BO you would find a sentence such as:

[...] BO usually requires an auxiliary global optimizer in each iteration to optimize the acquisition function. It is customary in the BO literature to use DIvided RECTangles (DIRECT) to accomplish such a task. Other global optimization algorithms like CMA-ES could also be applied.

• This is actually kind of surprising to me! Can you point me to a representative Bayesian optimization paper that you have in mind that uses DIRECT or CMA-ES? Thanks. – stackoverflax Mar 2 '16 at 12:35
• Why is it surprising? This is the norm -- you'll find references to DIRECT or other global optimizers in pretty much every BO paper. It is probably so well known in the community that a few papers don't even bother to mention -- it's just given for granted. I added a few references in my main comment above. – lacerbi Mar 2 '16 at 14:27
• This is not necessarily a good solution, but I've found it can be cheaper to just evaluate EI at a set of points sampled using Latin Hypercubes in cases where you only need to be nearby the minimum but not necessarily on top of it. – Reinstate Monica Mar 2 '16 at 14:56
• @user777: Yes, if speed is at stake I've been using both LH and Sobol quasi-random sequences as initial design (finding a slight advantage with the latter, but it might be problem-dependent), and then running a local optimizer such as BFGS from the best point(s). I will add this to the main comment. – lacerbi Mar 2 '16 at 15:04
• One way to justify the ad hoc nature of the LHS approach is that finding the minimum of a stochastic function (surface) is unnecessary because the error in the estimation of the minimum will swamp minute gains in precision. This is a very good answer, though. I'm glad someone else here cares about BO. :-) – Reinstate Monica Mar 2 '16 at 17:12