Bayesian Optimization

Summary

Bayesian Optimization is a powerful and mathematically elegant approach for optimizing expensive-to-evaluate functions, particularly in the domain of hyperparameter tuning in machine learning. Unlike classical methods such as grid search or random search, Bayesian Optimization builds a probabilistic model of the objective function and uses it strategically to decide where to evaluate the function next. This approach finds high-quality solutions with a dramatically reduced number of evaluations, which is a decisive advantage when each evaluation requires hours of model training.

Mathematical Principle

The mathematical foundation of Bayesian Optimization rests on two essential concepts: modeling the objective function and the strategy for acquiring the next point.

Gaussian Process

The first step consists of modeling the unknown objective function f(x) as a Gaussian process. Mathematically, we assume that the posterior distribution of f given the observations X follows a multivariate normal distribution:

p(f|X) = N(μ, K)

where μ represents the prior mean vector and K is the covariance matrix computed using a kernel, typically the RBF (Radial Basis Function) kernel. The RBF kernel is expressed as follows:

k(x_i, x_j) = σ² · exp(-||x_i – x_j||² / (2 · l²))

where σ² is the variance and l is the length scale that controls the smoothness of the function. This Gaussian process assumption means that any finite linear combination of values of f follows a Gaussian distribution. This allows the posterior distribution to be computed analytically after each new observation.

When we have n observations (x₁, y₁), …, (xₙ, yₙ), the prediction at a new point x* is obtained by the following formulas:

μ(x) = kᵀ · (K + σₙ²I)⁻¹ · y

σ²(x) = k(x, x) – kᵀ · (K + σₙ²I)⁻¹ · k*

where k is the vector of covariances between x and the training points, K is the covariance matrix computed on the existing observations, and σₙ² represents the measurement noise.

Acquisition Function: Expected Improvement

Once the Gaussian process is established, we need to decide where to evaluate the function next. This is the role of the acquisition function. The most popular is Expected Improvement, denoted EI(x):

EI(x) = (μ(x) – f(x_best) – ξ)·Φ(Z) + σ(x)·φ(Z)

where:

Z = (μ(x) – f(x_best) – ξ) / σ(x)

In this formula, μ(x) is the mean predicted by the Gaussian process at point x, f(x_best) is the best value observed so far, σ(x) is the predicted standard deviation, ξ (xi) is an exploration parameter that controls the exploration-exploitation tradeoff, Φ denotes the cumulative distribution function of the standard normal distribution and φ is the probability density function of that same distribution.

We choose the point x that maximizes EI(x) for the next evaluation. This point represents the best tradeoff between exploring uncertain regions (high σ(x)) and exploiting promising regions (high μ(x)).

Other acquisition functions also exist:

• Upper Confidence Bound (UCB): μ(x) + κ · σ(x), where κ controls the exploration level
• Probability of Improvement (PI): Φ(Z), simpler but less effective in practice
• Expected Quantile Improvement: a robust variant in the presence of noise

Intuition

Bayesian Optimization is like a smart gold prospector. Imagine you need to find the richest gold vein in a vast territory, but each shovel strike costs dearly in time and effort.

With grid search, you dig systematically at regular intervals, like a grid pattern. You cover everything, but you spend an enormous amount of time in barren areas. It’s methodical but extremely inefficient.

With random search, you dig anywhere, hoping to stumble upon a good area. It’s better than grid search in high-dimensional spaces, but you learn nothing from your previous shovel strikes.

Bayesian Optimization, on the other hand, acts like an experienced prospector. After each shovel strike, it builds a map of promising areas based on everything it has learned so far:

• The Gaussian process is the map: from the few points evaluated, it estimates not only the probable value of gold in each unexplored location, but also the associated uncertainty. Areas close to good discoveries are marked as promising, while uncertain areas are identified as potentially interesting.
• Expected Improvement is the compass: it points toward the best tradeoff between exploring the unknown (areas where uncertainty is high) and exploiting the known (areas where the predicted mean is high). This compass guides you to where you are most likely to find gold, taking into account everything you already know.

At each iteration, you dig where the compass indicates, you update the map with the new result, and the compass redirects you to a new optimal location. Progressively, the map becomes more precise and the compass converges toward the best areas.

Python Implementation

Installing Libraries

pip install scikit-optimize numpy matplotlib scikit-learn optuna

Example 1: Minimizing a Test Function

