Optimisation Bayésienne

Résumé

L’Optimisation Bayésienne est une approche puissante et mathématiquement élégante pour optimiser des fonctions coûteuses à évaluer, notamment dans le domaine du réglage des hyperparamètres en apprentissage automatique. Contrairement aux méthodes classiques telles que la recherche par grille (grid search) ou la recherche aléatoire (random search), l’Optimisation Bayésienne construit un modèle probabiliste de la fonction objectif et l’utilise de manière stratégique pour décider où évaluer la fonction ensuite. Cette approche permet de trouver des solutions de haute qualité avec un nombre d’évaluations considérablement réduit, ce qui représente un avantage décisif lorsque chaque évaluation nécessite des heures d’entraînement de modèle.

Principe Mathématique

Le fondement mathématique de l’Optimisation Bayésienne repose sur deux concepts essentiels : la modélisation de la fonction objectif et la stratégie d’acquisition du point suivant.

Processus Gaussien

La première étape consiste à modéliser la fonction objectif inconnue f(x) comme un processus gaussien. Mathématiquement, on suppose que la distribution a posteriori de f sachant les observations X suit une loi normale multidimensionnelle :

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

où μ représente le vecteur des moyennes a priori et K est la matrice de covariance calculée à l’aide d’un noyau (kernel), généralement le noyau RBF (fonction à base radiale). Le noyau RBF s’exprime de la manière suivante :

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

où σ² est la variance et l est la longueur caractéristique (length scale) qui contrôle la régularité de la fonction. Cette hypothèse de processus gaussien signifie que toute combinaison linéaire finie de valeurs de f suit une distribution gaussienne. Cela permet de calculer de manière analytique la distribution a posteriori après chaque nouvelle observation.

Lorsqu’on dispose de n observations (x₁, y₁), …, (xₙ, yₙ), la prédiction en un nouveau point x* s’obtient par les formules suivantes :

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

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

où k est le vecteur des covariances entre x et les points d’entraînement, K est la matrice de covariance calculée sur les observations existantes, et σₙ² représente le bruit de mesure.

Fonction d’Acquisition : Expected Improvement

Une fois le processus gaussien établi, il faut décider où évaluer la fonction ensuite. C’est le rôle de la fonction d’acquisition. La plus populaire est l’Expected Improvement (amélioration attendue), notée EI(x) :

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

où :

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

Dans cette formule, μ(x) est la moyenne prédite par le processus gaussien au point x, f(x_best) est la meilleure valeur observée jusqu’à présent, σ(x) est l’écart-type prédit, ξ (xi) est un paramètre d’exploration qui contrôle le compromis exploration-exploitation, Φ désigne la fonction de répartition de la loi normale standard et φ est la densité de probabilité de cette même loi.

On choisit le point x qui maximise EI(x) pour la prochaine évaluation. Ce point représente le meilleur compromis entre explorer des régions incertaines (σ(x) élevé) et exploiter des régions prometteuses (μ(x) élevé).

D’autres fonctions d’acquisition existent également :

• Upper Confidence Bound (UCB) : μ(x) + κ · σ(x), où κ contrôle le niveau d’exploration
• Probability of Improvement (PI) : Φ(Z), plus simple mais moins efficace en pratique
• Expected Quantile Improvement : variante robuste au bruit

Intuition

L’Optimisation Bayésienne est comme un chercheur d’or intelligent. Imaginez que vous devez trouver la veine d’or la plus riche dans un vaste territoire, mais chaque coup de pelle coûte cher en temps et en effort.

Avec la recherche par grille (grid search), vous creusez systématiquement à intervalles réguliers, comme un quadrillage. Vous couvrez tout, mais vous passez énormément de temps dans des zones stériles. C’est méthodique mais extrêmement inefficace.

Avec la recherche aléatoire (random search), vous creusez n’importe où, en espérant tomber sur une bonne zone. C’est mieux que la grille dans les espaces de grande dimension, mais vous ne tirez aucune leçon de vos coups de pelle précédents.

L’Optimisation Bayésienne, elle, agit comme un prospecteur expérimenté. Après chaque coup de pelle, elle construit une carte des zones prometteuses basée sur tout ce qu’elle a appris jusqu’à présent :

