Bayesian Optimization : Guide complet — Principes, Exemples et Implémentation Python

Résumé

La Bayesian Optimization (optimisation bayésienne) est une méthode d’optimisation séquentielle conçue pour les fonctions coûteuses à évaluer, typiquement en l’absence de gradient. Contrairement à la recherche aléatoire (Random Search) qui explore l’espace des paramètres aveuglément, la Bayesian Optimization construit un modèle probabiliste — un processus gaussien (Gaussian Process) — qui approxime la fonction objectif inconnue. À chaque itération, une fonction d’acquisition (Expected Improvement, UCB, Probability of Improvement) détermine le point le plus prometteur à évaluer. Ce cycle « modéliser → acquérir → mettre à jour » converge bien plus rapidement vers l’optimum que les méthodes aveugles, en particulier quand chaque évaluation de la fonction prend plusieurs minutes ou heures (entraînement d’un réseau de neurones, simulation physique, etc.).

Principe mathématique

L’optimisation bayésienne repose sur deux piliers mathématiques : le modèle de substitution (surrogate model) et la fonction d’acquisition.

1. Le surrogate model : Processus Gaussien

L’idée centrale est de modéliser la fonction objectif inconnue f(x) comme un processus gaussien :

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

où :

• m(x) est la fonction moyenne (souvent choisie comme m(x) = 0 après centrage des données).
• k(x, x’) est la fonction de covariance (ou kernel). Le choix le plus courant est le RBF kernel (Radial Basis Function), aussi appelé noyau gaussien :

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

Le paramètre ℓ (longueur de corrélation) détermine l’échelle à laquelle deux points sont considérés comme corrélés, tandis que σ² contrôle la variance globale.

Après avoir observé n points D_n = {(x_i, y_i)}{i=1}^n où y_i = f(x_i) + ε_i (avec un bruit gaussien ε_i ∼ N(0, σ_n²)), la distribution a posteriori de f en un nouveau point x* est une loi normale :

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

avec :

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

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

où K est la matrice de covariance K_{ij} = k(x_i, x_j), k_ est le vecteur des covariances entre x_ et les points observés, et Y est le vecteur des observations.

Ces deux grandeurs constituent le cœur du processus : μ(x_) est une prédiction de la valeur de la fonction, et σ(x_) mesure l’incertitude de cette prédiction.

2. La fonction d’acquisition

La fonction d’acquisition α(x) transforme la distribution prédictive du GP en un score unique qui guide le choix du prochain point à évaluer. Les trois fonctions d’acquisition les plus utilisées sont :

Expected Improvement (EI) :

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

où f_min désigne le meilleur minimum observé jusqu’ici. Sous la prédictive gaussienne, cette espérance s’exprime analytiquement :

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

où Z = (f_min − μ(x)) / σ(x), Φ est la fonction de répartition de la loi normale standard, et φ sa densité.

Le terme (f_min − μ(x)) Φ(Z) favorise l’exploitation : il croît lorsque la prédiction μ(x) est inférieure au meilleur minimum connu. Le terme σ(x) φ(Z) favorise l’exploration : il croît lorsque l’incertitude est élevée.

Upper Confidence Bound (UCB) :

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

Le paramètre κ contrôle le compromis exploration/exploitation : une grande valeur de κ favorise l’exploration (régions incertaines), tandis qu’une petite valeur favorise l’exploitation (régions déjà prometteuses).

Probability of Improvement (PI) :

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

PI tend à être plus exploitative que EI et est donc moins utilisée en pratique.

3. La boucle d’itération

L’algorithme se déroule comme suit :

1. Initialisation : évaluer f en n_initial points tirés aléatoirement (ou via un plan d’expériences de type Latin Hypercube Sampling).
2. Ajustement du GP : mettre à jour le processus gaussien avec les observations disponibles.
3. Optimisation de l’acquisition : trouver x_next = argmax_x α(x). Cette sous-optimisation se fait généralement avec un optimiseur local multi-redémarrages comme L-BFGS-B.
4. Évaluation : calculer y_next = f(x_next) (c’est l’étape coûteuse).
5. Mise à jour : ajouter (x_next, y_next) à D_n et retourner à l’étape 2.
6. Arrêt : après n_calls évaluations ou sur convergence.

Intuition : le géologue qui cherche du pétrole

Imaginez un géologue qui doit déterminer où forer pour trouver du pétrole. Chaque forage coûte des millions d’euros, il faut donc être stratégique.

La méthode brute serait de creuser à des endroits purement aléatoires (Random Search) ou sur une grille régulière (Grid Search). C’est extrêmement coûteux et inefficace.

La Bayesian Optimization fonctionne différemment. Le géologue commence par quelques forages d’exploration à travers la région. Avec ces quelques sondages, il construit progressivement une carte du sous-sol (le processus gaussien) : certaines zones apparaissent riches en pétrole (faible μ si l’on minimise, ou élevé si l’on maximise), d’autres sont encore totalement inconnues (σ élevé).

Ensuite, à chaque étape, le géologue utilise la fonction d’acquisition pour décider où forer ensuite. Il peut choisir :

• D’exploiter : forer là où la carte prédit le plus de pétrole (μ favorable).
• D’explorer : forer là où l’incertitude est grande (σ élevé), car on y découvrira potentiellement quelque chose d’inattendu.
• De combiner les deux via EI ou UCB, qui pondèrent intelligemment ces deux objectifs.

C’est exactement ce que fait l’optimisation bayésienne pour les hyperparamètres : chaque entraînement de modèle est un forage coûteux et le GP est la carte qui guide les choix futurs.

Implémentation Python

Nous présentons deux implémentations : avec scikit-optimize pour une approche pédagogique directe, et avec Optuna pour un usage moderne en production.

