Fork me on GitHub Top

Bayesian optimization with skopt

Gilles Louppe, Manoj Kumar July 2016.

import numpy as np
np.random.seed(123)

%matplotlib inline
import matplotlib.pyplot as plt

Problem statement

We are interested in solving $$x^* = \arg \min_x f(x)$$ under the constraints that

Disclaimer. If you do not have these constraints, then there is certainly a better optimization algorithm than Bayesian optimization.

Bayesian optimization loop

For $t=1:T$:

  1. Given observations $(x_i, y_i=f(x_i))$ for $i=1:t$, build a probabilistic model for the objective $f$. Integrate out all possible true functions, using Gaussian process regression.

  2. optimize a cheap acquisition/utility function $u$ based on the posterior distribution for sampling the next point. $$x_{t+1} = \arg \min_x u(x)$$ Exploit uncertainty to balance exploration against exploitation.

  3. Sample the next observation $y_{t+1}$ at $x_{t+1}$.

Acquisition functions

Acquisition functions $\text{u}(x)$ specify which sample $x$ should be tried next:

where $x_t^+$ is the best point observed so far.

In most cases, acquisition functions provide knobs (e.g., $\kappa$) for controlling the exploration-exploitation trade-off. - Search in regions where $\mu_{GP}(x)$ is high (exploitation) - Probe regions where uncertainty $\sigma_{GP}(x)$ is high (exploration)

Toy example

Let assume the following noisy function $f$:

noise_level = 0.1

def f(x, noise_level=noise_level):
    return np.sin(5 * x[0]) * (1 - np.tanh(x[0] ** 2)) + np.random.randn() * noise_level

Note. In skopt, functions $f$ are assumed to take as input a 1D vector $x$ represented as an array-like and to return a scalar $f(x)$.

# Plot f(x) + contours
x = np.linspace(-2, 2, 400).reshape(-1, 1)
fx = [f(x_i, noise_level=0.0) for x_i in x]
plt.plot(x, fx, "r--", label="True (unknown)")
plt.fill(np.concatenate([x, x[::-1]]),
         np.concatenate(([fx_i - 1.9600 * noise_level for fx_i in fx], 
                         [fx_i + 1.9600 * noise_level for fx_i in fx[::-1]])),
         alpha=.2, fc="r", ec="None")
plt.legend()
plt.grid()
plt.show()

png

Bayesian optimization based on gaussian process regression is implemented in skopt.gp_minimize and can be carried out as follows:

from skopt import gp_minimize

res = gp_minimize(f,                  # the function to minimize
                  [(-2.0, 2.0)],      # the bounds on each dimension of x
                  acq_func="EI",      # the acquisition function
                  n_calls=15,         # the number of evaluations of f 
                  n_random_starts=5,  # the number of random initialization points
                  noise=0.1**2,       # the noise level (optional)
                  random_state=123)   # the random seed

Accordingly, the approximated minimum is found to be:

"x^*=%.4f, f(x^*)=%.4f" % (res.x[0], res.fun)
'x^*=-0.2573, f(x^*)=-1.0715'

For further inspection of the results, attributes of the res named tuple provide the following information:

print(res)
          fun: -1.0714985908064383
    func_vals: array([-0.18788762, -0.96229886, -0.34643484, -0.46587165, -0.1773319 ,
       -0.80809029, -1.07149859, -0.92120939, -0.1608395 , -0.80998886,
       -0.32494081,  0.0329113 , -0.02279617,  0.05849308,  0.69296704])
       models: [GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5) + WhiteKernel(noise_level=0.01),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734)]
 random_state: <mtrand.RandomState object at 0x7f156b9e2bd0>
        space: Space([Real(low=-2.0, high=2.0, prior=uniform, transform=normalize)])
        specs: {'args': {'n_jobs': 1, 'kappa': 1.96, 'xi': 0.01, 'n_restarts_optimizer': 5, 'n_points': 10000, 'callback': None, 'verbose': False, 'random_state': <mtrand.RandomState object at 0x7f156b9e2bd0>, 'y0': None, 'x0': None, 'acq_optimizer': 'auto', 'acq_func': 'EI', 'n_random_starts': 5, 'n_calls': 15, 'base_estimator': GaussianProcessRegressor(alpha=1e-10, copy_X_train=True,
             kernel=1**2 * Matern(length_scale=1, nu=2.5),
             n_restarts_optimizer=2, noise=0.010000000000000002,
             normalize_y=True, optimizer='fmin_l_bfgs_b',
             random_state=843828734), 'dimensions': Space([Real(low=-2.0, high=2.0, prior=uniform, transform=normalize)]), 'func': <function f at 0x7f157857c158>}, 'function': 'base_minimize'}
            x: [-0.2572822321543613]
      x_iters: [[0.85182128209290919], [-0.28611629525269677], [0.76353942010744724], [0.87660124061909261], [-0.035524266269610472], [-0.30626086081540094], [-0.2572822321543613], [-0.22378966562926261], [1.0553495096977938], [-0.25562401679960356], [-0.52442654243976894], [-1.9992867763849413], [1.999866281741757], [-1.2339044799440428], [0.45960222721663557]]

