Asynchronous Bayesian optimisation

In this example, NUBO is used for asynchronous optimisation. This means that the optimisation loop is continued while some points are still being evaluated from the objective function. This is particularly useful for situations in which some evaluations take a longer time to complete but you do not want to waste time by waiting for these pending observations. In the script below, we randomly sample a pending point x_pending and assume that we are still waiting for its output. While waiting, we continue the optimisation loop for 10 iterations with a batch size of 4 each (a total of 40 evaluations) and find a solution close to the true optimum of -3.3224. The Hartmann6D synthetic test function acts as a substitute for a black-box objective function, such as an experiment or a simulation. Notice that we provide the pending point x_pending to the acquisition function MCUpperConfidenceBound as an argument. For asynchronous optimisation, Monte Carlo acquisition functions have to be used as this process is in general intractable for analytical functions.

import torch
from nubo.acquisition import MCExpectedImprovement, MCUpperConfidenceBound
from nubo.models import GaussianProcess, fit_gp
from nubo.optimisation import multi_joint
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)

# point pending evaluation
x_pending = torch.rand((1, dims))
print(f"Point pending evaluation: {x_pending.numpy().reshape(dims).round(4)}")

# 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), x_pending=x_pending, samples=512)
    acq = MCUpperConfidenceBound(gp=gp, beta=1.96**2, x_pending=x_pending, samples=512)

    # optimise acquisition function
    x_new, _ = multi_joint(func=acq, method="Adam", batch_size=4, bounds=bounds, lr=0.1, steps=200, 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}")

Point pending evaluation: [0.1874 0.0246 0.4963 0.5481 0.5602 0.5434]
New best at evaluation 35:       Inputs: [0.4185 0.9953 0.9987 0.4276 0.5028 0.0036],    Outputs: -2.2505
New best at evaluation 48:       Inputs: [0.403  0.9975 0.0053 0.4634 0.0037 0.0044],    Outputs: -2.4121
New best at evaluation 53:       Inputs: [4.010e-01 8.998e-01 9.766e-01 4.932e-01 1.000e-04 5.000e-04],          Outputs: -2.926
New best at evaluation 56:       Inputs: [4.135e-01 8.917e-01 9.949e-01 5.736e-01 7.400e-03 2.000e-04],          Outputs: -3.1253
New best at evaluation 61:       Inputs: [0.4098 0.8883 0.9959 0.5746 0.01   0.0463],    Outputs: -3.1903
Evaluation: 61   Solution: -3.1903