Comment optimiser une fonction dont chaque évaluation prend six heures de simulation ou détruit une éprouvette ? La descente de gradient est inutilisable — pas de gradient, et un budget de quelques dizaines d'essais. C'est le domaine de l'optimisation bayésienne, qui construit un modèle probabiliste de la fonction et s'en sert pour décider où essayer ensuite. Pour un profil d'ingénieur, c'est probablement le module dont l'application est la plus immédiate.

1. Le problème : une boîte noire coûteuse

On cherche \(x^\star = \arg\max_{x\in\mathcal{X}} f(x)\) où \(f\) est :

• coûteuse — une évaluation demande des heures de calcul, un essai physique, ou une campagne d'entraînement complète ;
• opaque — pas d'expression analytique, donc pas de gradient exploitable ;
• possiblement bruitée — deux mesures au même point ne donnent pas le même résultat.

Le budget se compte en dizaines d'évaluations, pas en millions. Une recherche par grille gaspille tout sur des zones sans intérêt ; une descente de gradient numérique dépenserait tout son budget à estimer une seule direction (il faut \(d+1\) évaluations pour un seul gradient en dimension \(d\)).

La stratégie gagnante consiste à modéliser \(f\) à partir des points déjà évalués, puis à choisir le prochain point en exploitant ce modèle. Ce modèle de substitution s'appelle un surrogate, et le choix naturel est un processus gaussien parce qu'il fournit, en plus d'une prédiction, une incertitude — exactement l'information nécessaire pour savoir où l'on ne sait pas encore.

2. Processus gaussiens : une loi sur des fonctions

L'idée demande un petit effort d'abstraction. Une loi normale ordinaire tire des nombres. Une loi normale multivariée tire des vecteurs, avec des corrélations entre coordonnées. Un processus gaussien tire des fonctions : on peut le voir comme une loi normale multivariée de dimension infinie, dont chaque « coordonnée » est la valeur de la fonction en un point.

La définition opérationnelle est plus simple : quelle que soit la collection finie de points \(x_1,\dots,x_n\) que l'on choisit, le vecteur \(\big(f(x_1),\dots,f(x_n)\big)\) suit une loi normale multivariée. Le processus est entièrement déterminé par une fonction moyenne (souvent prise nulle après centrage) et une fonction de covariance, le noyau :

Le noyau dit à quel point deux points se « ressemblent » : \(k(x,x')\) grand signifie que \(f(x)\) et \(f(x')\) sont fortement corrélées, donc proches. C'est la seule hypothèse qu'on injecte, et elle suffit.

  A PRIORI : une nuée de fonctions plausibles      APRÈS 3 OBSERVATIONS
  ╱‾╲    ╱╲                                        ╱‾╲    ╱╲
 ╱   ╲╱╲╱  ╲   toutes passent "n'importe où",    ●╱   ╲╱‾╲   ← toutes passent par
╱  ╱╲ ╲  ╱╲ ╲  l'incertitude est uniforme         ╲╱      ╲●    les points observés
 ╲╱  ╲ ╲╱   ╲                                      ╱╲    ╱  ╲   et l'incertitude se
                                                  ●    ╲╱    ╲  resserre AUTOUR d'eux
  ░░░░░░░░░░░░ bande d'incertitude large       ░░░▒▒░░░▒▒▒░░░░  et reste large LOIN

3. Ce que le conditionnement donne

Conditionner une loi normale par une partie de ses coordonnées donne encore une loi normale, dont moyenne et covariance sont explicites. Avec \(K = k(X,X) + \sigma_n^2 I\) la matrice de covariance des points observés et \(k_* = k(X,x_*)\) le vecteur de covariance entre les observations et le point à prédire :

Deux lectures importantes. D'abord, la moyenne prédictive est une combinaison linéaire des observations pondérée par la similarité — c'est exactement une interpolation par fonctions de base radiales, le sujet de mon projet mdisd, mais assortie d'une variance. Ensuite, et c'est la propriété décisive : la variance prédictive ne dépend pas des valeurs observées, seulement de la position des points. L'incertitude est donc grande là où l'on n'a pas mesuré et petite près des mesures, indépendamment de ce qu'on y a trouvé. C'est précisément ce dont l'optimisation bayésienne a besoin pour décider où explorer.