Approche 1 : scikit-optimize sur une fonction de test

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

# Fonction de test (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

# Définition de l'espace de recherche
space = [
    Real(-5.0, 10.0, name="x1"),
    Real(0.0, 15.0, name="x2")
]

# Optimisation bayésienne
result = gp_minimize(
    func=branin,
    dimensions=space,
    n_calls=50,
    n_initial_points=10,
    acq_func="EI",
    random_state=42,
    verbose=True
)

print(f"Meilleure valeur trouvée : {result.fun:.6f}")
print(f"Meilleurs paramètres : x1={result.x[0]:.4f}, x2={result.x[1]:.4f}")

# Visualisation de la convergence
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()

# Comparaison avec Random Search
from skopt.dummy import dummy_minimize

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

print("\nComparaison terminée. La BO converge en général")
print("5 à 10 fois plus vite qu'une recherche aléatoire.")

Approche 2 : Optuna pour l’optimisation d’un Random Forest

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

# Chargement des données
X, y = load_breast_cancer(return_X_y=True)

# Fonction objectif pour Optuna
def objective(trial):
    # Hyperparamètres à optimiser
    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]
    )

    # Création et évaluation du modèle
    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

# Création de l'étude avec un sampler TPE (Tree-structured
# Parzen Estimator) — l'approche bayésienne d'Optuna
study = optuna.create_study(
    direction="maximize",
    sampler=optuna.samplers.TPESampler(seed=42)
)

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

print(f"\nMeilleure accuracy : {study.best_value:.4f}")
print(f"Meilleurs hyperparamètres : {study.best_params}")

Hyperparamètres et configuration

Les paramètres clés de la Bayesian Optimization influencent directement sa performance :

Le choix du base_estimator est particulièrement important. Le Gaussian Process est le choix par défaut et le plus élégant mathématiquement, mais il a une complexité en O(n³) qui le rend coûteux au-delà de quelques centaines de points d’observation. Pour les espaces de grande dimension (d > 20) ou un grand nombre d’évaluations, les Random Forests (base_estimator="RF") ou les Gradient Boosting Regression Trees (base_estimator="GBRT") sont des alternatives plus scalables.

Avantages et limites

Avantages

• Efficacité remarquable sur les fonctions coûteuses : converge en 10 à 50 évaluations là où le Grid Search en demande des milliers.
• Ne nécessite pas de gradient : fonctionne avec n’importe quelle fonction boîte noire (black-box).
• Incorporation naturelle de l’incertitude : la quantifie explicitement via la variance du GP.
• Gère les paramètres continus, catégoriels et entiers de manière native (avec Optuna).
• Théorie probabiliste solide : contrairement aux méthodes heuristiques, on dispose de garanties de convergence.
• Parallélisation possible : certaines variantes (q-EI, batch BO) permettent d’évaluer plusieurs points simultanément.

Limites

• Complexité cubique : l’inversion de la matrice de covariance du GP est en O(n³), ce qui limite le nombre d’évaluations pratiques à quelques centaines avec un GP classique.
• Sensible au choix du kernel : un mauvais noyau (ou de mauvais hyperparamètres de noyau) peut dégrader les performances.
• Performance réduite en très haute dimension : au-delà de 20-30 dimensions, la malédiction de la dimensionalité frappe. Des variantes spécialisées comme REMBO (Random Embeddings for BO) ou TuRBO existent mais sont plus complexes.
• Optimisation de l’acquisition coûteuse : maximiser α(x) est lui-même un problème d’optimisation non convexe, qui peut piéger dans des optimums locaux.
• Moins efficace que les méthodes basées gradient quand les gradients sont disponibles (ce qui est rarement le cas pour l’optimisation d’hyperparamètres).

4 cas d’usage concrets

1. Optimisation d’hyperparamètres de modèles de Machine Learning

C’est l’application la plus classique. Régler les hyperparamètres d’un XGBoost, d’un Random Forest ou d’un réseau de neurones profond revient à chercher un optimum dans un espace où chaque évaluation coûte plusieurs minutes à plusieurs heures d’entraînement. La Bayesian Optimization est particulièrement adaptée car elle nécessite beaucoup moins d’évaluations que le Grid Search ou le Random Search. Optuna, basé sur le sampler TPE (un cousin du GP), est devenu l’outil standard de l’industrie pour ce cas d’usage.

2. Optimisation de simulations physiques et numériques

Les simulations de dynamique des fluides (CFD), de collision de particules, ou de réactions chimiques sont extrêmement coûteuses (parfois des jours sur un supercalculateur). L’optimisation bayésienne permet de trouver les paramètres optimaux (forme d’une aile, température de réaction, composition d’un alliage) en un nombre minimal de simulations.

3. Conception d’expériences scientifiques (Design of Experiments)

En chimie, en biologie ou en science des matériaux, chaque expérience réelle coûte du temps et de l’argent. L’optimisation bayésienne guide progressivement l’expérimentateur vers les conditions optimales : concentration de réactif, température, pH, pression. C’est le fondement de l’approche Bayesian Experimental Design.

4. Réglage de systèmes et tuning de performances

Le réglage des paramètres d’une base de données (PostgreSQL, MySQL), d’un moteur de recherche (Elasticsearch), ou d’un cluster Kubernetes implique de tester des configurations et de mesurer des métriques de performance. Chaque test peut prendre des heures de benchmark. La Bayesian Optimization trouve la configuration optimale en quelques dizaines d’essais seulement. Des outils comme OtterTune utilisent précisément cette approche pour le database tuning automatique.

Voir aussi

• Explorez la Conjecture de Goldbach en Python : Guide Pratique pour les Passionnés de Programmation
• Résolution des Équations Diophantiennes III avec Python : Guide Complet et Astuces