• Le processus gaussien est la carte : à partir des quelques points évalués, il estime non seulement la valeur probable de l’or en chaque endroit non exploré, mais aussi l’incertitude associée. Les zones proches des bonnes découvertes sont marquées comme prometteuses, tandis que les zones incertaines sont identifiées comme potentiellement intéressantes.
• L’Expected Improvement est le compas : il pointe vers le meilleur compromis entre explorer l’inconnu (les zones où l’incertitude est grande) et exploiter le connu (les zones où la moyenne prédite est élevée). Ce compas vous guide vers l’endroit où vous avez le plus de chances de trouver de l’or, en tenant compte de tout ce que vous savez déjà.

À chaque itération, vous creusez là où le compas indique, vous mettez à jour la carte avec le nouveau résultat, et le compas vous redirige vers un nouvel endroit optimal. Progressivement, la carte se précise et le compas converge vers les meilleures zones.

Implémentation Python

Installation des bibliothèques

pip install scikit-optimize numpy matplotlib scikit-learn optuna

Exemple 1 : Minimisation d’une fonction de test

Commençons par l’exemple classique de la fonction de Branin, une fonction à deux variables souvent utilisée comme benchmark pour l’optimisation bayésienne :

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

# Fonction de Branin (benchmark classique)
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

# Optimisation bayésienne avec gp_minimize
resultat = gp_minimize(
    branin,                       # fonction à minimiser
    [(-5.0, 10.0), (0.0, 15.0)], # bornes des variables
    n_calls=50,                   # nombre d'évaluations
    random_state=42,              # reproductibilité
    acq_func="EI",                # Expected Improvement
    noise="gaussian"              # modèle de bruit
)

print(f"Minimum trouvé : {resultat.fun:.4f}")
print(f"Point optimal    : {resultat.x}")
print(f"Nombre d'appels  : {resultat.n_iters}")

Résultats typiques de cette optimisation :

Minimum trouvé : 0.3979
Point optimal    : [3.142, 2.275]
Nombre d'appels  : 50

Le minimum global de la fonction de Branin vaut environ 0,3979. L’Optimisation Bayésienne le trouve en seulement 50 évaluations, alors qu’une recherche par grille fine en nécessiterait plusieurs centaines.

Visualisation de la convergence

# Suivi de la convergence
fig, ax = plt.subplots(1, 2, figsize=(14, 5))

# Tracé de convergence
plot_convergence(resultat, ax=ax[0])
ax[0].set_title("Convergence de l'Optimisation Bayésienne", fontsize=12)
ax[0].set_xlabel("Nombre d'itérations")
ax[0].set_ylabel("Valeur optimale trouvée")

# Carte de la fonction et des points évalués
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='Valeur de la fonction')
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='Points évalués')
ax[1].scatter(resultat.x[0], resultat.x[1], c='yellow', s=200, marker='*', edgecolors='black', linewidths=2, zorder=6, label='Minimum trouvé')
ax[1].set_title("Espace de recherche et points évalués", 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()

Exemple 2 : Optimisation d’hyperparamètres RandomForest

Examinons maintenant l’application pratique la plus courante : le réglage des hyperparamètres d’un classifieur RandomForest sur un ensemble de données :

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')

# Chargement des données
donnees = load_breast_cancer()
X, y = donnees.data, donnees.target

# Fonction objectif à minimiser (on minimise 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
    )

    # Validation croisée à 5 plis
    scores = cross_val_score(modele, X, y, cv=5, scoring='accuracy')
    return 1.0 - scores.mean()  # On minimise l'inverse de l'accuracy

# Définition de l'espace de recherche
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("Démarrage de l'Optimisation Bayésienne pour 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
)

# Affichage des résultats
print("\n" + "=" * 60)
print("RÉSULTATS DE L'OPTIMISATION BAYÉSIENNE")
print("=" * 60)
print(f"Meilleure accuracy : {1.0 - resultat_rf.fun:.4f}")
print(f"\nMeilleurs hyperparamètres :")
for i, dim in enumerate(espace_recherche):
    print(f"  {dim.name}: {resultat_rf.x[i]}")

