Noisy Bayesian optimization

Basic Bayesian Optimization¶

In this tutorial we demonstrate the use of Xopt to preform Bayesian Optimization on a simple test problem.

Define the test problem¶

Here we define a simple optimization problem, where we attempt to minimize the sin function in the domian [0,2*pi]. Note that the function used to evaluate the objective function takes a dictionary as input and returns a dictionary as the output.

In [1]:

import torch

from xopt.evaluator import Evaluator
from xopt.generators.bayesian import UpperConfidenceBoundGenerator
from xopt.generators.bayesian.models.standard import StandardModelConstructor
from xopt import Xopt

from xopt.vocs import VOCS

import numpy as np
import math

# define variables and function objectives
vocs = VOCS(
    variables={"x": [0, 2 * math.pi]},
    objectives={"f": "MINIMIZE"},
)

import torch from xopt.evaluator import Evaluator from xopt.generators.bayesian import UpperConfidenceBoundGenerator from xopt.generators.bayesian.models.standard import StandardModelConstructor from xopt import Xopt from xopt.vocs import VOCS import numpy as np import math # define variables and function objectives vocs = VOCS( variables={"x": [0, 2 * math.pi]}, objectives={"f": "MINIMIZE"}, )

In [2]:

# define a test function to optimize
def sin_function(input_dict):
    # return multiple noisy measurements
    return {
        "f": np.sin(input_dict["x"]) + 0.5 * np.random.randn(5),
        "dummy": np.random.randn(5),
    }

# define a test function to optimize def sin_function(input_dict): # return multiple noisy measurements return { "f": np.sin(input_dict["x"]) + 0.5 * np.random.randn(5), "dummy": np.random.randn(5), }

Create Xopt objects¶

Create the evaluator to evaluate our test function and create a generator that uses the Upper Confidence Bound acquisition function to perform Bayesian Optimization.

In [3]:

evaluator = Evaluator(function=sin_function)
gp_constructor = StandardModelConstructor(use_low_noise_prior=False)
generator = UpperConfidenceBoundGenerator(vocs=vocs, gp_constructor=gp_constructor)
X = Xopt(evaluator=evaluator, generator=generator)

evaluator = Evaluator(function=sin_function) gp_constructor = StandardModelConstructor(use_low_noise_prior=False) generator = UpperConfidenceBoundGenerator(vocs=vocs, gp_constructor=gp_constructor) X = Xopt(evaluator=evaluator, generator=generator)

Generate and evaluate initial points¶

To begin optimization, we must generate some random initial data points. The first call to X.step() will generate and evaluate a number of randomly points specified by the generator. Note that if we add data to xopt before calling X.step() by assigning the data to X.data, calls to X.step() will ignore the random generation and proceed to generating points via Bayesian optimization.

In [4]:

# call X.random_evaluate() to generate + evaluate 3 initial points
X.random_evaluate(2)

# inspect the gathered data
X.data

# call X.random_evaluate() to generate + evaluate 3 initial points X.random_evaluate(2) # inspect the gathered data X.data

Out[4]:

	x	f	dummy	xopt_runtime	xopt_error
0	2.009904	0.096715	-0.099704	0.000039	False
1	2.009904	1.160537	-1.793205	0.000039	False
2	2.009904	0.742171	1.155227	0.000039	False
3	2.009904	0.372807	-0.010179	0.000039	False
4	2.009904	0.862657	-1.127788	0.000039	False
5	3.771894	0.065682	-0.891836	0.000011	False
6	3.771894	-0.573799	-0.395535	0.000011	False
7	3.771894	-0.81518	-0.29085	0.000011	False
8	3.771894	-0.751859	-0.755548	0.000011	False
9	3.771894	-0.313088	0.655376	0.000011	False

Do bayesian optimization steps¶

To perform optimization we simply call X.step() in a loop. This allows us to do intermediate tasks in between optimization steps, such as examining the model and acquisition function at each step (as we demonstrate here).

In [5]:

n_steps = 5

# test points for plotting
test_x = torch.linspace(*torch.tensor(X.vocs.bounds).flatten(), 50).double()

for i in range(n_steps):
    # get the Gaussian process model from the generator
    model = X.generator.train_model()

    # visualize model
    fig, ax = X.generator.visualize_model(n_grid=len(test_x))

    # plot true function
    true_f = torch.sin(test_x)
    ax[0, 0].plot(test_x, true_f, "C1--")

    # do the optimization step
    X.step()

n_steps = 5 # test points for plotting test_x = torch.linspace(*torch.tensor(X.vocs.bounds).flatten(), 50).double() for i in range(n_steps): # get the Gaussian process model from the generator model = X.generator.train_model() # visualize model fig, ax = X.generator.visualize_model(n_grid=len(test_x)) # plot true function true_f = torch.sin(test_x) ax[0, 0].plot(test_x, true_f, "C1--") # do the optimization step X.step()

In [6]:

ax[0]

ax[0]

Out[6]:

array([<Axes: ylabel='f'>], dtype=object)

In [7]:

# access the collected data
X.data

# access the collected data X.data

Out[7]:

	x	f	dummy	xopt_runtime	xopt_error
