Backtesting·Strategy Optimization·Intermediate

Bayesian Optimization

Implement Bayesian optimization for strategy parameter tuning using Gaussian process surrogate models to efficiently search high-dimensional parameter spaces, intelligently balancing exploration of uncertain regions with exploitation of known high-performing parameter zones.

backtestingoptimization

Bayesian Optimization: A Comprehensive Guide

Bayesian Optimization is a powerful global optimization strategy for objective functions that are expensive to evaluate, noisy, or lack a known analytical form (e.g., black-box functions). It's particularly useful in machine learning for hyperparameter tuning, where evaluating a single set of hyperparameters can take hours or days.

Why Bayesian Optimization?

Traditional optimization methods often struggle with black-box functions. For instance:

Gradient-based methods require gradients, which might not be available.
Random search can be inefficient as it doesn't learn from past evaluations.
Grid search becomes computationally infeasible quickly with increasing dimensions.

Bayesian Optimization, in contrast, builds a probabilistic model of the objective function (the surrogate model) based on previous evaluations. It then uses this model to intelligently decide where to sample next, aiming to find the global optimum efficiently.

Core Components

Bayesian Optimization primarily consists of two main components:

Surrogate Model (Probabilistic Model): A statistical model that approximates the expensive objective function. It provides both an estimate of the function's value and the uncertainty around that estimate.
Acquisition Function: A criterion that uses the surrogate model's predictions to determine the next most promising point to evaluate. It balances exploration (sampling where uncertainty is high) and exploitation (sampling where the model predicts a high objective value).

1. Surrogate Model: Gaussian Processes

A common choice for the surrogate model is a Gaussian Process (GP). A GP is a collection of random variables, any finite number of which have a joint Gaussian distribution. It defines a distribution over functions, allowing us to not only predict the mean value of the objective function at any point but also the uncertainty (variance) of that prediction.

Key features of GPs for Bayesian Optimization:

Probabilistic predictions: Provides a mean and variance for each prediction.
Uncertainty quantification: The variance naturally quantifies how confident the model is in its predictions, which is crucial for balancing exploration and exploitation.
Non-parametric: Can model complex functions without assuming a specific functional form.

[2]

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern

def objective_function(x):
    """
    Evaluates the objective function at point x.

    Inputs:
    x (np.array): A 1D array of points where the function is evaluated.

    Outputs:
    np.array: The function values at the given points.
    """
    return -np.sin(x) * (x - 2)**2 + x/2

# Plot the true objective function for reference
x_true = np.linspace(0, 10, 100).reshape(-1, 1)
y_true = objective_function(x_true)

plt.figure(figsize=(10, 6))
plt.plot(x_true, y_true, label='True Objective Function', color='blue')
plt.title('True Objective Function to be Optimized')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
plt.show()

2. Acquisition Function: Expected Improvement (EI)

Acquisition functions guide the search by proposing the next point to evaluate. They are designed to balance two competing goals:

Exploitation: Sampling where the surrogate model predicts a high objective value (to improve the current best).
Exploration: Sampling where the surrogate model is uncertain (to reduce uncertainty and potentially discover new optima).

Expected Improvement (EI) is a popular acquisition function. It quantifies the expected amount of improvement over the current best observed objective value, considering both the mean prediction and the uncertainty of the surrogate model. It's calculated as:

$$ EI(x) = E[max(0, f(x) - f_{best})] $$

Where:

$f(x)$ is the predicted value of the objective function at point $x$.
$f_{best}$ is the current best observed objective value.
$E[\cdot]$ denotes the expectation, which for a Gaussian Process can be computed analytically.

[3]