Together these attributes can be used to visually inspect the results of the minimization, such as the convergence trace or the acquisition function at the last iteration:

from skopt.plots import plot_convergence
plot_convergence(res);

png

Let us now visually examine

  1. The approximation of the fit gp model to the original function.
  2. The acquistion values that determine the next point to be queried.
from skopt.acquisition import gaussian_ei

plt.rcParams["figure.figsize"] = (8, 14)

x = np.linspace(-2, 2, 400).reshape(-1, 1)
x_gp = res.space.transform(x.tolist())
fx = np.array([f(x_i, noise_level=0.0) for x_i in x])

# Plot the 5 iterations following the 5 random points
for n_iter in range(5):
    gp = res.models[n_iter]
    curr_x_iters = res.x_iters[:5+n_iter]
    curr_func_vals = res.func_vals[:5+n_iter]

    # Plot true function.
    plt.subplot(5, 2, 2*n_iter+1)
    plt.plot(x, fx, "r--", label="True (unknown)")
    plt.fill(np.concatenate([x, x[::-1]]),
             np.concatenate([fx - 1.9600 * noise_level, 
                             fx[::-1] + 1.9600 * noise_level]),
             alpha=.2, fc="r", ec="None")

    # Plot GP(x) + contours
    y_pred, sigma = gp.predict(x_gp, return_std=True)
    plt.plot(x, y_pred, "g--", label=r"$\mu_{GP}(x)$")
    plt.fill(np.concatenate([x, x[::-1]]),
             np.concatenate([y_pred - 1.9600 * sigma, 
                             (y_pred + 1.9600 * sigma)[::-1]]),
             alpha=.2, fc="g", ec="None")

    # Plot sampled points
    plt.plot(curr_x_iters, curr_func_vals,
             "r.", markersize=8, label="Observations")

    # Adjust plot layout
    plt.grid()

    if n_iter == 0:
        plt.legend(loc="best", prop={'size': 6}, numpoints=1)

    if n_iter != 4:
        plt.tick_params(axis='x', which='both', bottom='off', 
                        top='off', labelbottom='off') 

    # Plot EI(x)
    plt.subplot(5, 2, 2*n_iter+2)
    acq = gaussian_ei(x_gp, gp, y_opt=np.min(curr_func_vals))
    plt.plot(x, acq, "b", label="EI(x)")
    plt.fill_between(x.ravel(), -2.0, acq.ravel(), alpha=0.3, color='blue')

    next_x = res.x_iters[5+n_iter]
    next_acq = gaussian_ei(res.space.transform([next_x]), gp, y_opt=np.min(curr_func_vals))
    plt.plot(next_x, next_acq, "bo", markersize=6, label="Next query point")

    # Adjust plot layout
    plt.ylim(0, 0.1)
    plt.grid()

    if n_iter == 0:
        plt.legend(loc="best", prop={'size': 6}, numpoints=1)

    if n_iter != 4:
        plt.tick_params(axis='x', which='both', bottom='off', 
                        top='off', labelbottom='off') 

plt.show()

png

The first column shows the following:

  1. The true function.
  2. The approximation to the original function by the gaussian process model
  3. How sure the GP is about the function.

The second column shows the acquisition function values after every surrogate model is fit. It is possible that we do not choose the global minimum but a local minimum depending on the minimizer used to minimize the acquisition function.

At the points closer to the points previously evaluated at, the variance dips to zero.

Finally, as we increase the number of points, the GP model approaches the actual function. The final few points are clustered around the minimum because the GP does not gain anything more by further exploration:

plt.rcParams["figure.figsize"] = (6, 4)

# Plot f(x) + contours
x = np.linspace(-2, 2, 400).reshape(-1, 1)
x_gp = res.space.transform(x.tolist())

fx = [f(x_i, noise_level=0.0) for x_i in x]
plt.plot(x, fx, "r--", label="True (unknown)")
plt.fill(np.concatenate([x, x[::-1]]),
         np.concatenate(([fx_i - 1.9600 * noise_level for fx_i in fx], 
                         [fx_i + 1.9600 * noise_level for fx_i in fx[::-1]])),
         alpha=.2, fc="r", ec="None")

# Plot GP(x) + contours
gp = res.models[-1]
y_pred, sigma = gp.predict(x_gp, return_std=True)

plt.plot(x, y_pred, "g--", label=r"$\mu_{GP}(x)$")
plt.fill(np.concatenate([x, x[::-1]]),
         np.concatenate([y_pred - 1.9600 * sigma, 
                         (y_pred + 1.9600 * sigma)[::-1]]),
         alpha=.2, fc="g", ec="None")

# Plot sampled points
plt.plot(res.x_iters, 
         res.func_vals, 
         "r.", markersize=15, label="Observations")

plt.title(r"$x^* = %.4f, f(x^*) = %.4f$" % (res.x[0], res.fun))
plt.legend(loc="best", prop={'size': 8}, numpoints=1)
plt.grid()

plt.show()

png