Acquisition functions¶
An acquisition function dictates the exploration policy of the optimizer. In practice, there are many different acquisition policies and which one is the best might depend on your situation. This package comes with several acquisition functions out-of-the-box, but it is also very easy to write your acquisition function. This notebook showcases how to do that using various examples.
Let’s begin with a few import statements and some function definitions we will use later.
[1]:
from bayes_opt import BayesianOptimization
from bayes_opt import acquisition
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
[2]:
def target(x):
return np.exp(-(x - 2)**2) + np.exp(-(x - 6)**2/10) + 1/ (x**2 + 1)
x = np.linspace(-2, 10, 1000).reshape(-1, 1)
y = target(x)
def posterior(optimizer, grid):
mu, sigma = optimizer._gp.predict(grid, return_std=True)
return mu, sigma
def plot_gp(optimizer, x, y):
acquisition_function_ = optimizer.acquisition_function
fig = plt.figure(figsize=(16, 10))
steps = len(optimizer.space)
fig.suptitle(
'Gaussian Process and Utility Function After {} Steps'.format(steps),
fontsize=30
)
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
axis = plt.subplot(gs[0])
acq = plt.subplot(gs[1])
x_obs = np.array([[res["params"]["x"]] for res in optimizer.res])
y_obs = np.array([res["target"] for res in optimizer.res])
acquisition_function_._fit_gp(optimizer._gp, optimizer._space)
mu, sigma = posterior(optimizer, x)
axis.plot(x, y, linewidth=3, label='Target')
axis.plot(x_obs.flatten(), y_obs, 'D', markersize=8, label=u'Observations', color='r')
axis.plot(x, mu, '--', color='k', label='Prediction')
axis.fill(np.concatenate([x, x[::-1]]),
np.concatenate([mu - 1.9600 * sigma, (mu + 1.9600 * sigma)[::-1]]),
alpha=.6, fc='c', ec='None', label='95% confidence interval')
axis.set_xlim((-2, 10))
axis.set_ylim((None, None))
axis.set_ylabel('f(x)', fontdict={'size':20})
axis.set_xlabel('x', fontdict={'size':20})
utility = -1 * acquisition_function_._get_acq(gp=optimizer._gp)(x)
x = x.flatten()
acq.plot(x, utility, label='Utility Function', color='purple')
acq.plot(x[np.argmax(utility)], np.max(utility), '*', markersize=15,
label=u'Next Best Guess', markerfacecolor='gold', markeredgecolor='k', markeredgewidth=1)
acq.set_xlim((-2, 10))
#acq.set_ylim((0, np.max(utility) + 0.5))
acq.set_ylabel('Utility', fontdict={'size':20})
acq.set_xlabel('x', fontdict={'size':20})
axis.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)
acq.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)
return fig, fig.axes
Most (simple) acquisition operate directly on the estimates of mean and standard deviation of the GP. For this, all it takes is to implement a .base_acq function, which combines the mean/std values into an acquisition policy. One very simple such function could simple disregard the std and always opt to investigate the maximum of the GP mean.
[3]:
class GreedyAcquisition(acquisition.AcquisitionFunction):
def __init__(self, random_state=None):
super().__init__(random_state)
def base_acq(self, mean, std):
return mean # disregard std
[4]:
acquisition_function = GreedyAcquisition(random_state=42)
optimizer = BayesianOptimization(target, {'x': (-2, 10)}, acquisition_function=acquisition_function, random_state=173)
optimizer.maximize(init_points=3, n_iter=7)
plot_gp(optimizer, x, y);
| iter | target | x |
-------------------------------------
| 1 | 0.5721 | -0.8809 |
| 2 | 0.9973 | 6.511 |
| 3 | 0.7985 | 4.267 |
| 4 | 0.9958 | 6.525 |
| 5 | 1.024 | 6.132 |
| 6 | 1.027 | 6.024 |
| 7 | 1.027 | 5.975 |
| 8 | 1.027 | 5.964 |
| 9 | 1.027 | 5.962 |
| 10 | 1.027 | 5.961 |
=====================================
As we can see, the greedy policy is too greedy and does not explore enough to find the maximum.
What if our acquisition function does not depend solely on mean and std? In that case, we can intervene at a deeper level, and overwrite the ._get_acq function. The built-in version of ._get_acq is a higher-order function which returns a function accepting an array-like x and performing the following: 1. Evaluate the GP mean \(\mu\) and std \(\sigma\) at points x. 2. If applicapable, evaluate the constraint fulfilment probability \(p_c\) at x. 3. Return
-1 * base_acq(mean, std) or -1 * base_acq(mean, std) * p_c
An example of such an acquisition function is Thompson Sampling. If you consider a Gaussian Process as a prior over functions, Thompson Sampling works by sampling a function from this prior and then selecting the argmax as next point of interest. This is a somewhat noisy version of the greedy acquisition, which will encourage some exploration.
While we usually find the argmax of the acquisition function through a combination of random sampling and gradient-based optimization, we will skip the gradient-based optimization here, as it is quite expensive and would require us to fix the multivariate normal. We can do this by additionally overwriting .suggest, specifically by changing the default argument of n_l_bfgs_b to 0.
[5]:
import warnings
class ThompsonSampling(acquisition.AcquisitionFunction):
def __init__(self, random_state=None):
super().__init__(random_state)
def base_acq(self, y_mean, y_cov):
assert y_cov.shape[0] == y_cov.shape[1], "y_cov must be a square matrix."
return self.random_state.multivariate_normal(y_mean, y_cov)
def _get_acq(self, gp, constraint=None):
if constraint is not None:
msg = (
f"Received constraints, but acquisition function {type(self)} "
+ "does not support constrained optimization."
)
raise acquisition.ConstraintNotSupportedError(msg)
# overwrite the base method since we require cov not std
dim = gp.X_train_.shape[1]
def acq(x):
x = x.reshape(-1, dim)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
mean, cov = gp.predict(x, return_cov=True)
return -1 * self.base_acq(mean, cov)
return acq
def suggest(self, gp, target_space, n_random=1_000, n_l_bfgs_b=0, fit_gp: bool = True):
# reduce n_random and n_l_bfgs_b to reduce the computational load
return super().suggest(gp, target_space, n_random, n_l_bfgs_b, fit_gp)
[6]:
acquisition_function = ThompsonSampling(random_state=13)
optimizer = BayesianOptimization(target, {'x': (-2, 10)}, acquisition_function=acquisition_function, random_state=173)
optimizer.maximize(init_points=3, n_iter=7)
fig, axs = plot_gp(optimizer, x, y);
print("Adding GP samples to the plot... this can take up to several minutes.")
x_reduced = np.linspace(-2, 10, 100).reshape(-1, 1)
for i in range(99):
y_sample = -1 * acquisition_function._get_acq(gp=optimizer._gp)(x_reduced)
y_sample_argmax = (x_reduced[np.argmax(y_sample)], np.max(y_sample))
axs[1].plot(x_reduced, y_sample, 'k', alpha=0.1, label='other samples from the GP' if i == 0 else None)
axs[1].plot(*y_sample_argmax, 'r.', alpha=0.4, markersize=10, label='other argmaxes' if i == 0 else None)
axs[1].legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.);
| iter | target | x |
-------------------------------------
| 1 | 0.5721 | -0.8809 |
| 2 | 0.9973 | 6.511 |
| 3 | 0.7985 | 4.267 |
| 4 | 0.7342 | 7.82 |
| 5 | 0.9227 | 0.8877 |
| 6 | 0.5426 | 8.523 |
| 7 | 0.2361 | 9.857 |
| 8 | 1.027 | 6.039 |
| 9 | 1.025 | 5.801 |
| 10 | 1.027 | 5.934 |
=====================================
Adding GP samples to the plot... this can take up to several minutes.
Let’s push things even further!
In some situations, you might be working in a setting where evaluations are not equally expensive. For example, let’s say that you’re running an expensive mechanical experiment, and performing the experiment is cheap, but adjusting its parameters is an expensive process. In particular, the further away it is from the current set-up, the more expensive it becomes. In such a situation, you could factor in the expense by running a specialized version of Expected Improvement (EI) called EIpu (EI per unit cost). Here we divide the EI by a specific cost function.
[7]:
class ExpectedImprovementPerUnitCost(acquisition.ExpectedImprovement):
def __init__(self, xi, exploration_decay=None, exploration_decay_delay=None, random_state=None) -> None:
super().__init__(xi, exploration_decay, exploration_decay_delay, random_state)
self.last_x = None
def cost(self, x):
if self.last_x is None:
return 1
return np.mean((self.last_x - np.atleast_2d(x))**2, axis=1) + 1.
def _get_acq(self, gp, constraint=None):
super_acq = super()._get_acq(gp, constraint)
acq = lambda x: super_acq(x) / self.cost(x)
return acq
def suggest(self, gp, target_space, n_random=10000, n_l_bfgs_b=10, fit_gp: bool = True):
# let's get the most recently evaluated point from the target_space
self.last_x = target_space.params[-1]
return super().suggest(gp, target_space, n_random, n_l_bfgs_b, fit_gp)
acquisition_function = ExpectedImprovementPerUnitCost(1e-4, random_state=42)
optimizer = BayesianOptimization(target, {'x': (-2, 10)}, acquisition_function=acquisition_function, random_state=173)
optimizer.maximize(init_points=3, n_iter=7)
fig, axs = plot_gp(optimizer, x, y)
ax2 = axs[1].twinx()
ax2.plot(x, acquisition_function.cost(x), label='Cost to access', color='orange')
axs[1].plot(x, -1 * super(acquisition.ExpectedImprovement, acquisition_function)._get_acq(gp=optimizer._gp)(x), 'k--', label='Expected Improvement')
lines, labels = axs[1].get_legend_handles_labels()
labels[0] = 'Expected Improvement per Unit Cost'
lines2, labels2 = ax2.get_legend_handles_labels()
axs[1].legend().remove()
ax2.legend(lines + lines2, labels + labels2, loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.);
| iter | target | x |
-------------------------------------
| 1 | 0.5721 | -0.8809 |
| 2 | 0.9973 | 6.511 |
| 3 | 0.7985 | 4.267 |
| 4 | 1.016 | 5.617 |
| 5 | 1.027 | 5.932 |
| 6 | 0.45 | 8.875 |
| 7 | 1.026 | 6.076 |
| 8 | 1.027 | 5.954 |
| 9 | 1.027 | 5.955 |
| 10 | 1.027 | 5.953 |
=====================================
[ ]:
[ ]: