14  Optimal Computational Budget Allocation in Spot

This chapter demonstrates how noisy functions can be handled with Optimal Computational Budget Allocation (OCBA) by Spot.

14.1 Example: Spot, OCBA, and the Noisy Sphere Function

import numpy as np
from math import inf
from spotPython.fun.objectivefunctions import analytical
from spotPython.spot import spot
import matplotlib.pyplot as plt
from spotPython.utils.init import fun_control_init, get_spot_tensorboard_path
from spotPython.utils.init import fun_control_init, design_control_init, surrogate_control_init

PREFIX = "09"

14.1.1 The Objective Function: Noisy Sphere

The spotPython package provides several classes of objective functions. We will use an analytical objective function with noise, i.e., a function that can be described by a (closed) formula: \[f(x) = x^2 + \epsilon\]

Since sigma is set to 0.1, noise is added to the function:

fun = analytical().fun_sphere
fun_control = fun_control_init(
    PREFIX=PREFIX,
    sigma=0.1)
Created spot_tensorboard_path: runs/spot_logs/09_maans14_2024-04-22_00-27-42 for SummaryWriter()

A plot illustrates the noise:

x = np.linspace(-1,1,100).reshape(-1,1)
y = fun(x, fun_control=fun_control)
plt.figure()
plt.plot(x,y, "k")
plt.show()

Spot is adopted as follows to cope with noisy functions:

  1. fun_repeats is set to a value larger than 1 (here: 2)
  2. noise is set to true. Therefore, a nugget (Lambda) term is added to the correlation matrix
  3. init size (of the design_control dictionary) is set to a value larger than 1 (here: 2)
spot_1_noisy = spot.Spot(fun=fun,
                   fun_control=fun_control_init( 
                   lower = np.array([-1]),
                   upper = np.array([1]),
                   fun_evals = 20,
                   fun_repeats = 2,
                   infill_criterion="ei",
                   noise = True,
                   tolerance_x=0.0,
                   ocba_delta = 1,                   
                   show_models=True),
                   design_control=design_control_init(init_size=3, repeats=2),
                   surrogate_control=surrogate_control_init(noise=True))
spot_1_noisy.run()