4. Noyaux et hyperparamètres

RBF / gaussien      k = σ²·exp(-‖x-x'‖²/2ℓ²)      fonctions infiniment dérivables :
                                                   très (souvent trop) lisse.
Matérn ν=5/2        polynôme·exp(-√5 r/ℓ)          deux fois dérivable. LE DÉFAUT en
                                                   optimisation bayésienne : réaliste
                                                   pour un système physique.
Matérn ν=3/2                                       une fois dérivable : plus rugueux.
Périodique          exp(-2sin²(π r/p)/ℓ²)          saisonnalité de période connue.
Linéaire            x·x'                           tendance affine.
Blanc               σ²·δ(x,x')                     bruit d'observation.
Somme / produit     k₁+k₂, k₁·k₂                   combiner : tendance + saisonnalité
                                                   + bruit, chacun avec ses paramètres.

La longueur caractéristique \(\ell\) est l'hyperparamètre décisif : elle dit à quelle distance deux points cessent de s'informer mutuellement. Petite, le GP devient nerveux et l'incertitude explose entre les points ; grande, il lisse tout et manque les variations réelles. Avec une longueur par dimension (Automatic Relevance Determination), le modèle apprend lui-même quelles variables comptent — une grande longueur signifie que la sortie est insensible à cette entrée, ce qui constitue au passage une analyse de sensibilité gratuite (module 21).

Les hyperparamètres (\(\ell\), \(\sigma^2\), \(\sigma_n^2\)) s'ajustent en maximisant la vraisemblance marginale :

Le premier terme récompense l'ajustement aux données, le second pénalise la complexité du modèle — une matrice de covariance « permissive », qui autorise beaucoup de fonctions, a un grand déterminant. Le rasoir d'Occam est ici automatique, sans validation croisée : c'est l'un des grands attraits du cadre bayésien (module 10).

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, ConstantKernel, WhiteKernel

noyau = (ConstantKernel(1.0)                          # amplitude σ²
         * Matern(length_scale=np.ones(3), nu=2.5)    # ARD : une longueur par dimension
         + WhiteKernel(noise_level=1e-3))             # bruit d'observation

gp = GaussianProcessRegressor(kernel=noyau, normalize_y=True, n_restarts_optimizer=10)
gp.fit(X, y)                       # optimise les hyperparamètres par vraisemblance

mu, sigma = gp.predict(X_test, return_std=True)
print(gp.kernel_)                  # les longueurs apprises = sensibilité par variable
# une longueur très grande sur une dimension → cette variable n'influence pas la sortie

Le coût est la limite structurelle : la factorisation de \(K\) est en \(O(n^3)\) et la mémoire en \(O(n^2)\). Au-delà de quelques milliers de points, il faut des approximations creuses (points inducteurs, variational GP dans GPyTorch). Pour l'optimisation bayésienne, où \(n\) dépasse rarement quelques centaines, ce n'est pas un obstacle.

5. La boucle d'optimisation bayésienne

1. Évaluer f sur un plan initial (Latin Hypercube ou Sobol, ~10 points par dimension)
2. Répéter jusqu'à épuisement du budget :
      a. ajuster le GP sur TOUTES les observations disponibles
      b. maximiser la fonction d'acquisition α(x) sur le domaine
         (α est bon marché : on peut y lancer un optimiseur classique multi-départs)
      c. évaluer f au point trouvé — LA seule opération coûteuse
      d. ajouter (x, y) aux observations
3. Renvoyer le meilleur point observé

Le déplacement conceptuel est là : on remplace l'optimisation d'une fonction coûteuse par une suite d'optimisations d'une fonction bon marché, construite à partir du GP. Consommer quelques secondes de calcul pour choisir intelligemment un point vaut largement le coup quand ce point coûte six heures de simulation.

6. Fonctions d'acquisition

L'acquisition traduit le compromis exploitation (aller là où le modèle prédit une bonne valeur) / exploration (aller là où il ne sait rien). En notant \(f^+\) le meilleur résultat observé, \(Z = (\mu(x)-f^+-\xi)/\sigma(x)\), \(\Phi\) et \(\varphi\) les fonctions de répartition et de densité normales :

L'amélioration espérée (EI) est l'espérance du gain par rapport au meilleur point connu, sachant que l'on ne perd rien si le nouveau point est moins bon. Sa structure est parlante : le premier terme récompense une moyenne élevée, le second une forte incertitude. C'est le choix par défaut — pas de paramètre à régler (\(\xi\) vaut typiquement 0,01) et un bon équilibre empirique.

      μ(x) prédite        ╱‾‾╲                    ●  points déjà évalués
   ●            ╱‾╲     ╱     ╲       ●          ░  incertitude σ(x)
  ░░░░       ░░╱   ╲░░░╱       ╲░░░░░░░
             ▲              ▲
             │              └── EXPLORATION : zone mal connue, σ grand
             └── EXPLOITATION : moyenne élevée près d'un bon point
  EI équilibre les deux ; UCB laisse régler le curseur par β.

UCB       μ(x) + β·σ(x)     exploration réglable explicitement. Garanties théoriques
                            de regret (GP-UCB). β ≈ 2 en pratique.
PI        P(f(x) > f⁺)      probabilité d'amélioration : trop exploitante, converge
                            souvent vers un optimum local.
Thompson  tirer une fonction du postérieur et l'optimiser. Naturellement
          parallélisable : k tirages donnent k points distincts.
Entropy search  maximise l'information gagnée SUR LA POSITION de l'optimum plutôt
                que sur sa valeur. Plus coûteux, souvent plus efficace.
Knowledge gradient  anticipe le gain sur la meilleure valeur atteignable après
                    l'évaluation.

7. BoTorch en pratique

BoTorch (Meta) est construit sur PyTorch et GPyTorch. Son atout : les acquisitions y sont différentiables et estimées par Monte-Carlo, ce qui permet de les maximiser par gradient et d'en construire des variantes (batch, contraintes, multi-objectif) sans dériver de nouvelles formules analytiques. Ax est la surcouche qui gère la campagne complète.

import torch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.acquisition import qLogExpectedImprovement
from botorch.optim import optimize_acqf
from botorch.utils.sampling import draw_sobol_samples
from gpytorch.mlls import ExactMarginalLogLikelihood

bornes = torch.tensor([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]], dtype=torch.double)

X = draw_sobol_samples(bounds=bornes, n=10, q=1).squeeze(1)   # plan initial
Y = simulation_couteuse(X).unsqueeze(-1)                      # à MAXIMISER

for iteration in range(40):
    gp  = SingleTaskGP(X, Y)                       # normalisation gérée en interne
    mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
    fit_gpytorch_mll(mll)                          # hyperparamètres par vraisemblance

    acq = qLogExpectedImprovement(model=gp, best_f=Y.max())
    candidat, _ = optimize_acqf(
        acq, bounds=bornes,
        q=4,                                       # 4 points EN PARALLÈLE
        num_restarts=20, raw_samples=512,          # multi-départs sur l'acquisition
    )

    X = torch.cat([X, candidat])
    Y = torch.cat([Y, simulation_couteuse(candidat).unsqueeze(-1)])

print("meilleur point :", X[Y.argmax()], "valeur :", Y.max().item())

Le paramètre q=4 illustre un cas fréquent en ingénierie : si quatre simulations peuvent tourner en parallèle sur un cluster, on veut quatre points complémentaires, pas quatre fois le même. Les acquisitions « q » de BoTorch traitent cela correctement, en pénalisant la redondance à l'intérieur du lot.

8. Contraintes, multi-objectif, multi-fidélité

• Contraintes coûteuses — quand la faisabilité elle-même doit être simulée (une contrainte mécanique, une température maximale), on modélise chaque contrainte par son propre GP et l'on pondère l'acquisition par la probabilité de faisabilité. Le point proposé est alors celui qui maximise le gain attendu et réalisable.
• Multi-objectif — rendement et coût, précision et latence. On ne cherche plus un optimum mais un front de Pareto : l'ensemble des compromis qu'aucun autre point ne domine sur tous les critères à la fois. L'acquisition de référence est qEHVI (expected hypervolume improvement), qui mesure le gain de volume dominé.
• Espaces mixtes — variables continues, entières et catégorielles mélangées (un choix de matériau, un nombre de couches). Traitées par des noyaux dédiés ou par relaxation continue.
• Multi-fidélité — souvent le gain le plus spectaculaire : on dispose d'un simulateur grossier rapide et d'un simulateur fin coûteux. Un GP multi-tâche apprend la corrélation entre les deux niveaux et n'appelle le coûteux que lorsque c'est décisif. Le même principe s'applique au réglage d'hyperparamètres, où « peu d'époques » est une fidélité basse — c'est l'idée de Hyperband, croisée avec la modélisation bayésienne dans BOHB.

Méthode              Évaluations nécessaires   Parallélisme   Quand l'utiliser
──────────────────────────────────────────────────────────────────────────────────
Grid search          très nombreuses           parfait        2-3 paramètres, pas plus
Random search        nombreuses                parfait        baseline honnête, robuste
Hyperband / ASHA     moyennes                  bon            entraînements interruptibles
TPE (Optuna)         peu                       bon            défaut pratique en ML
GP + EI (BoTorch)    très peu                  bon (batch)    évaluations RÉELLEMENT
                                                              coûteuses, dimension < 20

Un repère honnête : en dessous de 20 ou 30 évaluations disponibles et en dimension modérée, l'optimisation bayésienne par GP domine nettement. Pour régler les hyperparamètres d'un gradient boosting qui s'entraîne en dix secondes, Optuna (TPE) est plus simple et suffisant (module 14). Et la recherche aléatoire reste une baseline étonnamment solide, qu'il faut toujours inclure dans la comparaison.

9. Cas d'usage et limites

Plans d'expériences      choisir les essais physiques suivants (chimie, matériaux,
                         procédés) — chaque essai coûte des jours et des matières
