In the Introduction Bayesian Optimization GPyOpt we showed how GPyOpt can be used to solve optimization problems with some basic functionalities. The object
GPyOpt.methods.BayesianOptimization
is used to initialize the desired functionalities, such us the acquisition function, the initial design or the model. In some cases we want to have control over those objects and we may want to replace some element in the loop without having to integrate the new elements in the base code framework. This is now possible through the modular implementation of the package using the
GPyOpt.methods.ModularBayesianOptimization
class. In this notebook we are going to show how to use the backbone of GPyOpt to run a Bayesian optimization algorithm in which we will use our own acquisition function. In particular we are going to use the Expected Improvement integrated over the jitter parameter. That is
$$acqu_{IEI}(x;\{x_n,y_n\},\theta) = \int acqu_{EI}(x;\{x_n,y_n\},\theta,\psi) p(\psi;a,b)d\psi $$where $p(\psi;a,b)$ is, in this example, the distribution $Beta(a,b)$.
This acquisition is not available in GPyOpt, but we will implement and use in this notebook. The same can be done for other models, acquisition optimizers etc.
As usual, we start loading GPy and GPyOpt.
%pylab inline
import GPyOpt
import GPy
In this example we will use the Branin function as a test case.
# --- Function to optimize
func = GPyOpt.objective_examples.experiments2d.branin()
func.plot()
Because we are won't use the pre implemented wrapper, we need to create the classes for each element of the optimization. In total we need to create:
objective = GPyOpt.core.task.SingleObjective(func.f)
space = GPyOpt.Design_space(space =[{'name': 'var_1', 'type': 'continuous', 'domain': (-5,10)},
{'name': 'var_2', 'type': 'continuous', 'domain': (1,15)}])
model = GPyOpt.models.GPModel(optimize_restarts=5,verbose=False)
aquisition_optimizer = GPyOpt.optimization.AcquisitionOptimizer(space)
initial_design = GPyOpt.experiment_design.initial_design('random', space, 5)
from GPyOpt.acquisitions.base import AcquisitionBase
from GPyOpt.acquisitions.EI import AcquisitionEI
from numpy.random import beta
class jitter_integrated_EI(AcquisitionBase):
analytical_gradient_prediction = True
def __init__(self, model, space, optimizer=None, cost_withGradients=None, par_a=1, par_b=1, num_samples= 10):
super(jitter_integrated_EI, self).__init__(model, space, optimizer)
self.par_a = par_a
self.par_b = par_b
self.num_samples = num_samples
self.samples = beta(self.par_a,self.par_b,self.num_samples)
self.EI = AcquisitionEI(model, space, optimizer, cost_withGradients)
def acquisition_function(self,x):
acqu_x = np.zeros((x.shape[0],1))
for k in range(self.num_samples):
self.EI.jitter = self.samples[k]
acqu_x +=self.EI.acquisition_function(x)
return acqu_x/self.num_samples
def acquisition_function_withGradients(self,x):
acqu_x = np.zeros((x.shape[0],1))
acqu_x_grad = np.zeros(x.shape)
for k in range(self.num_samples):
self.EI.jitter = self.samples[k]
acqu_x_sample, acqu_x_grad_sample =self.EI.acquisition_function_withGradients(x)
acqu_x += acqu_x_sample
acqu_x_grad += acqu_x_grad_sample
return acqu_x/self.num_samples, acqu_x_grad/self.num_samples
Now we initialize the class for this acquisition and we plot the histogram of the used samples to integrate the acquisition.
acquisition = jitter_integrated_EI(model, space, optimizer=aquisition_optimizer, par_a=1, par_b=10, num_samples=200)
xx = plt.hist(acquisition.samples,bins=50)
# --- CHOOSE a collection method
evaluator = GPyOpt.core.evaluators.Sequential(acquisition)
With all the classes on place,including the one we have created for this example, we can now create the Bayesian optimization object.
bo = GPyOpt.methods.ModularBayesianOptimization(model, space, objective, acquisition, evaluator, initial_design)
And we run the optimization.
max_iter = 10
bo.run_optimization(max_iter = max_iter)
We plot the acquisition and the diagnostic plots.
bo.plot_acquisition()
bo.plot_convergence()