Typical problems#

This notebook covers some problems that are commonly encountered in black-box optimisation and how they can be optimised with the off-the-shelf optimise function. This function combines everything required for one optimisation step and returns one or multiple candidate points.

Single-point optimisation#

In this example, NUBO is used for sequential single-point optimisation. The Hartmann6D synthetic test function acts as a substitute for a black-box objective function, such as an experiment or a simulation. The optimise function uses the analytical ExpectedImprovement acquisition function and optimies it via the L-BFGS-B algorithm by default. 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.algorithms import optimise
from nubo.test_functions import Hartmann6D
from nubo.utils import gen_inputs


# 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 = 40

for iter in range(iters):

    # NUBO
    x_new = optimise(x_train, y_train, bounds=bounds)

    # 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 31:       Inputs: [0.4805 0.1888 0.6875 0.1961 0.2543 0.5498],    Outputs: [-1.761]
New best at evaluation 40:       Inputs: [0.363  0.1418 0.6915 0.2731 0.2711 0.6849],    Outputs: [-2.5578]
New best at evaluation 43:       Inputs: [0.2742 0.1642 0.496  0.2955 0.2636 0.7238],    Outputs: [-3.0372]
Evaluation: 43   Solution: -3.0372

Constrained multi-point optimisation#

In this example, 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 64 samples. Each batch of 4 is found sequentially (also known as greedy optimisation) by optimising the acquisition function usually with the stochastic Adam optimiser. However, we also consider two constraints on the input space 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 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.

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.

[2]:
import torch
from nubo.algorithms import optimise
from nubo.test_functions import Hartmann6D
from nubo.utils import gen_inputs


# 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.]])
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):

    # NUBO
    x_new = optimise(x_train, y_train,
                     bounds=bounds,
                     batch_size=4,
                     acquisition="UCB",
                     beta=5.0,
                     constraints=cons,
                     mc_samples=64)

    # 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 35:       Inputs: [0.3274 0.1726 0.3425 0.2882 0.2615 0.6945],    Outputs: -2.934
New best at evaluation 43:       Inputs: [0.2317 0.2549 0.3728 0.2457 0.2965 0.7019],    Outputs: -3.0234
New best at evaluation 48:       Inputs: [0.2224 0.1383 0.4157 0.2715 0.3251 0.6475],    Outputs: -3.2717
New best at evaluation 51:       Inputs: [0.1882 0.1534 0.4292 0.2797 0.3012 0.6633],    Outputs: -3.2934
New best at evaluation 59:       Inputs: [0.2108 0.1393 0.5028 0.2754 0.3058 0.6631],    Outputs: -3.3101
New best at evaluation 63:       Inputs: [0.1893 0.1625 0.4902 0.2713 0.3125 0.6604],    Outputs: -3.3168
Evaluation: 63   Solution: -3.3168

Noisy observations with continuous and discrete parameters#

In this example, NUBO is used for sequential single-point optimisation with continuous and discrete parameters and noisy observations. Additionally to the bounds, a dictionary containing the dimensions as keys and the possible values as values have to be specified for the discrete values. 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 by specifying acquisition="UCB" with a trade-off parameter beta=5.0. The dictionary of discrete values is provided to the optimise function and the noisy argument is set to True to allow he optimisation of a noisy acquisition function. The optimisation loop is run for 40 iterations and finds a solution close to the true optimum of -3.3224.

[3]:
import torch
from nubo.algorithms import optimise
from nubo.test_functions import Hartmann6D
from nubo.utils import gen_inputs


# test function
func = Hartmann6D(minimise=False, noise_std=0.05)
dims = 6

# specify bounds
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):

    # NUBO
    x_new = optimise(x_train, y_train,
                     bounds=bounds,
                     acquisition="UCB",
                     beta=5.0,
                     discrete=discrete,
                     noisy=True)

    # 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 38:       Inputs: [0.4    1.     0.     0.5369 0.3    0.    ],    Outputs: [-2.6125]
New best at evaluation 42:       Inputs: [0.4    0.922  0.     0.5468 0.3    0.    ],    Outputs: [-2.9671]
New best at evaluation 45:       Inputs: [0.4    0.9201 1.     0.5586 0.3    0.    ],    Outputs: [-3.0494]
New best at evaluation 46:       Inputs: [0.4    0.9158 1.     0.5582 0.3    0.0571],    Outputs: [-3.1341]
New best at evaluation 49:       Inputs: [0.4    0.8774 1.     0.561  0.3    0.0419],    Outputs: [-3.1727]
New best at evaluation 51:       Inputs: [0.4    0.8744 1.     0.5736 0.3    0.0454],    Outputs: [-3.1938]
New best at evaluation 52:       Inputs: [0.4    0.8617 1.     0.5805 0.3    0.0539],    Outputs: [-3.2136]
New best at evaluation 57:       Inputs: [0.4    0.8721 1.     0.575  0.3    0.0361],    Outputs: [-3.2376]
New best at evaluation 60:       Inputs: [0.4    0.8715 1.     0.5715 0.3    0.0421],    Outputs: [-3.2734]
Evaluation: 60   Solution: -3.2734