# Comparaison avec les valeurs par défaut
modele_defaut = RandomForestClassifier(random_state=42, n_jobs=-1)
scores_defaut = cross_val_score(modele_defaut, X, y, cv=5, scoring='accuracy')
print(f"\nAccuracy par défaut        : {scores_defaut.mean():.4f}")
print(f"Accuracy après optimisation : {1.0 - resultat_rf.fun:.4f}")
print(f"Amélioration               : {(1.0 - resultat_rf.fun - scores_defaut.mean())*100:.2f}%")

Suivi de convergence détaillé

import matplotlib.pyplot as plt
import numpy as np

# Historique des meilleures valeurs trouvées
meilleures_valeurs = np.minimum.accumulate([1 - val for val in resultat_rf.func_vals])

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

# Évolution de la meilleure valeur
axes[0].plot(meilleures_valeurs, 'b-', linewidth=2, marker='o', markersize=4)
axes[0].set_title("Progression de l'Optimisation Bayésienne", fontsize=12)
axes[0].set_xlabel("Itération")
axes[0].set_ylabel("Meilleure accuracy")
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=scores_defaut.mean(), color='red', linestyle='--', alpha=0.7, label='Baseline par défaut')
axes[0].legend()

# Distribution des valeurs explorées
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 des performances explorées", fontsize=12)
axes[1].set_xlabel("Accuracy")
axes[1].set_ylabel("Fréquence")
axes[1].axvline(x=1.0 - resultat_rf.fun, color='red', linestyle='--', linewidth=2, label='Meilleure valeur')
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()

Hyperparamètres Clés de l’Optimisation Bayésienne

La configuration de l’Optimisation Bayésienne implique plusieurs hyperparamètres qui influencent directement la qualité et la rapidité de la convergence :

n_calls

Le nombre total d’évaluations de la fonction objectif. C’est l’hyperparamètre le plus important et le plus déterminant :

• Valeur typique : 30 à 100 pour des problèmes de 3 à 10 dimensions
• Augmenter n_calls améliore la qualité du résultat mais augmente le coût computationnel
• Règle empirique : environ 10 × nombre de dimensions pour une bonne exploration initiale

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

base_estimator

Le modèle probabiliste utilisé pour approximer la fonction objectif :

• “GP” (Gaussian Process, par défaut) : le choix standard, excellent pour les espaces continus de petite dimension (≤ 20)
• “RF” (Random Forest) : robuste, gère mieux les variables catégorielles et les espaces mixtes
• “ET” (Extra Trees) : variante du Random Forest, parfois plus rapide
• “GBRT” (Gradient Boosted Regression Trees) : performant pour les espaces de grande dimension

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

acq_func

La fonction d’acquisition qui guide l’exploration de l’espace de recherche :

• “EI” (Expected Improvement, par défaut) : le meilleur compromis général, équilibre naturellement exploration et exploitation
• “LCB” (Lower Confidence Bound) : plus exploratoire, utile quand on craint de rester coincé dans un optimum local
• “PI” (Probability of Improvement) : plus exploiteur, converge vite mais risque de manquer l’optimum global
• “EIps” (Expected Improvement per Second) : variante qui tient compte du temps d’évaluation

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

kernel

Le noyau du processus gaussien détermine la structure de corrélation supposée entre les points de l’espace de recherche :

• Matérn 5/2 (par défaut dans scikit-optimize) : bon compromis, permet des fonctions non infiniment dérivables
• RBF (Radial Basis Function) : suppose des fonctions très lisses, converge vite si l’hypothèse est vérifiée
• Matérn 3/2 : plus rugueux, adapté aux fonctions avec des variations abruptes

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

n_random_starts

Le nombre de points échantillonnés aléatoirement avant de commencer l’optimisation bayésienne proprement dite. Ces points initiaux permettent au processus gaussien d’avoir une première estimation raisonnable de la fonction :

• Valeur typique : 5 à 20 (environ 10 par défaut)
• Trop peu : le modèle initial peut être trompeur
• Trop : on gaspille des évaluations sans guide bayésien

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

Avantages et Limites

Avantages de l’Optimisation Bayésienne

1. Efficacité sample : trouve des solutions compétitives avec 5 à 10 fois moins d’évaluations que la recherche aléatoire. C’est l’avantage principal et le plus décisif.
2. Prise en compte de l’historique : chaque évaluation informe les suivantes. Contrairement au grid search ou au random search qui sont des méthodes sans mémoire, l’Optimisation Bayésienne apprend de chaque essai et affine progressivement sa stratégie.
3. Gestion naturelle du compromis exploration-exploitation : la fonction d’acquisition équilibre automatiquement l’exploration de nouvelles régions et l’exploitation des zones prometteuses, sans intervention manuelle.
4. Estimation de l’incertitude : le processus gaussien fournit non seulement une prédiction mais aussi un intervalle de confiance, ce qui permet de quantifier notre certitude sur chaque région de l’espace de recherche.
5. Compatible avec les fonctions bruyantes : le paramètre de bruit du processus gaussien permet de gérer les évaluations non déterministes, fréquentes en apprentissage automatique.

Limites de l’Optimisation Bayésienne

1. Coût computationnel du processus gaussien : l’inversion de la matrice de covariance K a une complexité en O(n³) où n est le nombre d’observations. Cela devient prohibitif au-delà de quelques milliers de points.
2. Malédiction de la dimensionnalité : l’Optimisation Bayésienne fonctionne bien pour des espaces de 1 à ~20 dimensions. Au-delà, le processus gaussien peine à modéliser efficacement la fonction et les méthodes basées sur les forêts d’arbres (SMAC, TPE) deviennent préférables.
3. Hypothèse de lissage : le processus gaussien suppose que la fonction est relativement lisse. Si la fonction objectif présente des discontinuités ou des variations très abruptes, le modèle peut être trompeur.
4. Sensibilité aux hyperparamètres du noyau : le choix et l’optimisation des hyperparamètres du noyau influencent significativement les résultats. Un mauvais noyau peut conduire à une convergence prématurée vers un optimum local.

4 Cas d’Usage Concrets

Cas d’usage 1 : Réglage d’hyperparamètres de réseaux de neurones profonds

L’entraînement d’un réseau de neurones profond peut prendre des heures, voire des jours. Chaque évaluation d’une configuration d’hyperparamètres (taux d’apprentissage, nombre de couches, nombre de neurones, taux de dropout, etc.) représente un investissement temps considérable. L’Optimisation Bayésienne réduit drastiquement le nombre d’entraînements nécessaires pour trouver une configuration performante.

# Architecture typique pour un réseau de neurones
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')
]

Cas d’usage 2 : Calibration de modèles physiques et simulations numériques

En ingénierie, la calibration de modèles physiques (mécanique des fluides, thermodynamique, science des matériaux) nécessite des simulations numériques extrêmement coûteuses. L’Optimisation Bayésienne permet de trouver les paramètres physiques qui minimisent l’écart entre simulation et mesure expérimentale avec un nombre minimal de simulations. Des entreprises comme Airbus et BMW utilisent cette approche pour optimiser la forme de leurs composants aérodynamiques.

Cas d’usage 3 : Optimisation de pipelines de traitement de données

Les pipelines de données en production comportent de nombreux paramètres : seuils de filtrage, tailles de fenêtres glissantes, taux d’échantillonnage, paramètres de prétraitement. L’Optimisation Bayésienne permet d’optimiser l’ensemble du pipeline de bout en bout en maximisant une métrique business (taux de conversion, revenu par utilisateur, précision de détection).

Cas d’usage 4 : Conception de molécules et découverte de médicaments

En chimie computationnelle et en découverte de médicaments, chaque évaluation correspond à une simulation de dynamique moléculaire ou à une expérience de laboratoire coûteuse. L’Optimisation Bayésienne guide la recherche vers les molécules les plus prometteuses en maximisant des propriétés cibles (affinité de liaison, solubilité, toxicité réduite) tout en minimisant le nombre d’expériences nécessaires. Des plateformes comme la découverte de médicaments assistée par machine learning utilisent cette approche quotidiennement.

Voir aussi

• Comment implémenter une recherche binaire en Python : Un guide étape par étape
• Dynamiser vos Projets Python : Simuler des Amibes dans une Grille 3D avec Python