Gaussian Process Summer School, 13th September 2017. Written by Javier Gonzalez

The goal of this lab session is to illustrate the concepts seen during the tutorial in Gaussian processes for Global optimization. We will focus on two aspects of Bayesian Optimization (BO): (1) the choice of the model (2) the choice of the acquisition function.

The technical material associated to the methods used in this lab can be found in the slides of the tutorial.

`pip install gpyopt`

Some of the options of GPyOpt depend on other external packages: DIRECT, cma, pyDOE. Please be sure that this are installed if you want to use all the options. With everything installed, you are ready to start.

Now, just as in the previous lab, specify to include plots in the notebook and to import relevant libraries.

```
In [ ]:
```%pylab inline
import GPy
import GPyOpt
import numpy as np
from numpy.random import seed
seed(12345)

Before starting with the lab, remember that (BO) is an heuristic for global optimization of black-box functions. Let $f: {\mathcal X} \to R$ be a 'well behaved' continuous function defined on a compact subset ${\mathcal X} \subseteq R^d$. Our goal is to solve the global optimization problem of finding $$ x_{M} = \arg \min_{x \in {\mathcal X}} f(x). $$

We assume that $f$ is a *black-box* from which only perturbed evaluations of the type $y_i = f(x_i) + \epsilon_i$, with $\epsilon_i \sim\mathcal{N}(0,\sigma^2)$, are available. The goal is to find $x_M$ by minimizing the number of evaluations of $f$. To do this, we need to determine two crucial bits:

A

**Gaussian process**that will capture the our beliefs on $f$.An

**acquisition function**that based on the model will be useful to determine where to collect new evaluations of f.

Remember that every time a new data point is collected the model is updated and the acquisition function optimized again.

We start with a one-dimensional example. Consider here the Forrester function

$$f(x) =(6x-2)^2 \sin(12x-4),$$defined on the interval $[0, 1]$.

The minimum of this function is located at $x_{min}=0.78$. We assume that the evaluations of $f$ to are perturbed by zero-mean Gaussian noise with standard deviation 0.1. The Forrester function is part of the benchmark of functions of GPyOpt. To create the true function, the perturbed version and the boundaries of the problem you need to run the following cell.

```
In [ ]:
```f_true = GPyOpt.objective_examples.experiments1d.forrester() # true function object
f_sim = GPyOpt.objective_examples.experiments1d.forrester(sd=.1) # noisy version
bounds = [{'name': 'var_1', 'type': 'continuous', 'domain': (0,1)}] # problem constraints
f_objective = f_sim.f # objective function

To plot the true $f$, simply write:

```
In [ ]:
```f_true.plot()

```
In [ ]:
```n = 8
x = np.random.rand(n).reshape(n,1)
x

```
In [ ]:
```f_objective(x)

```
In [ ]:
``````
bounds
```

To use BO to solve this problem, we need to create a GPyOpt object in which we need to specify the following elements:

- The function to optimize.
- The box constraints of the problem.
- The model, by default we use a GP with a Mattern32 kernel.
- The acquisition function (and its parameters).

```
In [ ]:
```# Creation of the object that we will use to run BO.
seed(123)
myBopt = GPyOpt.methods.BayesianOptimization(f = f_objective, # function to optimize
domain = bounds, # box-constraints of the problem
acquisition_type ='EI') # acquisition = Expected improvement

```
In [ ]:
```GPyOpt.methods.BayesianOptimization?

```
In [ ]:
```myBopt.X

```
In [ ]:
```myBopt.Y

```
In [ ]:
```myBopt.model.model

```
In [ ]:
```# Run the optimization (may take a few senconds)
max_iter = 5 # evaluation budget
myBopt.run_optimization(max_iter) # run optimization

```
In [ ]:
```myBopt.plot_acquisition()

```
In [ ]:
```myBopt.plot_convergence()

```
In [ ]:
```myBopt.model.model

Use Bayesian optimization to find the minimum of the function $f(x)= x^2 + 10 \sin(x)$ in the interval [-10, 10].

1.1 Define the bounds of the problem, the function and check that it admits a numpy array of observations as input.

```
In [ ]:
``````
# Answer to 1.1 here:
```

```
In [ ]:
``````
# Answer to 1.2 here:
```

```
In [ ]:
``````
# Answer to 1.3 here:
```

```
In [ ]:
```from GPyOpt.acquisitions import AcquisitionEI, AcquisitionLCB, AcquisitionMPI

```
In [ ]:
```seed(1234)
n = 10
X = np.random.rand(n).reshape(n,1)
Y = f_objective(X)
m = GPy.models.GPRegression(X,Y)
m.optimize_restarts(10)
m.plot([0,1])

```
In [ ]:
```# GPyOpt model
model = GPyOpt.models.GPModel(optimize_restarts=5,verbose=False)
model.model = m
# GPyOpt design space
space = GPyOpt.Design_space(bounds)
# GPyOpt acquisition optimizer
aquisition_optimizer = GPyOpt.optimization.AcquisitionOptimizer(space)

Now we are ready to create the acquisitions.

```
In [ ]:
```acq_EI = AcquisitionEI(model,space,aquisition_optimizer,jitter=0)
acq_LCB = AcquisitionLCB(model,space,aquisition_optimizer,exploration_weight=2)
acq_MPI = AcquisitionMPI(model,space,aquisition_optimizer,jitter=0)

```
In [ ]:
```# Plot the three acquisition functions (factor 0.1 added in in the LCB for visualization)
X_grid = np.linspace(0,1,200)[:, None]
plt.figure(figsize=(10,6))
plt.title('Acquisition functions',size=25)
plt.plot(X_grid, - 0.1*acq_LCB.acquisition_function(X_grid),label='LCB')
plt.plot(X_grid, -acq_EI.acquisition_function(X_grid),label='EI')
plt.plot(X_grid, -acq_MPI.acquisition_function(X_grid),label='MPI')
plt.xlabel('x',size=15)
plt.ylabel('a(x)',size=15)
legend()

```
In [ ]:
``````
# Answer to 2.1 here:
```

```
In [ ]:
``````
# Answer to 2.2 here:
```

Consider the sixhumpcamel function defined as $$f(x_1,x_2) = \left(4-2.1x_1^2 + \frac{x_1^4}{3} \right)x_1^2 + x_1x_2 + (-4 +4x_2^2)x_2^2,$$

in $[-2,2]\times [-1,1]$. This function has two global minima, at $(0.0898,-0.7126)$ and $(-0.0898,0.7126)$. This function is also implemented in GPyOpt so, to load and visualize it simply run.

```
In [ ]:
```GPyOpt.objective_examples.experiments2d.sixhumpcamel().plot()
f_shc = GPyOpt.objective_examples.experiments2d.sixhumpcamel(sd = 0.1) # simulated version with some noise
f_shc_objective = f_shc.f
bounds_shc =[{'name': 'var_1', 'type': 'continuous', 'domain': f_shc.bounds[0]},
{'name': 'var_2', 'type': 'continuous', 'domain': f_shc.bounds[1]}]

```
In [ ]:
``````
# Answer to 3.1
```

3.2 In the three cases run the optimization for 30 iterations

```
In [ ]:
``````
# Answer to 3.2 here:
```

```
In [ ]:
``````
# Answer to 3.3 here:
```

```
In [ ]:
``````
# Answer to 3.4 here:
```