Calibration de modèle    identifier les paramètres d'un modèle physique à partir
                         de mesures (module 21)
Conception assistée      optimiser une géométrie évaluée par CFD ou éléments finis
Réglage coûteux          hyperparamètres d'un gros réseau (des heures par essai)
Contrôle / robotique     réglage de contrôleurs sur système réel, où chaque essai
                         a un coût et un risque
A/B testing adaptatif    allocation séquentielle, à la frontière des bandits (module 24)

What this module covers

• The expensive black-box optimisation problem, and why gradient methods do not apply.
• Gaussian processes as distributions over functions: the posterior mean and variance derived from conditioning, why the predictive variance ignores observed values, and the Cholesky-based implementation.
• Kernels (RBF, Matérn, periodic, composite), length scales and ARD as free sensitivity analysis, and hyperparameter fitting by marginal likelihood with its built-in Occam penalty.
• The Bayesian optimisation loop: replacing one expensive optimisation with many cheap ones.
• Acquisition functions: expected improvement derived and explained, UCB, probability of improvement, Thompson sampling, entropy search and knowledge gradient.
• A complete BoTorch loop, including batch acquisition for parallel simulations.
• Extensions: expensive constraints, multi-objective optimisation and Pareto fronts (qEHVI), mixed spaces, and multi-fidelity modelling.
• A comparison of tuning strategies (grid, random, Hyperband, TPE, GP-based) and an honest account of where Bayesian optimisation breaks down.

Course outline

1. Modern Python: typing, decorators, iterators and tooling
2. Concurrency in Python: the GIL, threads, processes, asyncio
3. Numerical toolbox: NumPy, pandas, polars, Dask
4. SQL and databases
5. Docker, Compose and Kubernetes
6. Serving a model: FastAPI, Pydantic, Jinja
7. Data engineering: formats, pipelines and streams
8. Math for ML: optimisation, linear algebra, SVD and PCA
9. Inferential statistics: estimate, test, decide
10. Probability, Bayesian inference and uncertainty quantification
11. Causal inference: measuring the effect of an action
12. Data science in practice: from problem to dataset
13. Learning with few labels
14. Classical machine learning
15. Learning theory: why generalisation happens
16. Deep learning
17. Time series, signals and time-to-event data
18. Non-tabular data: text, images, audio, graphs
19. Gaussian processes and Bayesian optimisation
20. Reinforcement learning
21. Machine learning for simulation and engineering
22. LLMs, embeddings, RAG and adaptation
23. Search, ranking and recommendation
24. From scores to decisions
25. Scaling up: GPUs, distributed training, performance
26. MLOps: taking a model to production
27. Fairness, robustness, privacy and security
28. End-to-end project: property listings, from scraping to price prediction