def expected_improvement(X, gp_model, current_best_f):
    """
    Calculates the Expected Improvement (EI) for a given set of points X (for minimization).
    """
    mu, sigma = gp_model.predict(X, return_std=True)
    mu = mu.flatten()
    sigma = sigma.flatten()

    # For minimization, we want f_best - f(x)
    # We define Z as (current_best_f - mu) / sigma
    # EI = (current_best_f - mu) * norm.cdf(Z) + sigma * norm.pdf(Z)

    with np.errstate(divide='warn'):
        Z = np.zeros_like(sigma)
        non_zero_sigma_indices = sigma > 1e-10
        Z[non_zero_sigma_indices] = (current_best_f - mu[non_zero_sigma_indices]) / sigma[non_zero_sigma_indices]

        ei = (current_best_f - mu) * norm.cdf(Z) + sigma * norm.pdf(Z)

    return np.maximum(0, ei)

Bayesian Optimization Algorithm Workflow

The Bayesian Optimization process iteratively refines its understanding of the objective function:

Initialization: Evaluate the objective function at a small number of randomly chosen points.
Model Training: Fit a Gaussian Process (surrogate model) to all observed data points (inputs and their corresponding objective values).
Acquisition Function Optimization: Use the fitted GP to calculate the acquisition function (e.g., Expected Improvement) across a dense set of candidate points. Find the point that maximizes the acquisition function.
Objective Evaluation: Evaluate the true objective function at the new point identified in step 3.
Update Data: Add the new point and its objective value to the set of observed data.
Repeat: Go back to step 2 and repeat until a stopping criterion is met (e.g., maximum number of iterations, convergence).

[4]

def bayesian_optimization(objective_func, bounds, n_iter, n_initial_points=3, random_seed=42):
    """
    Performs Bayesian Optimization to find the minimum of an objective function.
    """
    np.random.seed(random_seed)

    # 1. Initialization: Sample initial random points
    n_dims = len(bounds)
    X_initial = np.random.uniform([b[0] for b in bounds], [b[1] for b in bounds],
                                   size=(n_initial_points, n_dims))
    y_initial = np.array([objective_func(x) for x in X_initial]).reshape(-1, 1)

    X_observed = X_initial
    y_observed = y_initial

    # Keep track of the best observed value at each step
    best_f_history = [np.min(y_observed)]

    # Define the Gaussian Process Regressor
    kernel = Matern(length_scale=1.0, nu=2.5)
    gp_model = GaussianProcessRegressor(kernel=kernel, alpha=1e-6, normalize_y=True,
                                        n_restarts_optimizer=10, random_state=random_seed)

    for i in range(n_iter):
        # 2. Model Training: Fit GP to observed data
        gp_model.fit(X_observed, y_observed)

        # Current best observed function value
        current_best_f = np.min(y_observed)

        # 3. Acquisition Function Optimization
        # Generate candidate points - ensure 2D shape
        x_candidates = np.linspace(bounds[0][0], bounds[0][1], 1000).reshape(-1, 1)

        # Calculate EI for all candidates
        ei_values = expected_improvement(x_candidates, gp_model, current_best_f)

        # Find the point with the maximum EI
        next_x_idx = np.argmax(ei_values)
        next_x = x_candidates[next_x_idx].reshape(1, -1)  # Fixed: reshape to 2D

        # 4. Objective Evaluation
        next_f = objective_func(next_x).reshape(1, -1)  # Fixed: ensure 2D shape

        # 5. Update Data - ensure both are 2D
        X_observed = np.vstack((X_observed, next_x))
        y_observed = np.vstack((y_observed, next_f))

        best_f_history.append(np.min(y_observed))

    best_idx = np.argmin(y_observed)
    best_x = X_observed[best_idx]
    best_f = y_observed[best_idx]

    return best_x, best_f, list(zip(X_observed, y_observed)), best_f_history

Demonstration and Visualization

Let's apply Bayesian Optimization to our objective_function and visualize the process over several iterations. We will see how the Gaussian Process updates its belief about the function and how the acquisition function guides the search.

[5]

# Define the search space bounds
bounds = [(0, 10)]

# Run Bayesian Optimization
n_iterations = 15
n_initial = 3
best_x, best_f, all_observations, best_f_history = bayesian_optimization(
    objective_function, bounds, n_iterations, n_initial_points=n_initial, random_seed=42
)

print(f"Optimal x found: {best_x.item():.4f}")
print(f"Minimum f(x) found: {best_f.item():.4f}")