Let’s start with the classic Branin function example, a two-variable function often used as a benchmark for Bayesian optimization:

import numpy as np
from skopt import gp_minimize
from skopt.plots import plot_convergence, plot_objective
import matplotlib.pyplot as plt

# Branin function (classic benchmark)
def branin(x):
    x1, x2 = x
    a = 1.0
    b = 5.1 / (4 * np.pi**2)
    c = 5.0 / np.pi
    r = 6.0
    s = 10.0
    t = 1.0 / (8 * np.pi)
    y = a * (x2 - b * x1**2 + c * x1 - r)**2 + s * (1 - t) * np.cos(x1) + s
    return y

# Bayesian optimization with gp_minimize
resultat = gp_minimize(
    branin,                       # function to minimize
    [(-5.0, 10.0), (0.0, 15.0)], # variable bounds
    n_calls=50,                   # number of evaluations
    random_state=42,              # reproducibility
    acq_func="EI",                # Expected Improvement
    noise="gaussian"              # noise model
)

print(f"Minimum found   : {resultat.fun:.4f}")
print(f"Optimal point   : {resultat.x}")
print(f"Number of calls : {resultat.n_iters}")

Typical results of this optimization:

Minimum found   : 0.3979
Optimal point   : [3.142, 2.275]
Number of calls : 50

The global minimum of the Branin function is approximately 0.3979. Bayesian Optimization finds it in only 50 evaluations, whereas a fine grid search would require several hundred.

Convergence Visualization

# Convergence tracking
fig, ax = plt.subplots(1, 2, figsize=(14, 5))

# Convergence plot
plot_convergence(resultat, ax=ax[0])
ax[0].set_title("Bayesian Optimization Convergence", fontsize=12)
ax[0].set_xlabel("Number of iterations")
ax[0].set_ylabel("Best value found")

# Function map and evaluated points
x1 = np.linspace(-5, 10, 200)
x2 = np.linspace(0, 15, 200)
X1, X2 = np.meshgrid(x1, x2)
Z = np.array([[branin([a, b]) for a, b in zip(r1, r2)] for r1, r2 in zip(X1, X2)])

contour = ax[1].contourf(X1, X2, Z, levels=50, cmap='viridis')
plt.colorbar(contour, ax=ax[1], label='Function value')
x_eval = np.array(resultat.x_iters)
ax[1].scatter(x_eval[:, 0], x_eval[:, 1], c='red', s=30, edgecolors='white', linewidths=1, zorder=5, label='Evaluated points')
ax[1].scatter(resultat.x[0], resultat.x[1], c='yellow', s=200, marker='*', edgecolors='black', linewidths=2, zorder=6, label='Minimum found')
ax[1].set_title("Search space and evaluated points", fontsize=12)
ax[1].set_xlabel("x₁")
ax[1].set_ylabel("x₂")
ax[1].legend()

plt.tight_layout()
plt.savefig("bayesian_optimization_convergence.png", dpi=150, bbox_inches="tight")
plt.show()

Example 2: RandomForest Hyperparameter Optimization

Now let’s examine the most common practical application: tuning the hyperparameters of a RandomForest classifier on a dataset:

from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import cross_val_score
from skopt.space import Integer, Real
from skopt import gp_minimize
import warnings
warnings.filterwarnings('ignore')

# Loading data
donnees = load_breast_cancer()
X, y = donnees.data, donnees.target

# Objective function to minimize (we minimize 1 - accuracy)
def objectif_rf(hyperparams):
    n_estimators, max_depth, min_samples_split, min_samples_leaf, max_features = hyperparams

    modele = RandomForestClassifier(
        n_estimators=int(n_estimators),
        max_depth=int(max_depth) if max_depth > 0 else None,
        min_samples_split=int(min_samples_split),
        min_samples_leaf=int(min_samples_leaf),
        max_features=max_features,
        random_state=42,
        n_jobs=-1
    )

    # 5-fold cross-validation
    scores = cross_val_score(modele, X, y, cv=5, scoring='accuracy')
    return 1.0 - scores.mean()  # We minimize the inverse of accuracy

# Defining the search space
espace_recherche = [
    Integer(50, 500, name='n_estimators'),
    Integer(3, 50, name='max_depth'),
    Integer(2, 20, name='min_samples_split'),
    Integer(1, 10, name='min_samples_leaf'),
    Real(0.1, 1.0, name='max_features')
]