spotPython tuning: 0.03475493366922229 [####------] 40.00% 
spotPython tuning: 0.0004463018568303854 [#####-----] 50.00% 
spotPython tuning: 0.0004463018568303854 [######----] 60.00% 
spotPython tuning: 0.0001590474610240226 [#######---] 70.00% 
spotPython tuning: 4.2454542934289965e-09 [########--] 80.00% 
spotPython tuning: 2.2370853591440457e-10 [#########-] 90.00% 
spotPython tuning: 2.2370853591440457e-10 [##########] 100.00% Done...

{'CHECKPOINT_PATH': 'runs/saved_models/',
 'DATASET_PATH': 'data/',
 'PREFIX': None,
 'RESULTS_PATH': 'results/',
 'TENSORBOARD_PATH': 'runs/',
 '_L_in': None,
 '_L_out': None,
 '_torchmetric': None,
 'accelerator': 'auto',
 'converters': None,
 'core_model': None,
 'core_model_name': None,
 'counter': 20,
 'data': None,
 'data_dir': './data',
 'data_module': None,
 'data_set': None,
 'data_set_name': None,
 'db_dict_name': None,
 'design': None,
 'device': None,
 'devices': 1,
 'enable_progress_bar': False,
 'eval': None,
 'fun_evals': 20,
 'fun_repeats': 2,
 'horizon': None,
 'infill_criterion': 'ei',
 'k_folds': 3,
 'log_graph': False,
 'log_level': 50,
 'loss_function': None,
 'lower': array([-1]),
 'max_surrogate_points': 30,
 'max_time': 1,
 'metric_params': {},
 'metric_river': None,
 'metric_sklearn': None,
 'metric_sklearn_name': None,
 'metric_torch': None,
 'model_dict': {},
 'n_points': 1,
 'n_samples': None,
 'n_total': None,
 'noise': True,
 'num_workers': 0,
 'ocba_delta': 1,
 'oml_grace_period': None,
 'optimizer': None,
 'path': None,
 'prep_model': None,
 'prep_model_name': None,
 'progress_file': None,
 'save_model': False,
 'scenario': None,
 'seed': 123,
 'show_batch_interval': 1000000,
 'show_models': True,
 'show_progress': True,
 'shuffle': None,
 'sigma': 0.0,
 'spot_tensorboard_path': None,
 'spot_writer': None,
 'target_column': None,
 'target_type': None,
 'task': None,
 'test': None,
 'test_seed': 1234,
 'test_size': 0.4,
 'tolerance_x': 0.0,
 'train': None,
 'upper': array([1]),
 'var_name': None,
 'var_type': ['num'],
 'verbosity': 0,
 'weight_coeff': 0.0,
 'weights': 1.0,
 'weights_entry': None}

14.3 Noise and Surrogates: The Nugget Effect

14.3.1 The Noisy Sphere

14.3.1.1 The Data

We prepare some data first:

import numpy as np
import spotPython
from spotPython.fun.objectivefunctions import analytical
from spotPython.spot import spot
from spotPython.design.spacefilling import spacefilling
from spotPython.build.kriging import Kriging
import matplotlib.pyplot as plt

gen = spacefilling(1)
rng = np.random.RandomState(1)
lower = np.array([-10])
upper = np.array([10])
fun = analytical().fun_sphere
fun_control = fun_control_init(    
    sigma=2,
    seed=125)
X = gen.scipy_lhd(10, lower=lower, upper = upper)
y = fun(X, fun_control=fun_control)
X_train = X.reshape(-1,1)
y_train = y

A surrogate without nugget is fitted to these data:

S = Kriging(name='kriging',
            seed=123,
            log_level=50,
            n_theta=1,
            noise=False)
S.fit(X_train, y_train)

X_axis = np.linspace(start=-13, stop=13, num=1000).reshape(-1, 1)
mean_prediction, std_prediction, ei = S.predict(X_axis, return_val="all")

plt.scatter(X_train, y_train, label="Observations")
plt.plot(X_axis, mean_prediction, label="mue")
plt.legend()
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
_ = plt.title("Sphere: Gaussian process regression on noisy dataset")

In comparison to the surrogate without nugget, we fit a surrogate with nugget to the data:

S_nug = Kriging(name='kriging',
            seed=123,
            log_level=50,
            n_theta=1,
            noise=True)
S_nug.fit(X_train, y_train)
X_axis = np.linspace(start=-13, stop=13, num=1000).reshape(-1, 1)
mean_prediction, std_prediction, ei = S_nug.predict(X_axis, return_val="all")
plt.scatter(X_train, y_train, label="Observations")
plt.plot(X_axis, mean_prediction, label="mue")
plt.legend()
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
_ = plt.title("Sphere: Gaussian process regression with nugget on noisy dataset")

The value of the nugget term can be extracted from the model as follows:

S.Lambda
S_nug.Lambda
8.374496269458742e-05

We see:

  • the first model S has no nugget,
  • whereas the second model has a nugget value (Lambda) larger than zero.

14.4 Exercises

14.4.1 Noisy fun_cubed

Analyse the effect of noise on the fun_cubed function with the following settings:

fun = analytical().fun_cubed
fun_control = fun_control_init(    
    sigma=10,
    seed=123)
lower = np.array([-10])
upper = np.array([10])

14.4.2 fun_runge

Analyse the effect of noise on the fun_runge function with the following settings:

lower = np.array([-10])
upper = np.array([10])
fun = analytical().fun_runge
fun_control = fun_control_init(    
    sigma=0.25,
    seed=123)

14.4.3 fun_forrester

Analyse the effect of noise on the fun_forrester function with the following settings:

lower = np.array([0])
upper = np.array([1])
fun = analytical().fun_forrester
fun_control = {"sigma": 5,
               "seed": 123}

14.4.4 fun_xsin

Analyse the effect of noise on the fun_xsin function with the following settings:

lower = np.array([-1.])
upper = np.array([1.])
fun = analytical().fun_xsin
fun_control = fun_control_init(    
    sigma=0.5,
    seed=123)