Bayesian Optimization: Complete Guide — Principles, Examples, and Python Implementation

Summary

Bayesian Optimization is a sequential optimization method designed for expensive-to-evaluate functions, typically in the absence of gradients. Unlike Random Search, which blindly explores the parameter space, Bayesian Optimization builds a probabilistic model — a Gaussian Process — that approximates the unknown objective function. At each iteration, an acquisition function (Expected Improvement, UCB, Probability of Improvement) determines the most promising point to evaluate. This “model → acquire → update” cycle converges much faster to the optimum than blind methods, especially when each function evaluation takes minutes or hours (training a neural network, physical simulation, etc.).

Mathematical Principles

Bayesian optimization rests on two mathematical pillars: the surrogate model and the acquisition function.

1. The Surrogate Model: Gaussian Process

The core idea is to model the unknown objective function f(x) as a Gaussian Process:

f(x) ∼ GP(m(x), k(x, x’))

where:

• m(x) is the mean function (often chosen as m(x) = 0 after centering the data).
• k(x, x’) is the covariance function (or kernel). The most common choice is the RBF kernel (Radial Basis Function), also called the Gaussian kernel:

k(x, x’) = σ² exp( −∥x − x’∥² / (2ℓ²) )

The parameter ℓ (correlation length) determines the scale at which two points are considered correlated, while σ² controls the overall variance.

After observing n points D_n = {(x_i, y_i)}{i=1}^n where y_i = f(x_i) + ε_i (with Gaussian noise ε_i ∼ N(0, σ_n²)), the posterior distribution of f at a new point x* is a normal distribution:

f(x_) | D_n, x_ ∼ N(μ(x_), σ²(x_))

with:

μ(x_) = k_ᵀ (K + σ_n² I)⁻¹ Y

σ²(x_) = k(x_, x_) − k_ᵀ (K + σn² I)⁻¹ k*

where K is the covariance matrix K_{ij} = k(x_i, x_j), k_ is the vector of covariances between x_ and the observed points, and Y is the vector of observations.

These two quantities are the heart of the process: μ(x_) is a prediction of the function value, and σ(x_) measures the uncertainty of that prediction.

2. The Acquisition Function

The acquisition function α(x) transforms the GP predictive distribution into a single score that guides the choice of the next point to evaluate. The three most commonly used acquisition functions are:

Expected Improvement (EI):

EI(x) = E[max(f_min − f(x), 0)]

where f_min denotes the best minimum observed so far. Under the Gaussian predictive distribution, this expectation has a closed-form expression:

EI(x) = (f_min − μ(x)) Φ(Z) + σ(x) φ(Z)

where Z = (f_min − μ(x)) / σ(x), Φ is the standard normal cumulative distribution function, and φ is its density.

The term (f_min − μ(x)) Φ(Z) favors exploitation: it grows when the prediction μ(x) is lower than the best known minimum. The term σ(x) φ(Z) favors exploration: it grows when uncertainty is high.

Upper Confidence Bound (UCB):

UCB(x) = μ(x) + κ · σ(x)

The parameter κ controls the exploration/exploration tradeoff: a large κ favors exploration (uncertain regions), while a small κ favors exploitation (already promising regions).

Probability of Improvement (PI):

PI(x) = P(f(x) < f_min) = Φ((f_min − μ(x)) / σ(x))

PI tends to be more exploitative than EI and is therefore less commonly used in practice.

3. The Iteration Loop

The algorithm proceeds as follows:

1. Initialization: evaluate f at n_initial randomly drawn points (or via a Latin Hypercube Sampling design of experiments).
2. GP fitting: update the Gaussian Process with the available observations.
3. Acquisition optimization: find x_next = argmax_x α(x). This sub-optimization is usually done with a multi-restart local optimizer like L-BFGS-B.
4. Evaluation: compute y_next = f(x_next) (this is the expensive step).
5. Update: add (x_next, y_next) to D_n and return to step 2.
6. Stopping: after n_calls evaluations or upon convergence.

Intuition: The Geologist Looking for Oil

Imagine a geologist who must decide where to drill to find oil. Each drill costs millions of euros, so strategy is essential.

The brute-force method would be to dig at purely random locations (Random Search) or on a regular grid (Grid Search). This is extremely costly and inefficient.

Bayesian Optimization works differently. The geologist starts with a few exploratory drills across the region. With these few probes, they progressively build a map of the subsurface (the Gaussian Process): some zones appear rich in oil (low μ if minimizing, or high if maximizing), while others are still completely unknown (high σ).

Then, at each step, the geologist uses the acquisition function to decide where to drill next. They can choose to:

• Exploit: drill where the map predicts the most oil (favorable μ).
• Explore: drill where uncertainty is high (high σ), because something unexpected might be discovered there.
• Combine both via EI or UCB, which intelligently balance these two objectives.

This is exactly what Bayesian optimization does for hyperparameters: each model training run is an expensive drill and the GP is the map that guides future choices.

Python Implementation

We present two implementations: with scikit-optimize for a direct pedagogical approach, and with Optuna for modern production use.

Approach 1: scikit-optimize on a test function

import numpy as np
import matplotlib.pyplot as plt
from skopt import gp_minimize
from skopt.space import Real

# Test function (Branin)
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)
    val = a * (x2 - b * x1**2 + c * x1 - r)**2
    val += s * (1 - t) * np.cos(x1) + s
    return val

# Search space definition
space = [
    Real(-5.0, 10.0, name="x1"),
    Real(0.0, 15.0, name="x2")
]