print("Starting Bayesian Optimization for RandomForest...")
print("=" * 60)

resultat_rf = gp_minimize(
    objectif_rf,
    espace_recherche,
    n_calls=40,
    n_random_starts=10,
    random_state=42,
    acq_func="EI",
    verbose=True
)

# Displaying results
print("\n" + "=" * 60)
print("BAYESIAN OPTIMIZATION RESULTS")
print("=" * 60)
print(f"Best accuracy: {1.0 - resultat_rf.fun:.4f}")
print(f"\nBest hyperparameters:")
for i, dim in enumerate(espace_recherche):
    print(f"  {dim.name}: {resultat_rf.x[i]}")

# Comparison with default values
modele_defaut = RandomForestClassifier(random_state=42, n_jobs=-1)
scores_defaut = cross_val_score(modele_defaut, X, y, cv=5, scoring='accuracy')
print(f"\nDefault accuracy        : {scores_defaut.mean():.4f}")
print(f"Optimized accuracy      : {1.0 - resultat_rf.fun:.4f}")
print(f"Improvement             : {(1.0 - resultat_rf.fun - scores_defaut.mean())*100:.2f}%")

Detailed Convergence Tracking

import matplotlib.pyplot as plt
import numpy as np

# History of best values found
meilleures_valeurs = np.minimum.accumulate([1 - val for val in resultat_rf.func_vals])

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Evolution of the best value
axes[0].plot(meilleures_valeurs, 'b-', linewidth=2, marker='o', markersize=4)
axes[0].set_title("Bayesian Optimization Progress", fontsize=12)
axes[0].set_xlabel("Iteration")
axes[0].set_ylabel("Best accuracy")
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=scores_defaut.mean(), color='red', linestyle='--', alpha=0.7, label='Default baseline')
axes[0].legend()

# Distribution of explored values
axes[1].hist([1 - val for val in resultat_rf.func_vals], bins=20, color='steelblue', edgecolor='white', alpha=0.8)
axes[1].set_title("Distribution of explored performances", fontsize=12)
axes[1].set_xlabel("Accuracy")
axes[1].set_ylabel("Frequency")
axes[1].axvline(x=1.0 - resultat_rf.fun, color='red', linestyle='--', linewidth=2, label='Best value')
axes[1].axvline(x=scores_defaut.mean(), color='green', linestyle='--', linewidth=2, label='Baseline')
axes[1].legend()

plt.tight_layout()
plt.savefig("bayesian_rf_convergence.png", dpi=150, bbox_inches="tight")
plt.show()

Key Hyperparameters of Bayesian Optimization

Configuring Bayesian Optimization involves several hyperparameters that directly influence the quality and speed of convergence:

n_calls

The total number of evaluations of the objective function. This is the most important and most determining hyperparameter:

• Typical value: 30 to 100 for problems with 3 to 10 dimensions
• Increasing n_calls improves the result quality but increases computational cost
• Rule of thumb: approximately 10 × number of dimensions for good initial exploration

resultat = gp_minimize(fonction, espace_recherche, n_calls=60)

base_estimator

The probabilistic model used to approximate the objective function:

• “GP” (Gaussian Process, default): the standard choice, excellent for low-dimensional continuous spaces (≤ 20)
• “RF” (Random Forest): robust, handles categorical variables and mixed spaces better
• “ET” (Extra Trees): a variant of Random Forest, sometimes faster
• “GBRT” (Gradient Boosted Regression Trees): performant for high-dimensional spaces

resultat = gp_minimize(fonction, espace_recherche, base_estimator="RF")

acq_func

The acquisition function that guides exploration of the search space:

• “EI” (Expected Improvement, default): the best general tradeoff, naturally balances exploration and exploitation
• “LCB” (Lower Confidence Bound): more exploratory, useful when worried about getting stuck in a local optimum
• “PI” (Probability of Improvement): more exploitative, converges quickly but risks missing the global optimum
• “EIps” (Expected Improvement per Second): a variant that accounts for evaluation time

resultat = gp_minimize(fonction, espace_recherche, acq_func="EI", acq_optimizer="auto")

kernel

The Gaussian process kernel determines the assumed correlation structure between points in the search space:

• Matérn 5/2 (default in scikit-optimize): good tradeoff, allows non-infinitely differentiable functions
• RBF (Radial Basis Function): assumes very smooth functions, converges fast if the hypothesis holds
• Matérn 3/2: rougher, suited to functions with abrupt variations

from skopt.learning import GaussianProcessRegressor
from scipy.optimize import minimize
noyau = "Matern52"  # or "RBF", "Matern32"
resultat = gp_minimize(fonction, espace_recherche, kernel=noyau)

n_random_starts

The number of randomly sampled points before starting the actual Bayesian optimization. These initial points allow the Gaussian process to have a reasonable first estimate of the function:

• Typical value: 5 to 20 (approximately 10 by default)
• Too few: the initial model can be misleading
• Too many: evaluations are wasted without Bayesian guidance

resultat = gp_minimize(fonction, espace_recherche, n_random_starts=10, n_calls=50)

Advantages and Limitations

Advantages of Bayesian Optimization

1. Sample efficiency: finds competitive solutions with 5 to 10 times fewer evaluations than random search. This is the primary and most decisive advantage.
2. History awareness: each evaluation informs the next. Unlike grid search or random search which are memoryless methods, Bayesian Optimization learns from each trial and progressively refines its strategy.
3. Natural exploration-exploitation tradeoff management: the acquisition function automatically balances exploring new regions and exploiting promising areas, without manual intervention.
4. Uncertainty estimation: the Gaussian process provides not only a prediction but also a confidence interval, allowing us to quantify our certainty about each region of the search space.
5. Compatible with noisy functions: the Gaussian process noise parameter handles non-deterministic evaluations, which are common in machine learning.

Limitations of Bayesian Optimization

1. Gaussian process computational cost: inverting the covariance matrix K has O(n³) complexity where n is the number of observations. This becomes prohibitive beyond a few thousand points.
2. Curse of dimensionality: Bayesian Optimization works well for spaces of 1 to ~20 dimensions. Beyond that, the Gaussian process struggles to model the function effectively and tree-based methods (SMAC, TPE) become preferable.
3. Smoothness assumption: the Gaussian process assumes the function is relatively smooth. If the objective function has discontinuities or very abrupt variations, the model can be misleading.
4. Sensitivity to kernel hyperparameters: the choice and optimization of kernel hyperparameters significantly influence the results. A poor kernel can lead to premature convergence to a local optimum.

4 Concrete Use Cases

Use Case 1: Deep Neural Network Hyperparameter Tuning

Training a deep neural network can take hours or even days. Each evaluation of a hyperparameter configuration (learning rate, number of layers, number of neurons, dropout rate, etc.) represents a considerable time investment. Bayesian Optimization drastically reduces the number of training runs needed to find a performant configuration.

# Typical architecture for a neural network
espace_recherche_reseau = [
    Real(1e-4, 1e-1, prior='log-uniform', name='learning_rate'),
    Integer(32, 512, name='hidden_units'),
    Integer(2, 8, name='num_layers'),
    Real(0.1, 0.8, name='dropout_rate'),
    Integer(16, 256, name='batch_size'),
    Real(1e-5, 1e-2, prior='log-uniform', name='weight_decay')
]

Use Case 2: Physical Model Calibration and Numerical Simulations

In engineering, calibrating physical models (fluid mechanics, thermodynamics, materials science) requires extremely expensive numerical simulations. Bayesian Optimization finds the physical parameters that minimize the gap between simulation and experimental measurement with a minimal number of simulations. Companies like Airbus and BMW use this approach to optimize the shape of their aerodynamic components.

Use Case 3: Data Processing Pipeline Optimization

Production data pipelines have many parameters: filtering thresholds, sliding window sizes, sampling rates, preprocessing parameters. Bayesian Optimization allows optimizing the entire pipeline end-to-end by maximizing a business metric (conversion rate, revenue per user, detection accuracy).

Use Case 4: Molecule Design and Drug Discovery

In computational chemistry and drug discovery, each evaluation corresponds to a molecular dynamics simulation or an expensive laboratory experiment. Bayesian Optimization guides the search toward the most promising molecules by maximizing target properties (binding affinity, solubility, reduced toxicity) while minimizing the number of experiments needed. Platforms like machine learning-assisted drug discovery use this approach daily.

See Also

• How to Implement Binary Search in Python: A Step-by-Step Guide
• Energize Your Python Projects: Simulating Amoebas in a 3D Grid with Python