# Extract observed points for plotting
X_obs = np.array([obs[0] for obs in all_observations])
y_obs = np.array([obs[1] for obs in all_observations]).reshape(-1, 1)

# Prepare for plotting the GP and Acquisition function
kernel = Matern(length_scale=1.0, nu=2.5)
gp_final = GaussianProcessRegressor(kernel=kernel, alpha=1e-6, normalize_y=True,
                                    n_restarts_optimizer=10, random_state=42)
gp_final.fit(X_obs, y_obs)

x_plot = np.linspace(bounds[0][0], bounds[0][1], 500).reshape(-1, 1)
mu_plot, sigma_plot = gp_final.predict(x_plot, return_std=True)

current_best_f_final = np.min(y_obs)
ei_plot = expected_improvement(x_plot, gp_final, current_best_f_final)

# Visualization 1: Objective function, GP model, and Acquisition Function
plt.figure(figsize=(12, 8))

# Plot true objective function
plt.plot(x_true, y_true, 'blue', label='True Objective Function')

# Plot GP mean and uncertainty
plt.plot(x_plot, mu_plot, 'r--', label='GP Mean Prediction')
plt.fill_between(x_plot.flatten(), mu_plot - 1.96 * sigma_plot,
                 mu_plot + 1.96 * sigma_plot, alpha=0.2, color='red',
                 label='95% Confidence Interval')

# Plot observed points
plt.scatter(X_obs[:n_initial], y_obs[:n_initial], color='green',
           marker='o', s=100, label='Initial Samples')
plt.scatter(X_obs[n_initial:], y_obs[n_initial:], color='purple',
           marker='x', s=100, label='BO Samples')
plt.scatter(best_x, best_f, color='cyan', marker='*', s=300, label='Best Found Point')

# Plot acquisition function
ax2 = plt.gca().twinx()
ax2.plot(x_plot, ei_plot, color='orange', linestyle=':', label='Expected Improvement')
ax2.set_ylabel('Expected Improvement')

plt.title('Bayesian Optimization Progress: Objective, GP, and EI')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.grid(True)
plt.show()

# Visualization 2: Convergence Plot
plt.figure(figsize=(10, 6))
plt.plot(range(len(best_f_history)), best_f_history, marker='o', linestyle='-', color='red')
plt.title('Convergence of Bayesian Optimization (Best Objective Value)')
plt.xlabel('Iteration')
plt.ylabel('Best Observed f(x)')
plt.grid(True)
plt.show()

/usr/local/lib/python3.12/dist-packages/sklearn/gaussian_process/_gpr.py:660: ConvergenceWarning: lbfgs failed to converge (status=2):
ABNORMAL: .

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
  _check_optimize_result("lbfgs", opt_res)

Optimal x found: 8.1582
Minimum f(x) found: -32.1029

Interpretation of Visualization 1:

The True Objective Function (blue line) is what we aim to optimize. In real-world scenarios, this function is unknown.
The GP Mean Prediction (red dashed line) shows the Gaussian Process's current estimate of the objective function's shape.
The 95% Confidence Interval (red shaded area) represents the uncertainty of the GP's prediction. Wider areas indicate higher uncertainty, which promotes exploration.
Initial Samples (green circles) are the points randomly chosen to kickstart the process.
BO Samples (purple crosses) are the points chosen by the Bayesian Optimization algorithm, guided by the acquisition function.
Best Found Point (cyan star) marks the best objective value discovered so far.
Expected Improvement (orange dotted line) shows where the algorithm believes the most improvement can be gained. The peaks of this function indicate the next most likely points to sample. Notice how the algorithm tends to sample near these peaks, balancing between exploiting known good regions and exploring uncertain ones.

Conclusion

Bayesian Optimization provides an efficient framework for optimizing expensive, black-box functions. By leveraging a probabilistic surrogate model and an acquisition function, it intelligently explores the search space, balancing exploration and exploitation to quickly converge to optimal solutions. It is widely applied in various fields, especially for hyperparameter optimization in machine learning, materials discovery, and experimental design, where function evaluations are costly.