Parallel multi-point Bayesian optimisation with fixed base samples#
This notebook shows how NUBO can perform parallel multi-point optimisation with fixed base samples. This enables the use of deterministic optimisers, such as L-BFGS-B and SLSQP, and parallel and constrained optimisation.
In the first example below, NUBO is used to perform multi-point optimisation that allows the candidates to be evaluated from the objective function in parallel. Multi-point optimisation is implemented in NUBO through Monte Carlo acquisition functions. The script below uses the MCUpperConfidenceBound acquisition function with 512 samples and fixed base samples. Each batch of 4 is found jointly with the multi_joint() function by optimising the acquisition function with the deterministic
L-BFGS-B optimiser. The Hartmann6D synthetic test function acts as a surrogate for a black-box objective function, such as an experiment or a simulation. The optimisation loop is run for 10 iterations returning batches of 4 each (a total of 40 evaluations) and finds a solution close to the true optimum of -3.3224.
[3]:
import torch
from nubo.acquisition import MCExpectedImprovement, MCUpperConfidenceBound
from nubo.models import GaussianProcess, fit_gp
from nubo.optimisation import multi_joint, multi_sequential
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
bounds = torch.tensor([[0., 0., 0., 0., 0., 0.], [1., 1., 1., 1., 1., 1.]])
# training data
x_train = gen_inputs(num_points=dims*5,
num_dims=dims,
bounds=bounds)
y_train = func(x_train)
# Bayesian optimisation loop
iters = 10
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 = MCExpectedImprovement(gp=gp, y_best=torch.max(y_train), samples=512, fix_base_samples=True)
acq = MCUpperConfidenceBound(gp=gp, beta=1.96**2, samples=512, fix_base_samples=True)
# optimise acquisition function
x_new, _ = multi_joint(func=acq, method="L-BFGS-B", batch_size=4, bounds=bounds, num_starts=5)
# x_new, _ = multi_sequential(func=acq, method="L-BFGS-B", batch_size=4, bounds=bounds, 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 torch.max(y_new) > torch.max(y_train[:-y_new.size(0)]):
best_eval = torch.argmax(y_train)
print(f"New best at evaluation {best_eval+1}: \t Inputs: {x_train[best_eval, :].numpy().reshape(dims).round(4)}, \t Outputs: {-y_train[best_eval].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 36: Inputs: [0. 0. 0. 0.1476 0.276 0.6154], Outputs: -1.5197
New best at evaluation 45: Inputs: [0. 0.1924 0.1315 0.2455 0.2446 0.6663], Outputs: -2.1614
New best at evaluation 48: Inputs: [0. 0.1833 0.6835 0.2924 0.2919 0.6498], Outputs: -2.4447
New best at evaluation 55: Inputs: [0.3088 0.1569 0.6135 0.2713 0.2833 0.6521], Outputs: -2.9539
New best at evaluation 65: Inputs: [0.3024 0.1871 0.46 0.2855 0.304 0.6988], Outputs: -3.1494
New best at evaluation 70: Inputs: [0.218 0.1662 0.4164 0.2816 0.3224 0.6443], Outputs: -3.2737
Evaluation: 70 Solution: -3.2737
Constrained optimisation#
In the second example, NUBO is used to maximise a function where the input space is bounded and constrained. The whole process is not too different from the unconstrained case. We only need to choose a different optimiser that allows the use of constraints when maximising the acquisition function MCUpperConfidenceBound with fixed base samples. At the moment parallel constrained optimisation is only supported by the sequential greedy optimisation strategy using multi_sequential(). NUBO
uses the SLSQP optimiser that can be provided with a dictionary or a tuple of dictionaries that specify one or multiple constraints. We specify two constraints to showcase the two different options: equality constraints and inequality constraints. Equality constraints require the constraint to be 0 while the result is non-negative for inequality constraints. Our first constraint {'type': 'ineq', 'fun': lambda x: 0.5 - x[0] - x[1]} is an inequality constraint and requires the sum of the first
two inputs to be smaller or equal to 0.5. The second constraint {'type': 'eq', 'fun': lambda x: 1.2442 - x[3] - x[4] - x[5]} is an equality constraint specifying that the sum of the last three inputs needs to be equal to 1.2442. The Hartmann6D synthetic test function acts as a substitute for a black-box objective funtion, such as an experiment or a simulation. The optimisation loop is run for 40 iterations and finds a solution close the true optimum of -3.3224. Important: Generating
initial input points with a Latin hypercube might not work for real problems as they will not consider the constraints but only the bounds. In these situations, other methods or selecting initial points by hand might be preferable. The purpose of this example is solely the demonstration of how NUBO handles constraints and constrained optimisation.
[4]:
import torch
from nubo.acquisition import MCExpectedImprovement, MCUpperConfidenceBound
from nubo.models import GaussianProcess, fit_gp
from nubo.optimisation import multi_sequential
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 constraints
bounds = torch.tensor([[0., 0., 0., 0., 0., 0.], [1., 1., 1., 1., 1., 1.]])
cons = ({'type': 'ineq', 'fun': lambda x: 0.5 - x[0] - x[1]},
{'type': 'eq', 'fun': lambda x: 1.2442 - x[3] - x[4] - x[5]})
# training data
x_train = gen_inputs(num_points=dims*5,
num_dims=dims,
bounds=bounds)
y_train = func(x_train)
# Bayesian optimisation loop
iters = 10
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 = MCExpectedImprovement(gp=gp, y_best=torch.max(y_train), samples=512, fix_base_samples=True)
acq = MCUpperConfidenceBound(gp=gp, beta=1.96**2, samples=512, fix_base_samples=True)
# optimise acquisition function
x_new, _ = multi_sequential(func=acq, method="SLSQP", batch_size=4, bounds=bounds, constraints=cons, 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 torch.max(y_new) > torch.max(y_train[:-y_new.size(0)]):
best_eval = torch.argmax(y_train)
print(f"New best at evaluation {best_eval+1}: \t Inputs: {x_train[best_eval, :].numpy().reshape(dims).round(4)}, \t Outputs: {-y_train[best_eval].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 31: Inputs: [0.2252 0. 0.5246 0.3175 0.2892 0.6376], Outputs: -2.9618
New best at evaluation 47: Inputs: [0.3662 0.1338 0.4741 0.312 0.2949 0.6373], Outputs: -2.9727
New best at evaluation 53: Inputs: [0.259 0.241 0.5032 0.2755 0.3047 0.664 ], Outputs: -3.1813
New best at evaluation 55: Inputs: [0.2064 0.1458 0.4293 0.2751 0.3058 0.6633], Outputs: -3.3008
New best at evaluation 59: Inputs: [0.206 0.1526 0.4764 0.2806 0.3169 0.6467], Outputs: -3.3169
New best at evaluation 63: Inputs: [0.21 0.1482 0.486 0.2736 0.3104 0.6602], Outputs: -3.3201
Evaluation: 63 Solution: -3.3201