0	2.009904	0.096715	-0.099704	0.000039	False
1	2.009904	1.160537	-1.793205	0.000039	False
2	2.009904	0.742171	1.155227	0.000039	False
3	2.009904	0.372807	-0.010179	0.000039	False
4	2.009904	0.862657	-1.127788	0.000039	False
5	3.771894	0.065682	-0.891836	0.000011	False
6	3.771894	-0.573799	-0.395535	0.000011	False
7	3.771894	-0.81518	-0.29085	0.000011	False
8	3.771894	-0.751859	-0.755548	0.000011	False
9	3.771894	-0.313088	0.655376	0.000011	False
10	4.615591	-1.183509	0.61423	0.000041	False
11	4.615591	-1.130828	-1.552392	0.000041	False
12	4.615591	-0.445194	0.085567	0.000041	False
13	4.615591	-0.82815	2.157936	0.000041	False
14	4.615591	-0.588861	0.934762	0.000041	False
15	6.232973	-0.458577	0.557369	0.000041	False
16	6.232973	-0.568781	-0.488964	0.000041	False
17	6.232973	-0.079223	1.588527	0.000041	False
18	6.232973	1.487086	0.779645	0.000041	False
19	6.232973	-1.161111	0.814123	0.000041	False
20	4.825666	-0.335034	-0.729702	0.000040	False
21	4.825666	-1.311131	-2.038999	0.000040	False
22	4.825666	-1.394021	0.358916	0.000040	False
23	4.825666	-1.139042	-0.111842	0.000040	False
24	4.825666	0.065738	0.470048	0.000040	False
25	4.803151	-1.053018	1.340062	0.000039	False
26	4.803151	-1.371725	-0.733488	0.000039	False
27	4.803151	-0.684791	-1.942095	0.000039	False
28	4.803151	-1.936039	0.623396	0.000039	False
29	4.803151	-0.930968	0.423361	0.000039	False
30	4.871621	-0.500851	0.430609	0.000041	False
31	4.871621	-1.302392	0.333442	0.000041	False
32	4.871621	-0.594924	2.112162	0.000041	False
33	4.871621	-1.556439	0.339615	0.000041	False
34	4.871621	-0.980312	0.773768	0.000041	False

Getting the optimization result¶

To get the best point (without evaluating it) we ask the generator to predict the optimum based on the posterior mean.

In [8]:

X.generator.get_optimum()

X.generator.get_optimum()

Out[8]:

	x
0	4.758693

Customizing optimization¶

Each generator has a set of options that can be modified to effect optimization behavior

In [9]:

X.generator.dict()

X.generator.dict()

/tmp/ipykernel_3436/1542263183.py:1: PydanticDeprecatedSince20: The `dict` method is deprecated; use `model_dump` instead. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.13/migration/
  X.generator.dict()

Out[9]:

{'returns_id': False,
 'supports_single_objective': True,
 'supports_constraints': True,
 'vocs': {'variables': {'x': {'dtype': None,
    'default_value': None,
    'domain': [0.0, 6.283185307179586],
    'type': 'ContinuousVariable'}},
  'objectives': {'f': {'dtype': None, 'type': 'MinimizeObjective'}},
  'constraints': {},
  'constants': {},
  'observables': {}},
 'model': ModelListGP(
   (models): ModuleList(
     (0): SingleTaskGP(
       (likelihood): GaussianLikelihood(
         (noise_covar): HomoskedasticNoise(
           (noise_prior): GammaPrior()
           (raw_noise_constraint): GreaterThan(1.000E-04)
         )
       )
       (mean_module): ConstantMean()
       (covar_module): RBFKernel(
         (lengthscale_prior): LogNormalPrior()
         (raw_lengthscale_constraint): GreaterThan(2.500E-02)
       )
       (outcome_transform): Standardize()
       (input_transform): Normalize()
     )
   )
   (likelihood): LikelihoodList(
     (likelihoods): ModuleList(
       (0): GaussianLikelihood(
         (noise_covar): HomoskedasticNoise(
           (noise_prior): GammaPrior()
           (raw_noise_constraint): GreaterThan(1.000E-04)
         )
       )
     )
   )
 ),
 'n_monte_carlo_samples': 128,
 'turbo_controller': None,
 'use_cuda': False,
 'gp_constructor': {'name': 'standard',
  'use_low_noise_prior': False,
  'covar_modules': {},
  'mean_modules': {},
  'trainable_mean_keys': [],
  'transform_inputs': True,
  'custom_noise_prior': None,
  'use_cached_hyperparameters': False,
  'train_method': 'lbfgs',
  'train_model': True,
  'train_config': None,
  'train_kwargs': None},
 'numerical_optimizer': {'name': 'LBFGS',
  'n_restarts': 20,
  'max_iter': 2000,
  'max_time': 5.0},
 'max_travel_distances': None,
 'fixed_features': None,
 'computation_time':    training  acquisition_optimization
 0  0.709882                  1.352164
 1  1.122954                  0.475173
 2  1.493260                  0.368842
 3  0.678186                  0.407492
 4  1.480506                  0.332324,
 'custom_objective': None,
 'n_interpolate_points': None,
 'n_candidates': 1,
 'beta': 2.0,
 'shift': 0.0}