Bayesian optimisation with continuous and discrete parameters#
In this example, NUBO is used for sequential single-point optimisation with continuous and discrete parameters. Additionally to the bounds, a dictionary containing the dimensions as keys and the possible values as values have to be specified. The Hartmann6D synthetic test function acts as a substitute for a black-box objective function, such as an experiment or a simulation. We use the analytical acquisiton function UpperConfidenceBound with \(\beta = 1.96^2\) corresponding to the
95% confidence interval of the Gaussian distribution. We optimise this acquisition function with the L-BFGS-B algorithm with 5 starts to avoid getting stuck in a local maximum. The optimisation loop is run for 40 iterations and finds a solution close to the true optimum of -3.3224.
[1]:
import torch
from nubo.acquisition import ExpectedImprovement, UpperConfidenceBound
from nubo.models import GaussianProcess, fit_gp
from nubo.optimisation import single
from nubo.test_functions import Hartmann6D
from nubo.utils import gen_inputs
from gpytorch.likelihoods import GaussianLikelihood
# test function
func = Hartmann6D(minimise=False)
dims = 6
# specify bounds and discrete values
bounds = torch.tensor([[0., 0., 0., 0., 0., 0.], [1., 1., 1., 1., 1., 1.]])
discrete = {0: [0.2, 0.4, 0.6, 0.8], 4: [0.3, 0.6, 0.9]}
# training data
x_train = gen_inputs(num_points=dims*5,
num_dims=dims,
bounds=bounds)
y_train = func(x_train)
# Bayesian optimisation loop
iters = 40
for iter in range(iters):
# specify Gaussian process
likelihood = GaussianLikelihood()
gp = GaussianProcess(x_train, y_train, likelihood=likelihood)
# fit Gaussian process
fit_gp(x_train, y_train, gp=gp, likelihood=likelihood, lr=0.1, steps=200)
# specify acquisition function
# acq = ExpectedImprovement(gp=gp, y_best=torch.max(y_train))
acq = UpperConfidenceBound(gp=gp, beta=1.96**2)
# optimise acquisition function
x_new, _ = single(func=acq, method="L-BFGS-B", bounds=bounds, discrete=discrete, num_starts=5)
# evaluate new point
y_new = func(x_new)
# add to data
x_train = torch.vstack((x_train, x_new))
y_train = torch.hstack((y_train, y_new))
# print new best
if y_new > torch.max(y_train[:-1]):
print(f"New best at evaluation {len(y_train)}: \t Inputs: {x_new.numpy().reshape(dims).round(4)}, \t Outputs: {-y_new.numpy().round(4)}")
# results
best_iter = int(torch.argmax(y_train))
print(f"Evaluation: {best_iter+1} \t Solution: {-float(y_train[best_iter]):.4f}")
New best at evaluation 41: Inputs: [0.4 0.328 1. 0.323 0.3 1. ], Outputs: [-0.984]
New best at evaluation 42: Inputs: [0.2 0.3743 1. 0.3532 0.3 0.9297], Outputs: [-1.165]
New best at evaluation 43: Inputs: [0.2 0.3774 1. 0.3439 0.3 0.8455], Outputs: [-1.2931]
New best at evaluation 44: Inputs: [0.2 0.317 1. 0.3622 0.3 0.8081], Outputs: [-1.3629]
New best at evaluation 46: Inputs: [0.2 0.3075 1. 0.3253 0.3 0.7696], Outputs: [-1.4464]
New best at evaluation 47: Inputs: [0.2 0.274 1. 0.2838 0.3 0.7108], Outputs: [-1.5098]
New best at evaluation 48: Inputs: [0.2 0.2371 0.7806 0.2867 0.3 0.7069], Outputs: [-2.4508]
New best at evaluation 49: Inputs: [0.2 0.1603 0.7432 0.2521 0.3 0.7525], Outputs: [-2.526]
New best at evaluation 50: Inputs: [0.2 0.2488 0.6924 0.2493 0.3 0.766 ], Outputs: [-2.642]
New best at evaluation 51: Inputs: [0.2 0.1938 0.6262 0.2829 0.3 0.7344], Outputs: [-2.9647]
New best at evaluation 52: Inputs: [0.2 0.1336 0.5351 0.283 0.3 0.5795], Outputs: [-3.1219]
New best at evaluation 54: Inputs: [0.2 0.2051 0.4751 0.3168 0.3 0.6334], Outputs: [-3.2154]
New best at evaluation 56: Inputs: [0.2 0.1743 0.4045 0.2843 0.3 0.6907], Outputs: [-3.239]
New best at evaluation 60: Inputs: [0.2 0.1582 0.4587 0.2535 0.3 0.6505], Outputs: [-3.2954]
New best at evaluation 61: Inputs: [0.2 0.1474 0.4688 0.2743 0.3 0.6572], Outputs: [-3.315]
New best at evaluation 62: Inputs: [0.2 0.1487 0.4694 0.2756 0.3 0.6573], Outputs: [-3.3152]
Evaluation: 62 Solution: -3.3152