# Bayesian optimization
result = gp_minimize(
    func=branin,
    dimensions=space,
    n_calls=50,
    n_initial_points=10,
    acq_func="EI",
    random_state=42,
    verbose=True
)

print(f"Best value found: {result.fun:.6f}")
print(f"Best parameters: x1={result.x[0]:.4f}, x2={result.x[1]:.4f}")

# Convergence visualization
from skopt.plots import plot_convergence
fig, ax = plt.subplots(figsize=(8, 5))
plot_convergence(result)
plt.savefig("bayesian_optimization_convergence.png", dpi=150)
plt.close()

# Comparison with Random Search
from skopt.dummy import dummy_minimize

result_random = dummy_minimize(
    func=branin,
    dimensions=space,
    n_calls=50,
    random_state=42
)

print("\nComparison complete. BO generally converges")
print("5 to 10 times faster than random search.")

Approach 2: Optuna for optimizing a Random Forest

import optuna
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import cross_val_score

# Load data
X, y = load_breast_cancer(return_X_y=True)

# Objective function for Optuna
def objective(trial):
    # Hyperparameters to optimize
    n_estimators = trial.suggest_int("n_estimators", 50, 500)
    max_depth = trial.suggest_int("max_depth", 3, 30)
    min_samples_split = trial.suggest_int("min_samples_split", 2, 20)
    min_samples_leaf = trial.suggest_int("min_samples_leaf", 1, 10)
    max_features = trial.suggest_categorical(
        "max_features", ["sqrt", "log2", None]
    )

    # Create and evaluate model
    model = RandomForestClassifier(
        n_estimators=n_estimators,
        max_depth=max_depth,
        min_samples_split=min_samples_split,
        min_samples_leaf=min_samples_leaf,
        max_features=max_features,
        random_state=42,
        n_jobs=-1
    )

    score = cross_val_score(
        model, X, y, cv=5, scoring="accuracy"
    ).mean()
    return score

# Create study with TPE sampler (Tree-structured
# Parzen Estimator) — Optuna's Bayesian approach
study = optuna.create_study(
    direction="maximize",
    sampler=optuna.samplers.TPESampler(seed=42)
)

study.optimize(objective, n_trials=50, show_progress_bar=True)

print(f"\nBest accuracy: {study.best_value:.4f}")
print(f"Best hyperparameters: {study.best_params}")

Hyperparameters and Configuration

The key parameters of Bayesian Optimization directly influence its performance:

The choice of base_estimator is particularly important. The Gaussian Process is the default and most mathematically elegant choice, but it has O(n³) complexity that makes it costly beyond a few hundred observation points. For high-dimensional spaces (d > 20) or a large number of evaluations, Random Forests (base_estimator="RF") or Gradient Boosting Regression Trees (base_estimator="GBRT") are more scalable alternatives.

Advantages and Limitations

Advantages

• Remarkable efficiency on expensive functions: converges in 10 to 50 evaluations where Grid Search requires thousands.
• No gradient required: works with any black-box function.
• Natural uncertainty incorporation: explicitly quantifies it via GP variance.
• Natively handles continuous, categorical, and integer parameters (with Optuna).
• Solid probabilistic theory: unlike heuristic methods, convergence guarantees exist.
• Parallelization possible: some variants (q-EI, batch BO) allow evaluating multiple points simultaneously.

Limitations

• Cubic complexity: inverting the GP covariance matrix is O(n³), limiting the number of practical evaluations to a few hundred with a standard GP.
• Sensitive to kernel choice: a poor kernel (or poor kernel hyperparameters) can degrade performance.
• Reduced performance in very high dimensions: beyond 20-30 dimensions, the curse of dimensionality strikes. Specialized variants like REMBO (Random Embeddings for BO) or TuRBO exist but are more complex.
• Costly acquisition optimization: maximizing α(x) is itself a non-convex optimization problem that can get trapped in local optima.
• Less efficient than gradient-based methods when gradients are available (which is rarely the case for hyperparameter optimization).

4 Concrete Use Cases

1. Machine Learning Model Hyperparameter Optimization

This is the most classic application. Tuning the hyperparameters of an XGBoost, Random Forest, or deep neural network involves searching for an optimum in a space where each evaluation costs minutes to hours of training. Bayesian Optimization is particularly well suited because it requires far fewer evaluations than Grid Search or Random Search. Optuna, based on the TPE sampler (a cousin of the GP), has become the industry standard tool for this use case.

2. Physical and Numerical Simulation Optimization

Fluid dynamics simulations (CFD), particle collisions, or chemical reactions are extremely expensive (sometimes days on a supercomputer). Bayesian optimization finds optimal parameters (wing shape, reaction temperature, alloy composition) in a minimal number of simulations.

3. Scientific Experiment Design (Design of Experiments)

In chemistry, biology, or materials science, each real-world experiment costs time and money. Bayesian optimization progressively guides the experimenter toward optimal conditions: reactant concentration, temperature, pH, pressure. This is the foundation of the Bayesian Experimental Design approach.

4. System Tuning and Performance Tuning

Tuning database parameters (PostgreSQL, MySQL), search engines (Elasticsearch), or Kubernetes clusters involves testing configurations and measuring performance metrics. Each test can take hours of benchmarking. Bayesian Optimization finds the optimal configuration in just a few dozen trials. Tools like OtterTune use precisely this approach for automatic database tuning.

See Also

• Explore the Goldbach Conjecture in Python: Practical Guide for Programming Enthusiasts
• Solving Diophantine Equations III with Python: Complete Guide and Tips