Probabilités, statistiques et transport optimal

Objectifs

Un modèle génératif ne produit pas « un » champ magnétique : il produit une loi de probabilité sur des champs, que l’on observe à travers des échantillons. Toutes les questions du stage MilaDy, qu’il s’agisse d’apprendre cette loi, de la contraindre ou de l’évaluer, s’écrivent donc dans la langue des probabilités. Ce module en donne les outils, du niveau L2 jusqu’aux notions qui fondent le flow matching et les contraintes « en loi ».

À la fin des vingt-cinq heures de travail prévues, vous saurez :

1. manipuler variables et vecteurs aléatoires, densités, espérance, variance et covariance, et calculer avec la loi gaussienne multivariée (échantillonnage par Cholesky, lois marginales et conditionnelles) ;
2. calculer la loi image d’une variable aléatoire par une application inversible grâce à la formule du jacobien, et expliquer pourquoi cette formule est la base des flots normalisants ;
3. estimer une espérance par Monte-Carlo avec un intervalle de confiance, justifier l’erreur en 1/\sqrt{N}, utiliser l’échantillonnage d’importance autonormalisé et en diagnostiquer la qualité par la taille d’échantillon effective (ESS) ;
4. démontrer la positivité de la divergence de Kullback–Leibler, obtenir les familles exponentielles par le principe du maximum d’entropie, démontrer par dualité que la loi la plus proche d’une loi q sous une contrainte de moments est une inclinaison exponentielle de q, et la calculer numériquement par la méthode de Newton ;
5. calculer les distances de Wasserstein W_1 et W_2 en dimension 1 par les quantiles, W_2 entre deux gaussiennes, la distance d’énergie et la MMD sur des échantillons, et choisir une métrique adaptée pour comparer des champs générés ;
6. simuler une équation différentielle stochastique par le schéma d’Euler–Maruyama, énoncer la formule d’Itô, établir l’équation de Fokker–Planck et sa solution stationnaire pour le processus d’Ornstein–Uhlenbeck ;
7. démontrer qu’un champ de vitesse transporte une loi selon l’équation de continuité \partial_t p + \nabla\cdot(p\,v) = 0, et en déduire pourquoi la vitesse « marginale » du flow matching transporte le bruit vers les données.

Prérequis

Le point de départ est le niveau de fin de deuxième année de licence (L2) scientifique, mais un étudiant de L1 motivé peut suivre les chapitres 1 à 3 en revoyant au passage les intégrales multiples. Il faut savoir calculer une intégrale simple et une intégrale double, effectuer un changement de variables dans une intégrale double (jacobien), manipuler des matrices symétriques (valeurs propres, déterminant, inverse) et des séries. Une première rencontre avec les probabilités discrètes (lancers de dés, loi binomiale) aide, sans être indispensable : le chapitre 1 reprend tout depuis les définitions.

Côté programmation, il suffit de savoir manipuler des tableaux NumPy et écrire une fonction en Python. SciPy sert pour quelques fonctions spéciales et PyTorch uniquement pour la différentiation automatique au chapitre 2.

Le module Calcul vectoriel et analyse de Fourier est utile pour le chapitre 6 (divergence, théorème de Green–Ostrogradski) et pour la section 5.6 (spectres). Le module Systèmes dynamiques, bifurcations et mesures invariantes peut être lu en parallèle : la notion de mesure invariante y est introduite du point de vue des trajectoires, alors qu’elle apparaît ici, au chapitre 6, du point de vue de l’équation de Fokker–Planck. Ce module prépare directement Apprentissage profond pour la physique, Modèles génératifs : diffusion et flow matching et Modèles génératifs sous contraintes physiques.

Plan de travail

Le niveau indique la difficulté : 1 pour les fondations, 2 pour le cœur du module, 3 pour ce qui touche déjà aux questions de recherche du stage. Il ne redescend jamais d’une séance à la suivante.

Chaque séance se termine par les exercices indiqués, à faire avant de lire leur solution. Les blocs de code sont courts (moins de 30 secondes chacun sur un seul cœur de processeur) : recopiez-les, exécutez-les, puis modifiez-les. Les schémas interactifs se manipulent directement dans la page ; ils reprennent les exemples du texte avec des curseurs.

1. Variables aléatoires, lois et moments

Un ensemble de champs MHD simulés à Rm donné n’est pas une collection de fichiers sans structure : c’est un échantillon d’une loi \mu_{\boldsymbol\lambda}, et ce que l’on demande à un modèle génératif, c’est de reproduire cette loi. Ce chapitre fixe le vocabulaire (variable aléatoire, loi, densité, moments) et se termine par la loi gaussienne multivariée, point de départ de tous les modèles de diffusion et de flow matching.

1.1 Espace probabilisé, variable aléatoire, loi

On note l’univers \mathcal{U} plutôt que \Omega, la lettre habituelle, pour ne pas le confondre avec le domaine physique \Omega \subset \mathbb{R}^d des autres modules. En pratique, on ne manipule presque jamais \mathcal{U} directement : on manipule des variables aléatoires et leurs lois.

La définition 1.2 dit déjà quelque chose d’important pour la suite : la loi de X est l’image de la probabilité \mathbb{P} par l’application X. Le chapitre 2 généralise cette opération « mesure image » à toute application, et c’est elle qu’un modèle génératif réalise.

1.2 Espérance, variance et deux inégalités

La formule de transfert est le principe de tous les calculs de ce module : pour calculer la moyenne d’une observable g(X), on n’a pas besoin de connaître la loi de g(X), seulement celle de X.

L’égalité \mathbb{E}[(X - \mathbb{E}X)^2] = \mathbb{E}[X^2] - \mathbb{E}[X]^2 s’obtient en développant le carré et en utilisant la linéarité. Elle montre au passage, avec Jensen appliqué à \varphi(x) = x^2, que \mathbb{E}[X]^2 \le \mathbb{E}[X^2].

L’inégalité de Tchebychev est grossière, mais elle suffit pour démontrer la loi des grands nombres au chapitre 3.

1.3 Vecteurs aléatoires, covariance, espérance conditionnelle

Un champ discrétisé \mathbf{x} \in \mathbb{R}^D est un vecteur aléatoire : chaque point de grille est une coordonnée, et toute l’information physique (corrélations spatiales, spectre, hélicité) est dans la loi jointe des coordonnées, pas dans les lois des coordonnées prises une à une.

Le point (iii) est un piège classique quand on valide un modèle génératif : reproduire les covariances des données ne garantit pas de reproduire leurs dépendances non linéaires, comme celles qui lient l’amplitude d’un mode à son signe dans un régime de renversements.

Ce théorème, élémentaire, est le ressort mathématique du flow matching : on y entraîne un réseau à prédire par moindres carrés une quantité aléatoire, et l’optimum est une espérance conditionnelle (proposition 6.15).

1.4 La loi gaussienne multivariée

Cette proposition est utilisée dans tous les modèles génératifs du parcours : on tire \mathbf{Z} coordonnée par coordonnée, puis on transforme. La définition 1.14 a aussi une conséquence propre aux gaussiennes, qui fait exception à la mise en garde du point (iii) de la proposition 1.11 : pour un vecteur gaussien, des coordonnées non corrélées sont indépendantes, car une matrice \Sigma diagonale rend la densité factorisable. L’hypothèse porte sur la loi jointe : deux variables qui sont chacune gaussiennes, sans former ensemble un vecteur gaussien, peuvent être non corrélées et dépendantes.

La moyenne conditionnelle est affine en \mathbf{x}_2 et la variance conditionnelle ne dépend pas de \mathbf{x}_2 : c’est la régression linéaire, et c’est aussi le calcul qui donne la vitesse du flow matching entre gaussiennes (exercice 15).

import numpy as np

rng = np.random.default_rng(0)
m = np.array([1.0, -2.0])
S = np.array([[2.0, 1.2],
              [1.2, 1.0]])
L = np.linalg.cholesky(S)                 # S = L L^T, L triangulaire inférieure
N = 400_000
Z = rng.standard_normal((N, 2))           # Z ~ N(0, I)
X = m + Z @ L.T                           # X = m + L Z ~ N(m, S)
print("moyenne empirique    :", X.mean(axis=0).round(3))
print("covariance empirique :", np.cov(X.T).round(3).tolist())

# Loi conditionnelle de X1 sachant X2 = -1 (proposition 1.16)
x2 = -1.0
mu_c = m[0] + S[0, 1] / S[1, 1] * (x2 - m[1])
var_c = S[0, 0] - S[0, 1] ** 2 / S[1, 1]
sel = np.abs(X[:, 1] - x2) < 0.01         # tranche fine autour de X2 = -1
print(f"théorie   : E[X1 | X2=-1] = {mu_c:.3f}, Var = {var_c:.3f}")
print(f"empirique : E[X1 | X2=-1] = {X[sel, 0].mean():.3f}, Var = {X[sel, 0].var():.3f}  ({sel.sum()} points)")

# Covariance nulle n'entraîne pas l'indépendance
Y = rng.standard_normal(N)
print("Cov(Y, Y^2) =", round(np.cov(Y, Y ** 2)[0, 1], 3),
      "; Corr(|Y|, Y^2) =", round(np.corrcoef(np.abs(Y), Y ** 2)[0, 1], 3))
# sortie :
# moyenne empirique    : [ 1.003 -1.998]
# covariance empirique : [[2.002, 1.199], [1.199, 0.999]]
# théorie   : E[X1 | X2=-1] = 2.200, Var = 0.560
# empirique : E[X1 | X2=-1] = 2.205, Var = 0.567  (1959 points)
# Cov(Y, Y^2) = 0.009 ; Corr(|Y|, Y^2) = 0.936

Les écarts entre théorie et empirique sont de l’ordre de 10^{-2} pour la tranche conditionnelle, qui ne contient que 1959 points : le chapitre 3 explique pourquoi cet ordre de grandeur est celui attendu. La covariance de Y et Y^2 est nulle aux fluctuations près, alors que |Y| et Y^2 sont corrélés à 94 %.

2. Changement de variables et mesure image

Un modèle génératif est, dans sa forme la plus simple, une application T_\theta qui transforme un bruit gaussien \mathbf{Z} en un champ T_\theta(\mathbf{Z}). La loi générée est donc la mesure image de \mathcal{N}(0, I) par T_\theta. Ce chapitre définit cette opération, donne la formule du jacobien qui calcule la densité image, et montre comment elle fonde les flots normalisants ; le chapitre 6 en donnera la version en temps continu, qui fonde le flow matching.

2.1 Mesure image

Cette formule est la raison pour laquelle on peut évaluer un modèle génératif sans jamais écrire sa densité : pour calculer la moyenne d’une observable (énergie, hélicité, divergence) sous la loi générée, il suffit d’évaluer l’observable sur des échantillons T_\theta(\mathbf{Z}_i).

En dimension 1, on sait donc transporter n’importe quelle loi sans atome sur n’importe quelle autre par une application croissante : si F est continue, F(X) est uniforme et G^{-1}\circ F envoie la loi de fonction de répartition F sur celle de fonction de répartition G. Cette application « des quantiles » réapparaîtra au chapitre 5 comme le transport optimal en dimension 1.

2.2 La formule du jacobien

L’interprétation est celle d’une conservation de la masse : la probabilité contenue dans un petit volume \mathrm{d}\mathbf{x} autour de \mathbf{x} se retrouve dans le volume |\det J_T|\,\mathrm{d}\mathbf{x} autour de T(\mathbf{x}). Si T dilate localement les volumes, la densité diminue d’autant. La figure 1 le montre sur un exemple en dimension 1.

Le schéma interactif 1 permet de faire varier une application monotone en dimension 1 et de voir les tranches de probabilité se déformer.

2.3 Flots normalisants

La section 4.2 montrera que maximiser la vraisemblance revient à minimiser la divergence de Kullback–Leibler \mathrm{KL}(\mu\,\|\,q_\theta). La revue de Papamakarios et al. (2021) décrit les architectures qui rendent le déterminant bon marché (couplages affines, transformations autorégressives, jacobiennes triangulaires).

Le code suivant vérifie la formule (2.1) en dimension 1, compare le log-déterminant du flot planaire à celui obtenu par différentiation automatique, puis contrôle la densité image en dimension 2 de deux façons indépendantes : par intégration de p_{\mathbf{Y}} sur une grille de \mathbf{y} (ce qui exige T^{-1}) et par comptage d’échantillons.

import numpy as np
import torch
torch.set_num_threads(1)
rng = np.random.default_rng(1)

# 1) Dimension 1 : Y = sinh(X) avec X ~ N(0, 1)
phi = lambda x: np.exp(-x ** 2 / 2) / np.sqrt(2 * np.pi)
p_Y = lambda y: phi(np.arcsinh(y)) / np.sqrt(1 + y ** 2)   # p_X(T^{-1}(y)) |(T^{-1})'(y)|
Y = np.sinh(rng.standard_normal(1_000_000))
hist, edges = np.histogram(Y, bins=80, range=(-4, 4))
centres = (edges[:-1] + edges[1:]) / 2
dens = hist / (len(Y) * (edges[1] - edges[0]))              # histogramme normalisé en densité
print("écart max histogramme / formule :", f"{np.abs(dens - p_Y(centres)).max():.4f}",
      "; p_Y(0) =", f"{p_Y(0.0):.4f}")

# 2) Dimension 2 : flot « planaire » T(x) = x + u tanh(w.x + b)
u = torch.tensor([1.5, -0.5]); w = torch.tensor([2.0, 1.0]); b = 0.3
def T(x):
    return x + u * torch.tanh(x @ w + b).unsqueeze(-1)
def logdet(x):                      # det(I + u psi^T) = 1 + u.psi avec psi = (1 - tanh^2) w
    return torch.log(torch.abs(1 + (1 - torch.tanh(x @ w + b) ** 2) * (u @ w)))
x = torch.tensor([0.4, -1.1])
J = torch.autograd.functional.jacobian(T, x)                # jacobienne 2 x 2 par autograd
print("log|det J| : autograd =", f"{torch.logdet(J).item():.6f}", "; formule =", f"{logdet(x).item():.6f}")

def T_inv(y, it=60):
    # a = w.x + b vérifie w.y + b = a + (u.w) tanh(a) : équation scalaire croissante (u.w > -1)
    c = y @ w + b
    a = c.clone()
    for _ in range(it):                                     # Newton scalaire
        a = a - (a + (u @ w) * torch.tanh(a) - c) / (1 + (u @ w) * (1 - torch.tanh(a) ** 2))
    return y - u * torch.tanh(a).unsqueeze(-1)

# 3) Vérification : p_Y(y) = p_X(T^{-1}(y)) / |det J_T(T^{-1}(y))|
g = torch.linspace(-8, 8, 1601, dtype=torch.float64); h = (g[1] - g[0]).item()
u, w = u.double(), w.double()
GY = torch.stack(torch.meshgrid(g, g, indexing="ij"), -1).reshape(-1, 2)
GX = T_inv(GY)
print("erreur d'inversion max :", f"{(T(GX) - GY).abs().max().item():.1e}")
pY = torch.exp(-(GX ** 2).sum(-1) / 2 - logdet(GX)) / (2 * np.pi)
print("intégrale de p_Y sur la grille :", f"{pY.sum().item() * h * h:.4f}")
region = (GY[:, 0] > 0.5) & (GY[:, 1] < 0.0)
Xs = torch.randn(400_000, 2, generator=torch.Generator().manual_seed(0), dtype=torch.float64)
Ys = T(Xs)
print("P(Y1 > 0.5, Y2 < 0) : formule =", f"{(pY * region).sum().item() * h * h:.4f}",
      "; échantillons =", f"{((Ys[:, 0] > 0.5) & (Ys[:, 1] < 0)).double().mean().item():.4f}")
# sortie :
# écart max histogramme / formule : 0.0050 ; p_Y(0) = 0.3989
# log|det J| : autograd = 1.252763 ; formule = 1.252763
# erreur d'inversion max : 2.7e-15
# intégrale de p_Y sur la grille : 1.0000
# P(Y1 > 0.5, Y2 < 0) : formule = 0.2799 ; échantillons = 0.2794

L’écart de 0{,}005 entre histogramme et formule en dimension 1 est l’erreur statistique attendue : les classes centrales, de largeur 0{,}1, contiennent environ 4\times10^4 tirages chacune, si bien que la densité y est estimée à 0{,}002 près (un écart-type), et le maximum sur 80 classes atteint deux à trois écarts-types. En dimension 2, la densité image intègre à 1 et prédit la probabilité d’une région à 5\times10^{-4} près, ce qui est compatible avec l’erreur Monte-Carlo de 400\,000 tirages (environ 7\times10^{-4}, voir le chapitre 3).

3. Estimation Monte-Carlo

Toute métrique d’un modèle génératif (spectre moyen, hélicité moyenne, taux de dissipation, fréquence des renversements) est une espérance que l’on estime avec un nombre fini d’échantillons. Avant de conclure qu’un modèle « respecte le bilan d’énergie », il faut savoir quelle erreur on commet du seul fait de l’échantillonnage. Ce chapitre donne cette barre d’erreur, puis montre comment réutiliser des échantillons d’une loi pour estimer des espérances sous une autre loi, ce qui servira à l’inclinaison exponentielle du chapitre 4.

3.1 L’estimateur Monte-Carlo et la loi des grands nombres

La loi forte des grands nombres, admise ici, affirme davantage : \hat I_N \to I presque sûrement dès que \mathbb{E}_\mu|f| < \infty (Durrett, 2019). Elle justifie qu’une seule longue suite de tirages suffise.

3.2 Théorème central limite et intervalle de confiance

3.3 Échantillonnage d’importance

Les événements rares, comme un renversement de polarité dans un régime où ils sont exceptionnels, posent un problème : si p = \mathbb{P}(A) = 3\times10^{-5}, un Monte-Carlo naïf avec 10^5 tirages voit l’événement trois fois, et l’erreur relative est de l’ordre de 60 %. L’idée de l’échantillonnage d’importance est de tirer ailleurs, puis de corriger par des poids.

La proposition optimale est inutilisable telle quelle (elle contient l’inconnue I), mais elle indique la direction : tirer là où |f|\,p est grand. Pour \mathbb{P}(X > 4) avec X \sim \mathcal{N}(0, 1), la proposition q = \mathcal{N}(4, 1) place la moitié des tirages dans l’événement ; le poids vaut w(y) = e^{8 - 4y}, et le code ci-dessous obtient une erreur relative de 3\times10^{-3} avec 10^5 tirages.

Souvent, la densité cible n’est connue qu’à une constante près : p = \tilde p/Z avec Z inconnu. C’est le cas d’une loi inclinée q\,e^{\lambda g}/Z au chapitre 4.

L’interprétation usuelle est que la variance de \hat I_N^{\text{AN}} est à peu près celle d’une moyenne de \mathrm{ESS} tirages indépendants de p. C’est une heuristique, utile mais pas un théorème : Elvira, Martino et Robert (2022) en discutent les limites. L’exercice 7 montre que pour une cible \mathcal{N}(0, 1) et une proposition \mathcal{N}(a, 1), \mathrm{ESS}/N \to e^{-a^2} : décaler la proposition de 2 écarts-types ne laisse que 1,8 % de tirages utiles.

import numpy as np
from scipy.stats import norm
rng = np.random.default_rng(2)

# 1) Estimer I = E[cos X], X ~ N(0, 1) ; valeur exacte exp(-1/2)
I = np.exp(-0.5)
print("     N   estimation   demi-largeur IC 95 %   |erreur|   |erreur| * sqrt(N)")
for N in [10 ** 2, 10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6]:
    f = np.cos(rng.standard_normal(N))
    m, s = f.mean(), f.std(ddof=1)
    print(f"{N:>7d}   {m:.5f}      {1.96 * s / np.sqrt(N):.5f}             {abs(m - I):.5f}    {abs(m - I) * np.sqrt(N):.3f}")
# Couverture de l'intervalle de confiance sur 2000 répétitions avec N = 1000
f = np.cos(rng.standard_normal((2000, 1000)))
m, s = f.mean(1), f.std(1, ddof=1)
print("taux de couverture de l'IC 95 % :", np.mean(np.abs(m - I) <= 1.96 * s / np.sqrt(1000)))

# 2) Événement rare p = P(X > 4)
p = norm.sf(4)
N = 100_000
x = rng.standard_normal(N)
print(f"\np exact = {p:.3e} ; MC naïf : {np.mean(x > 4):.3e} ({np.sum(x > 4)} succès sur {N})")
# Échantillonnage d'importance : proposition q = N(4, 1), poids w = p(y)/q(y) = exp(8 - 4y)
y = rng.normal(4.0, 1.0, N)
w = np.exp(8 - 4 * y)
est = np.mean(w * (y > 4))
se = np.std(w * (y > 4), ddof=1) / np.sqrt(N)
print(f"IS : {est:.4e} +/- {1.96 * se:.1e} (IC 95 %) ; erreur relative {abs(est - p) / p:.1e}")

# 3) Poids autonormalisés et taille d'échantillon effective (ESS)
def ess(w):
    return w.sum() ** 2 / (w ** 2).sum()
for mu_q in [0.0, 1.0, 2.0, 3.0]:                     # proposition N(mu_q, 1), cible N(0, 1)
    y = rng.normal(mu_q, 1.0, N)
    logw = -y ** 2 / 2 + (y - mu_q) ** 2 / 2          # log p - log q
    w = np.exp(logw - logw.max())                     # stabilisation numérique
    m2 = np.sum(w * y ** 2) / np.sum(w)               # estimateur autonormalisé de E_p[X^2] = 1
    print(f"q = N({mu_q:.0f}, 1) : ESS/N = {ess(w) / N:.3g} (théorie exp(-mu^2) = {np.exp(-mu_q ** 2):.3g}) ; E_p[X^2] ~ {m2:.3f}")
# sortie :
#      N   estimation   demi-largeur IC 95 %   |erreur|   |erreur| * sqrt(N)
#     100   0.63497      0.08392             0.02844    0.284
#    1000   0.59778      0.02829             0.00875    0.277
#   10000   0.60326      0.00884             0.00328    0.328
#  100000   0.60740      0.00277             0.00087    0.275
# 1000000   0.60644      0.00088             0.00009    0.089
# taux de couverture de l'IC 95 % : 0.955
#
# p exact = 3.167e-05 ; MC naïf : 3.000e-05 (3 succès sur 100000)
# IS : 3.1588e-05 +/- 4.2e-07 (IC 95 %) ; erreur relative 2.6e-03
# q = N(0, 1) : ESS/N = 1 (théorie exp(-mu^2) = 1) ; E_p[X^2] ~ 1.002
# q = N(1, 1) : ESS/N = 0.367 (théorie exp(-mu^2) = 0.368) ; E_p[X^2] ~ 1.001
# q = N(2, 1) : ESS/N = 0.0158 (théorie exp(-mu^2) = 0.0183) ; E_p[X^2] ~ 1.114
# q = N(3, 1) : ESS/N = 0.00101 (théorie exp(-mu^2) = 0.000123) ; E_p[X^2] ~ 0.882

Trois enseignements. Le produit |\text{erreur}|\times\sqrt N reste de l’ordre de \sigma \approx 0{,}45 (pour une erreur gaussienne, sa moyenne vaut \sigma\sqrt{2/\pi} \approx 0{,}36) ; l’intervalle de confiance couvre la vraie valeur dans 95,5 % des répétitions. L’échantillonnage d’importance transforme un événement invisible en une estimation précise à 0{,}3 % près. Enfin, la dernière ligne contient un piège : pour une proposition décalée de a = 3 (variable mu_q du code), l’ESS estimé (10^{-3}) est huit fois plus grand que l’ESS limite (1{,}2\times10^{-4}). Les tirages n’ont pas encore visité la région qui porte les plus gros poids, si bien que l’ESS empirique est optimiste et que l’estimation 0{,}882 de \mathbb{E}_p[X^2] = 1 est fausse sans que rien ne l’annonce. Un ESS faible est un signal d’alarme fiable ; un ESS correct n’est pas une garantie.

4. Entropie, divergence de Kullback–Leibler et I-projection

Supposons qu’un modèle génératif entraîné produise une loi q presque juste, mais qui viole en moyenne le bilan d’énergie : l’injection moyenne ne compense pas la dissipation moyenne. Comment corriger q le moins possible pour qu’elle satisfasse cette contrainte statistique ? Ce chapitre répond : il faut mesurer l’écart entre lois par la divergence de Kullback–Leibler, et la loi corrigée est alors une inclinaison exponentielle de q. C’est le mécanisme de la contrainte « en loi » du bilan d’énergie testée dans l’expérience E2.

4.1 Entropie

4.2 Divergence de Kullback–Leibler

La divergence de Kullback–Leibler n’est pas une distance : elle n’est pas symétrique et ne vérifie pas l’inégalité triangulaire. L’exercice 8 donne la formule entre deux gaussiennes,

\mathrm{KL}\big(\mathcal{N}(m_1, s_1^2)\,\|\,\mathcal{N}(m_2, s_2^2)\big) = \log\frac{s_2}{s_1} + \frac{s_1^2 + (m_1 - m_2)^2}{2s_2^2} - \frac12,

dont l’asymétrie est visible dès que s_1 \neq s_2.

Cette invariance distingue la divergence de Kullback–Leibler des distances de Wasserstein du chapitre 5, qui dépendent de la géométrie de l’espace : changer d’unités ou de variables ne change pas la KL, mais change W_p.

4.3 Maximum d’entropie et familles exponentielles

Le principe du maximum d’entropie (Jaynes, 1957) dit : parmi toutes les lois compatibles avec ce que l’on sait (des valeurs moyennes), choisir la plus étalée, celle qui n’ajoute aucune information non justifiée.

4.4 I-projection : la loi la plus proche sous contrainte

On remplace maintenant la mesure de Lebesgue du maximum d’entropie par une loi de référence q, celle produite par le modèle génératif. On note les multiplicateurs \lambda \in \mathbb{R}^m ; il ne faut pas les confondre avec les paramètres de contrôle \boldsymbol\lambda (Re, Rm…) des autres modules.

La démonstration ci-dessus est complète, mais elle suppose que \lambda^\star existe ; c’est la forme du résultat que l’on utilise en pratique, et la section 4.5 dit ce qui se passe quand cette hypothèse tombe. Csiszár (1975) a développé la géométrie générale de ces projections (existence et caractérisation de la loi la plus proche dans un ensemble convexe de lois, inégalité de type Pythagore, qui devient l’égalité (iii) pour des contraintes linéaires comme ici) ; pour le cas des contraintes de moments, Cover et Thomas (2006, chapitres 11 et 12) donnent un exposé élémentaire. Khalifa, Elsahar et Dymetman (2021) appliquent cette construction à un modèle génératif de texte pré-entraîné pour imposer des contraintes de moments sur les textes produits.

Les poids \bar w_i sont exactement les poids autonormalisés de la définition 3.9 avec \tilde p = q\,e^{\lambda\cdot\mathbf{g}} et la proposition q : l’I-projection d’un modèle génératif se réalise en repondérant ses échantillons, et sa qualité se diagnostique par l’ESS du chapitre 3.

Le code suivant implémente la proposition 4.12 pour une loi de base bimodale, q = \frac12\mathcal{N}(-2, 0{,}5^2) + \frac12\mathcal{N}(2, 0{,}5^2), image simplifiée des deux polarités d’une dynamo à renversements. Il compare le multiplicateur estimé sur N = 20\,000 tirages au multiplicateur exact, calculable ici parce que l’inclinaison d’un mélange de gaussiennes est un mélange de gaussiennes (exercice 10).

import numpy as np
from scipy.optimize import brentq
rng = np.random.default_rng(3)

def tilt(G, c, iters=50, tol=1e-10):
    """I-projection empirique : poids w_i ∝ exp(lam . g(x_i)) tels que sum_i w_i g(x_i) = c.
    G : (N, m) valeurs des fonctions de contrainte ; c : (m,). Newton sur le dual
    psi(lam) = log mean exp(lam . g) - lam . c (convexe : gradient E_lam[g] - c, hessienne Cov_lam[g])."""
    N, m = G.shape
    lam = np.zeros(m)
    for k in range(iters):
        a = G @ lam
        w = np.exp(a - a.max()); w /= w.sum()           # poids autonormalisés
        mean = w @ G
        grad = mean - c
        if np.linalg.norm(grad) < tol:
            break
        H = (G - mean).T @ ((G - mean) * w[:, None])    # covariance pondérée
        step = np.linalg.solve(H + 1e-12 * np.eye(m), grad)
        t = 1.0                                         # recherche linéaire (Armijo) sur psi,
        psi = lambda l: np.log(np.mean(np.exp(G @ l - (G @ l).max()))) + (G @ l).max() - l @ c
        # avec une marge de 1e-12 pour les erreurs d'arrondi près de l'optimum
        while psi(lam - t * step) > psi(lam) - 1e-4 * t * grad @ step + 1e-12 and t > 1e-8:
            t /= 2
        lam = lam - t * step
    a = G @ lam; w = np.exp(a - a.max()); w /= w.sum()
    ess = 1.0 / np.sum(w ** 2)
    kl = np.sum(w * np.log(N * w + 1e-300))             # KL(p_lam || q_N) entre lois empiriques
    return lam, w, ess, kl, k

# Loi de base q : mélange 0,5 N(-2, 0,5^2) + 0,5 N(2, 0,5^2) (deux « polarités »)
N, s = 20_000, 0.5
x = np.where(rng.random(N) < 0.5, -2.0, 2.0) + s * rng.standard_normal(N)
def mean_exact(lam):        # moyenne de la loi inclinée exacte (calcul analytique, exercice 10)
    a = np.array([-2.0, 2.0]); logp = lam * a + lam ** 2 * s ** 2 / 2
    p = np.exp(logp - logp.max()); p /= p.sum()
    return p @ (a + lam * s ** 2)

print(" c     lambda (éch.)  lambda (exact)   ESS/N     KL     itérations")
for c in [0.0, 0.5, 1.5, 2.5, 3.0, 3.5]:
    lam, w, ess, kl, k = tilt(x[:, None], np.array([c]))
    lam_ex = brentq(lambda l: mean_exact(l) - c, -60, 60)
    print(f"{c:3.1f}   {lam[0]:9.4f}     {lam_ex:9.4f}     {ess / N:8.5f}   {kl:5.2f}     {k}")
lam, w, ess, kl, k = tilt(x[:, None], np.array([4.0]))
print(f"c = 4 : max x_i = {x.max():.3f} < 4, pas de convergence en {k + 1} itérations, lambda = {lam[0]:.1e}, ESS = {ess:.2f}")

# Deux contraintes à la fois : E[X] = 0,5 et E[X^2] = 4,5
G = np.stack([x, x ** 2], axis=1)
lam, w, ess, kl, k = tilt(G, np.array([0.5, 4.5]))
print("\ndeux contraintes : lambda =", lam.round(4), "; moments obtenus =", (w @ G).round(4),
      f"; ESS/N = {ess / N:.3f} ; KL = {kl:.3f}")
# sortie :
#  c     lambda (éch.)  lambda (exact)   ESS/N     KL     itérations
# 0.0      0.0011        0.0000      0.99999    0.00     2
# 0.5      0.1214        0.1198      0.94299    0.03     3
# 1.5      0.4328        0.4300      0.64044    0.29     5
# 2.5      2.0250        2.0053      0.18707    1.20     6
# 3.0      4.3169        4.0000      0.01898    2.76     6
# 3.5      9.3655        6.0000      0.00105    5.88     8
# c = 4 : max x_i = 3.711 < 4, pas de convergence en 50 itérations, lambda = 1.3e+13, ESS = 1.00
#
# deux contraintes : lambda = [0.1146 0.0555] ; moments obtenus = [0.5 4.5] ; ESS/N = 0.927 ; KL = 0.037

Le tableau se lit en deux régimes. Tant que c \le 1{,}5, l’inclinaison transfère de la masse d’un mode vers l’autre, les multiplicateurs estimé et exact coïncident à 3\times10^{-3} près et l’ESS reste au-dessus de 60 %. Pour c \ge 3, la contrainte exige de déplacer le mode de droite au-delà de ce que l’échantillon contient : le multiplicateur estimé s’écarte du vrai (9,4 au lieu de 6 pour c = 3{,}5), la divergence explose et l’ESS tombe à 0,1 %. Le schéma interactif 2 permet de parcourir continûment ces deux régimes.

4.5 Quand le support manque

Les conséquences pratiques sont importantes pour le stage. L’inclinaison exponentielle repondère les échantillons d’un modèle, elle ne crée pas de masse là où le modèle n’en met pas. Elle convient donc pour corriger un léger biais statistique (un bilan d’énergie violé de quelques pour cent), pas pour faire apparaître un régime absent, comme des renversements qu’un modèle n’aurait jamais produits. Quand l’ESS chute, trois remèdes existent dans la littérature récente. Le premier consiste à guider l’échantillonneur pendant la génération, par des méthodes de Monte-Carlo séquentiel qui ré-échantillonnent des particules au fil du débruitage (Wu et al., 2023 ; Singhal et al., 2025). Mun (2026) combine précisément ces ingrédients sur un problème proche du nôtre : pour un émulateur de simulations fondé sur un modèle de diffusion, il cherche la loi inclinée exponentiellement la plus proche en KL de la loi du modèle sous des contraintes de moments, et réalise cette inclinaison au moment de la génération par un tel guidage. Le deuxième remède consiste à affiner le modèle vers la loi inclinée, ce que Tauberschmidt et al. (2026) font pour un modèle de flow matching ; leur récompense est un résidu d’EDP en forme faible, calculé échantillon par échantillon, alors qu’une contrainte comme le bilan d’énergie moyen porte sur la loi entière (remarque 4.15). Le troisième remède est d’imposer la contrainte dès l’entraînement.

5. Distances entre lois et transport optimal

Pour dire qu’un modèle génératif est bon, il faut mesurer l’écart entre la loi générée q_{\boldsymbol\lambda} et la loi des données \mu_{\boldsymbol\lambda}, avec des échantillons des deux. La divergence de Kullback–Leibler s’y prête mal : elle exige des densités, et elle est infinie dès que les supports sont disjoints, même s’ils sont très proches. Les distances de transport optimal mesurent au contraire de combien il faut déplacer la masse pour passer d’une loi à l’autre ; elles tiennent compte de la géométrie de l’espace des états.

5.1 Les problèmes de Monge et de Kantorovich

Ce problème peut n’avoir aucune solution admissible : une application déterministe ne peut pas couper une masse ponctuelle en deux, si bien qu’il n’existe aucune T avec T_\#\delta_0 = \frac12\delta_{-1} + \frac12\delta_1. Kantorovich a levé cet obstacle en autorisant à répartir la masse.

Pour des lois discrètes \mu = \sum_i a_i\delta_{\mathbf{x}_i} et \nu = \sum_j b_j\delta_{\mathbf{y}_j}, un couplage est une matrice (\pi_{ij}) positive de sommes de lignes a_i et de sommes de colonnes b_j : le problème de Kantorovich est un programme linéaire. Sous des hypothèses faibles (coût continu et positif), il admet une solution (Villani, 2009, théorème 4.1).

On utilisera aussi la formule duale de Kantorovich–Rubinstein, admise : W_1(\mu, \nu) = \sup\{\mathbb{E}_\mu[f] - \mathbb{E}_\nu[f]\}, le supremum portant sur les fonctions f 1-lipschitziennes. Elle a une lecture très concrète pour l’évaluation : si W_1(q, \mu) est petit, la moyenne de toute observable lipschitzienne est proche sous les deux lois.

5.2 Le cas de la dimension 1 : les quantiles

En pratique, calculer W_p entre deux nuages 1D de même taille revient à les trier et à comparer les valeurs rang par rang, en O(N\log N).

5.3 Le cas gaussien en dimension quelconque

Cette formule a un usage très répandu : la FID (Fréchet Inception Distance) de Heusel et al. (2017), métrique standard des modèles génératifs d’images, est exactement W_2^2 entre deux gaussiennes ajustées sur des descripteurs d’images. Pour des champs physiques, on peut l’appliquer de la même façon à des descripteurs choisis (énergies par bande spectrale, hélicités), en gardant à l’esprit qu’elle ne compare que les deux premiers moments.

5.4 Wasserstein tranché

En dimension d \ge 2, calculer W_p entre deux nuages de N points demande de résoudre un programme linéaire de taille N\times N, ce qui devient coûteux au-delà de quelques milliers de points. L’idée du transport « tranché » est de ne résoudre que des problèmes 1D.

5.5 Distance d’énergie et MMD

On retrouve la même structure que W_1 = \int|F - G|, mais avec un carré : la distance d’énergie est moins sensible aux queues.

Le code suivant vérifie les formules de cette section : W_1 et W_2 par tri, comparés à SciPy et aux valeurs exactes ; la distance d’énergie (notée D_E dans le code) par la proposition 5.12, comparée à SciPy (attention : scipy.stats.energy_distance renvoie \sqrt{\mathcal{E}}) ; la formule gaussienne en dimension 2 ; enfin SW_2 et la MMD gaussienne, avec la valeur de référence obtenue entre deux échantillons indépendants de la même loi.

import numpy as np
from scipy.stats import wasserstein_distance, energy_distance, norm
from scipy.linalg import sqrtm
rng = np.random.default_rng(4)

def wasserstein_1d(x, y, p=1):
    """W_p entre deux nuages 1D de même taille : on apparie les quantiles (tri)."""
    return np.mean(np.abs(np.sort(x) - np.sort(y)) ** p) ** (1 / p)

def energy_1d(x, y):
    """Distance d'énergie D_E = 2E|X-Y| - E|X-X'| - E|Y-Y'| en O(N log N) via l'identité 1D
    E|X-Y| = intégrale de F(1-G) + G(1-F), soit D_E = 2 * intégrale de (F - G)^2."""
    z = np.sort(np.concatenate([x, y]))
    F = np.searchsorted(np.sort(x), z[:-1], side="right") / len(x)
    G = np.searchsorted(np.sort(y), z[:-1], side="right") / len(y)
    return 2 * np.sum((F - G) ** 2 * np.diff(z))

# E|Z| pour Z ~ N(m, s^2) : s sqrt(2/pi) exp(-m^2 / 2s^2) + m (1 - 2 Phi(-m/s))
Eabs = lambda m, s: s * np.sqrt(2 / np.pi) * np.exp(-m ** 2 / (2 * s ** 2)) + m * (1 - 2 * norm.cdf(-m / s))

N = 100_000
x = rng.standard_normal(N)                  # mu = N(0, 1)
y = 1.0 + 2.0 * rng.standard_normal(N)      # nu = N(1, 2^2)
print(f"W1 : quantiles = {wasserstein_1d(x, y):.4f} ; scipy = {wasserstein_distance(x, y):.4f} ; exact E|1+Z| = {Eabs(1, 1):.4f}")
print(f"W2 : quantiles = {wasserstein_1d(x, y, 2):.4f} ; formule gaussienne sqrt(1^2 + (2-1)^2) = {np.sqrt(2):.4f}")
print(f"D_E : 2*int (F-G)^2 = {energy_1d(x, y):.4f} ; scipy energy_distance^2 = {energy_distance(x, y) ** 2:.4f}")
print(f"D_E exacte = {2 * Eabs(1, np.sqrt(5)) - Eabs(0, np.sqrt(2)) - Eabs(0, np.sqrt(8)):.4f}")

# W2 entre gaussiennes en dimension d (formule de Fréchet, Dowson et Landau 1982)
def w2_gauss(m1, S1, m2, S2):
    r = sqrtm(S2)
    c = np.real(sqrtm(r @ S1 @ r))
    return np.sqrt(np.sum((m1 - m2) ** 2) + np.trace(S1 + S2 - 2 * c))
S1 = np.array([[1.0, 0.0], [0.0, 1.0]]); S2 = np.array([[2.0, 1.2], [1.2, 1.0]])
print(f"\nW2 gaussien 2D : {w2_gauss(np.zeros(2), S1, np.array([1.0, -2.0]), S2):.4f}")

# Wasserstein tranché (sliced) et MMD gaussienne en 2D, avec la référence « même loi »
def sliced_w2(X, Y, n_dir=200):
    th = rng.standard_normal((n_dir, X.shape[1])); th /= np.linalg.norm(th, axis=1, keepdims=True)
    return np.sqrt(np.mean([wasserstein_1d(X @ t, Y @ t, 2) ** 2 for t in th]))
def mmd2(X, Y, h=1.0):
    k = lambda A, B: np.exp(-((A[:, None, :] - B[None, :, :]) ** 2).sum(-1) / (2 * h ** 2))
    n, m = len(X), len(Y)
    kxx, kyy = k(X, X), k(Y, Y)
    return ((kxx.sum() - np.trace(kxx)) / (n * (n - 1)) + (kyy.sum() - np.trace(kyy)) / (m * (m - 1))
            - 2 * k(X, Y).mean())                   # estimateur sans biais
n = 2000
L = np.linalg.cholesky(S2)
A = rng.standard_normal((n, 2)); A2 = rng.standard_normal((n, 2))   # deux échantillons de N(0, I)
B = rng.standard_normal((n, 2)) @ L.T                                # N(0, S2)
print(f"SW2 : même loi = {sliced_w2(A, A2):.4f} ; N(0,I) vs N(0,S2) = {sliced_w2(A, B):.4f}")
print(f"MMD^2 : même loi = {mmd2(A, A2):.5f} ; N(0,I) vs N(0,S2) = {mmd2(A, B):.5f}")
# sortie :
# W1 : quantiles = 1.1702 ; scipy = 1.1702 ; exact E|1+Z| = 1.1666
# W2 : quantiles = 1.4193 ; formule gaussienne sqrt(1^2 + (2-1)^2) = 1.4142
# D_E : 2*int (F-G)^2 = 0.5379 ; scipy energy_distance^2 = 0.5379
# D_E exacte = 0.5341
#
# W2 gaussien 2D : 2.3998
# SW2 : même loi = 0.0475 ; N(0,I) vs N(0,S2) = 0.4485
# MMD^2 : même loi = -0.00043 ; N(0,I) vs N(0,S2) = 0.04556

Les estimations à partir de 10^5 tirages s’écartent des valeurs exactes de quelques 10^{-3}, comme attendu au chapitre 3. La dernière partie illustre une règle d’usage : une distance estimée n’a de sens que comparée à sa valeur entre deux échantillons indépendants de la même loi (ici SW_2 \approx 0{,}05 et \mathrm{MMD}^2 \approx 0), qui fixe le niveau de bruit.

5.6 Quelle métrique pour comparer des champs générés ?

Un champ MHD discrétisé vit en dimension D \sim 10^5 à 10^7, et l’on ne dispose que de quelques centaines ou milliers d’échantillons. Aucune distance ne peut, dans ces conditions, comparer finement les lois jointes complètes. La bonne pratique consiste à combiner des métriques ciblées, chacune sensible à un aspect physique précis, et à toujours les rapporter au niveau de bruit « même loi ». Le tableau résume ce que voit chaque famille.

Lienen et al. (2024) proposent dans ce même esprit une distance W_2 entre ensembles de spectres d’énergie cinétique turbulente, pour éviter les descripteurs appris de type FID. Le code ci-dessous rend ces limites concrètes sur des champs 1D périodiques gaussiens de spectre E(k) \propto k^{-5/3} (les « données »), comparés à quatre « modèles » : un modèle parfait, un modèle au spectre trop raide, un modèle au bon spectre mais à loi ponctuelle déformée, et le champ nul.

import numpy as np
rng = np.random.default_rng(5)
n, M = 128, 2000                                   # points de grille, nombre de champs
k = np.fft.rfftfreq(n, 1 / n)                      # nombres d'onde 0..64

def champs(alpha, M, transforme=False):
    """Champs gaussiens périodiques 1D de spectre E(k) ~ k^-alpha, normalisés à variance 1."""
    amp = np.zeros_like(k); amp[1:] = k[1:] ** (-alpha / 2)
    c = amp * (rng.standard_normal((M, len(k))) + 1j * rng.standard_normal((M, len(k))))
    u = np.fft.irfft(c, n)
    u /= u.std()
    if transforme:                                 # loi à un point non gaussienne
        u = np.sign(u) * np.abs(u) ** 1.5; u /= u.std()
    return u

def spectre(u):
    return np.mean(np.abs(np.fft.rfft(u, axis=1)) ** 2, axis=0)[1:]

def energy(X, Y):                                  # distance d'énergie entre nuages de R^n
    d = lambda A, B: np.sqrt(((A[:, None, :] - B[None, :, :]) ** 2).sum(-1))
    nx, ny = len(X), len(Y)
    return 2 * d(X, Y).mean() - d(X, X).sum() / (nx * (nx - 1)) - d(Y, Y).sum() / (ny * (ny - 1))

ref = champs(5 / 3, M)                             # « données »
modeles = {"A : même loi": champs(5 / 3, M),
           "B : spectre k^-3": champs(3.0, M),
           "C : loi à 1 point déformée": champs(5 / 3, M, transforme=True),
           "D : champ nul": np.zeros((M, n))}
E_ref = spectre(ref)
print(f"{'modèle':28s} {'MSE apparié':>11s} {'W1 1 point':>10s} {'écart log-spectre':>17s} {'D_E champs':>10s}")
for nom, u in modeles.items():
    mse = np.mean((u - ref) ** 2)
    w1 = np.mean(np.abs(np.sort(u.ravel()) - np.sort(ref.ravel())))
    dspec = np.mean(np.abs(np.log(spectre(u) + 1e-12) - np.log(E_ref)))
    de = energy(u[:500], ref[:500])
    print(f"{nom:28s} {mse:11.3f} {w1:10.3f} {dspec:17.3f} {de:10.3f}")
# sortie :
# modèle                       MSE apparié W1 1 point écart log-spectre D_E champs
# A : même loi                       2.002      0.006             0.028      0.010
# B : spectre k^-3                   1.996      0.004             3.792      0.228
# C : loi à 1 point déformée         2.028      0.179             0.119      0.012
# D : champ nul                      1.000      0.798            30.568      6.403

Chaque colonne a son angle mort. L’erreur quadratique appariée classe le champ nul (1,0) devant le modèle parfait (2,0) : pour deux champs indépendants de même loi, \mathbb{E}|\mathbf{x} - \mathbf{x}'|^2 = 2\,\mathrm{Var}, si bien que cette métrique pénalise la variabilité qu’un modèle génératif doit justement reproduire. La distance W_1 sur les valeurs ponctuelles ne distingue pas le spectre trop raide (0,004), puisque les deux lois à un point sont des \mathcal{N}(0, 1). L’écart de spectres repère B mais voit à peine C. La distance d’énergie sur les champs entiers détecte B mais confond C avec le modèle parfait (0,012 contre 0,010) : avec 500 champs en dimension 128, elle manque de puissance pour un défaut de loi marginale. Seule la combinaison de métriques ciblées classe correctement les quatre modèles.

6. Processus stochastiques, Fokker–Planck et transport de lois

Trois objets du stage relèvent de ce chapitre. Une dynamique physique forcée aléatoirement, \dot{\mathbf{x}} = F_{\boldsymbol\lambda}(\mathbf{x}) + \sigma\,\dot W, fait évoluer une loi selon l’équation de Fokker–Planck, dont la solution stationnaire est la loi invariante que l’on veut apprendre. Un modèle de diffusion est lui-même une équation différentielle stochastique (Song et al., 2021). Enfin, le flow matching transporte le bruit vers les données par une équation différentielle ordinaire, et la loi transportée obéit à l’équation de continuité. Ce chapitre relie ces trois points de vue.

6.1 Mouvement brownien

Le point (ii) est le cœur du calcul stochastique : un accroissement brownien sur un pas \mathrm{d}t est d’ordre \sqrt{\mathrm{d}t}, et son carré (\mathrm{d}W)^2 se comporte comme \mathrm{d}t, de façon déterministe. Il donne aussi l’intuition d’un résultat classique, admis ici : les trajectoires browniennes ne sont presque sûrement dérivables nulle part, puisque le taux d’accroissement \Delta W/\Delta t est d’ordre 1/\sqrt{\Delta t}.

6.2 Équations différentielles stochastiques et schéma d’Euler–Maruyama

Pour \theta = 2, \sigma = 1, la variance exacte est 0{,}25 ; avec \Delta t = 0{,}1, le schéma converge vers 0{,}278, soit 11 % d’erreur. Ce biais en O(\Delta t) est l’ordre faible 1 annoncé ; il faut l’avoir en tête quand on fabrique des données d’entraînement par intégration numérique.

6.3 Formule d’Itô

6.4 L’équation de Fokker–Planck

Pour Ornstein–Uhlenbeck, V(x) = \theta x^2/2 et p_\infty \propto e^{-\theta x^2/\sigma^2} : c’est bien \mathcal{N}(0, \sigma^2/(2\theta)). La forme faible (i) est la plus utile pour le stage, pour deux raisons. Elle ne demande pas de connaître p_\infty, seulement de savoir évaluer \mathcal{L}\varphi sur des échantillons : c’est une contrainte vérifiable sur une loi générée. Et elle reste valable quand la loi invariante n’a pas de densité, ce qui est le cas des attracteurs chaotiques dissipatifs sans bruit. Avec \varphi égale à l’énergie, elle exprime l’égalité en moyenne de l’injection et de la dissipation ; ce sont ces « identités de stationnarité » que l’expérience E1 teste sur la loi générée du modèle de Gissinger (voir le rapport).

Le code suivant simule Ornstein–Uhlenbeck par Euler–Maruyama, compare les moments empiriques à la solution exacte de Fokker–Planck à trois instants, puis vérifie le biais de pas de temps de la proposition 6.6.

import numpy as np
rng = np.random.default_rng(6)

# Ornstein–Uhlenbeck : dX = -theta X dt + sigma dW, X_0 = x0 ; loi stationnaire N(0, sigma^2 / (2 theta))
theta, sigma, x0 = 2.0, 1.0, 1.5
v_inf = sigma ** 2 / (2 * theta)

def euler_maruyama(dt, T, N, instants=()):
    X, sorties = np.full(N, x0), {}
    for n in range(1, int(round(T / dt)) + 1):
        X = X - theta * X * dt + sigma * np.sqrt(dt) * rng.standard_normal(N)
        for t in instants:
            if n == int(round(t / dt)):
                sorties[t] = X.copy()
    return X, sorties

N = 100_000
X, sorties = euler_maruyama(2e-3, 4.0, N, instants=(0.25, 1.0, 4.0))
for t, Xt in sorties.items():
    m_ex = x0 * np.exp(-theta * t)                       # solution exacte de Fokker–Planck :
    v_ex = v_inf * (1 - np.exp(-2 * theta * t))          # p(., t) = N(m_ex, v_ex)
    print(f"t = {t:4.2f} : moyenne {Xt.mean():.4f} (exacte {m_ex:.4f}), variance {Xt.var():.4f} (exacte {v_ex:.4f})")

# Histogramme à t = 4 comparé à la densité stationnaire p_inf
h, e = np.histogram(X, bins=60, range=(-2, 2), density=True)
c = (e[:-1] + e[1:]) / 2
p_inf = np.exp(-c ** 2 / (2 * v_inf)) / np.sqrt(2 * np.pi * v_inf)
print(f"écart max histogramme / p_inf : {np.abs(h - p_inf).max():.3f} (p_inf(0) = {p_inf.max():.3f})")

# Biais de pas de temps : la chaîne d'Euler–Maruyama a pour variance stationnaire
# sigma^2 dt / (1 - (1 - theta dt)^2) = sigma^2 / (theta (2 - theta dt)), et non sigma^2 / (2 theta)
for dt in [0.01, 0.1, 0.4]:
    X, _ = euler_maruyama(dt, 8.0, N)
    print(f"dt = {dt:4.2f} : variance {X.var():.4f} ; prédiction discrète {sigma ** 2 / (theta * (2 - theta * dt)):.4f} ; continue {v_inf:.4f}")
# sortie :
# t = 0.25 : moyenne 0.9105 (exacte 0.9098), variance 0.1575 (exacte 0.1580)
# t = 1.00 : moyenne 0.2041 (exacte 0.2030), variance 0.2443 (exacte 0.2454)
# t = 4.00 : moyenne 0.0006 (exacte 0.0005), variance 0.2506 (exacte 0.2500)
# écart max histogramme / p_inf : 0.017 (p_inf(0) = 0.796)
# dt = 0.01 : variance 0.2534 ; prédiction discrète 0.2525 ; continue 0.2500
# dt = 0.10 : variance 0.2800 ; prédiction discrète 0.2778 ; continue 0.2500
# dt = 0.40 : variance 0.4185 ; prédiction discrète 0.4167 ; continue 0.2500

6.5 Équation de continuité et transport de lois par un champ de vitesse

Supprimons maintenant le bruit : \sigma = 0. Les particules suivent l’équation différentielle ordinaire \dot{\mathbf{x}} = v(\mathbf{x}, t), et l’équation de Fokker–Planck se réduit à l’équation de continuité de la mécanique des fluides, où la densité de probabilité joue le rôle de la masse volumique.

C’est la version continue de la formule du jacobien (corollaire 2.6) : un champ qui dilate les volumes (\nabla\cdot v > 0) fait baisser la densité le long des trajectoires, et un champ à divergence nulle conserve la densité de chaque particule, comme un écoulement incompressible conserve la masse volumique d’une particule fluide. Chen et al. (2018) en tirent les flots normalisants continus, dont la vraisemblance s’obtient en intégrant \nabla\cdot v le long des trajectoires.

Ce principe est celui de Lipman et al. (2023), de Liu, Gong et Liu (2023, « rectified flow ») et d’Albergo et Vanden-Eijnden (2023, interpolants stochastiques). Les trajectoires individuelles \mathbf{x}_t = (1 - t)\mathbf{x}_0 + t\,\mathbf{x}_1 sont des droites qui se croisent ; les trajectoires du flot marginal ne se croisent pas (pour un champ lipschitzien, deux solutions distinctes d’une même EDO ne peuvent se rencontrer), mais elles ont les mêmes lois à chaque instant. Le code ci-dessous le vérifie pour un bruit \mathcal{N}(0, 1) et une cible bimodale, pour laquelle v se calcule exactement : sachant la composante k du mélange, le couple (\mathbf{x}_t, \mathbf{x}_1 - \mathbf{x}_0) est gaussien et la proposition 1.16 donne l’espérance conditionnelle.

import numpy as np
rng = np.random.default_rng(7)

# Source x0 ~ N(0, 1) ; cible x1 ~ 0,5 N(-2, 0,4^2) + 0,5 N(2, 0,4^2) ; x_t = (1 - t) x0 + t x1
pis, ms, ss = np.array([0.5, 0.5]), np.array([-2.0, 2.0]), np.array([0.4, 0.4])

def gauss(x, m, v):
    return np.exp(-(x - m) ** 2 / (2 * v)) / np.sqrt(2 * np.pi * v)

def p_t(x, t):                      # loi de x_t sachant la composante k : N(t m_k, (1-t)^2 + t^2 s_k^2)
    x = np.asarray(x)[..., None]
    return np.sum(pis * gauss(x, t * ms, (1 - t) ** 2 + t ** 2 * ss ** 2), axis=-1)

def v(x, t):
    """Vitesse marginale v(x, t) = E[x1 - x0 | x_t = x] (régression gaussienne par composante)."""
    x = np.asarray(x)[..., None]
    var = (1 - t) ** 2 + t ** 2 * ss ** 2               # Var(x_t | k)
    cov = t * ss ** 2 - (1 - t)                          # Cov(x_t, x1 - x0 | k)
    post = pis * gauss(x, t * ms, var)                   # P(k | x_t = x), non normalisée
    v_k = ms + cov / var * (x - t * ms)                  # E[x1 - x0 | x_t = x, k]
    return np.sum(post * v_k, axis=-1) / np.sum(post, axis=-1)

def transporte(z, n_pas=100):
    """Intègre dx/dt = v(x, t) de t = 0 à t = 1 (Runge–Kutta d'ordre 4)."""
    tau = 1.0 / n_pas
    for i in range(n_pas):
        t = i * tau
        k1 = v(z, t); k2 = v(z + tau / 2 * k1, t + tau / 2)
        k3 = v(z + tau / 2 * k2, t + tau / 2); k4 = v(z + tau * k3, t + tau)
        z = z + tau / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    return z

# 1) Vérification de l'équation de continuité d_t p + d_x (p v) = 0 par différences finies
x = np.linspace(-4, 4, 2001); h = x[1] - x[0]; dt = 1e-5
for t in [0.2, 0.5, 0.8]:
    dpdt = (p_t(x, t + dt) - p_t(x, t - dt)) / (2 * dt)
    div = np.gradient(p_t(x, t) * v(x, t), h)
    print(f"t = {t} : max|d_t p + d_x(p v)| = {np.abs(dpdt + div).max():.1e}  (max|d_t p| = {np.abs(dpdt).max():.2f})")

# 2) Transport de 50 000 particules du bruit vers la cible par l'EDO
z = transporte(rng.standard_normal(50_000))
cible = np.where(rng.random(50_000) < 0.5, -2.0, 2.0) + 0.4 * rng.standard_normal(50_000)
print(f"W1(particules transportées, cible) = {np.mean(np.abs(np.sort(z) - np.sort(cible))):.4f}")
print(f"fraction à droite : {np.mean(z > 0):.4f} ; moyenne de |x| : {np.abs(z).mean():.3f} (cible {np.abs(cible).mean():.3f})")

# 3) L'application x0 -> x1 obtenue est croissante : les trajectoires ne se croisent pas
z0 = np.linspace(-3, 3, 7)
print("x0 :", z0.round(2), "\nx1 :", transporte(z0).round(2))
# sortie :
# t = 0.2 : max|d_t p + d_x(p v)| = 2.8e-06  (max|d_t p| = 0.15)
# t = 0.5 : max|d_t p + d_x(p v)| = 2.3e-05  (max|d_t p| = 1.38)
# t = 0.8 : max|d_t p + d_x(p v)| = 7.4e-05  (max|d_t p| = 1.72)
# W1(particules transportées, cible) = 0.0173
# fraction à droite : 0.4976 ; moyenne de |x| : 1.998 (cible 2.002)
# x0 : [-3. -2. -1.  0.  1.  2.  3.]
# x1 : [-3.11 -2.68 -2.19  0.    2.19  2.68  3.11]

Le résidu de l’équation de continuité est de l’ordre de 10^{-5}, soit l’erreur de troncature des différences finies, pour des dérivées temporelles d’ordre 1. Les 5\times10^4 particules transportées sont à W_1 = 0{,}017 de la cible, et l’application \mathbf{x}_0 \mapsto \mathbf{x}(1) est croissante : en dimension 1, le flot marginal réalise ici une application monotone, donc l’application de transport optimal du théorème 5.6. Le point \mathbf{x}_0 = 0 reste en 0 par symétrie : c’est un point d’équilibre instable du flot, qui sépare les particules envoyées vers chacun des deux modes.

Exercices

Les exercices sont gradués : ★ pour une application directe du cours, ★★ pour un calcul ou une démonstration en plusieurs étapes, ★★★ pour un exercice proche des questions du stage. Cherchez chaque exercice avant de lire sa solution.

1. ★ Loi exponentielle. Soit X de densité a\,e^{-ax} sur [0, \infty), avec a > 0. Calculer \mathbb{P}(X > t), \mathbb{E}[X] et \mathrm{Var}(X), puis montrer l’absence de mémoire : \mathbb{P}(X > s + t \mid X > s) = \mathbb{P}(X > t).

2. ★ Covariance nulle sans indépendance. Soit X \sim \mathcal{N}(0, 1) et Y = X^2. Montrer que \mathrm{Cov}(X, Y) = 0, puis que X et Y ne sont pas indépendantes.

3. ★★ Conditionnement gaussien. Avec \mathbf{m} = (1, -2) et \Sigma = \begin{pmatrix}2 & 1{,}2\\ 1{,}2 & 1\end{pmatrix}, (a) vérifier que \Sigma est définie positive ; (b) montrer que R = X_1 - 1{,}2\,X_2 est indépendante de X_2 et calculer sa loi ; (c) en déduire la loi de X_1 sachant X_2 = -1 et retrouver l’exemple 1.17 ; (d) calculer le facteur de Cholesky L.

4. ★ Loi log-normale. Soit X \sim \mathcal{N}(m, s^2) et Y = e^X. Calculer la densité de Y par la formule du jacobien, puis \mathbb{E}[Y].

5. ★★ Méthode de Box–Muller. Soit U_1, U_2 indépendantes uniformes sur (0, 1), R = \sqrt{-2\log U_1} et \Theta = 2\pi U_2. Montrer que (Z_1, Z_2) = (R\cos\Theta, R\sin\Theta) suit \mathcal{N}(0, I_2).

6. ★ Monte-Carlo pour \pi. On tire N points uniformes dans [0, 1]^2 et l’on estime \pi par 4 fois la proportion de points tels que x^2 + y^2 \le 1. Combien de points faut-il pour que la demi-largeur de l’intervalle de confiance à 95 % soit 10^{-3} ?

7. ★★ ESS d’une proposition décalée, et fléau de la dimension. Cible p = \mathcal{N}(0, I_d), proposition q = \mathcal{N}(\mathbf{a}, I_d). (a) Calculer w = p/q, \mathbb{E}_q[w] et \mathbb{E}_q[w^2]. (b) En déduire la limite de \mathrm{ESS}/N. (c) Application : un modèle génératif de champs à d = 10^4 points de grille dont chaque valeur est biaisée de 0{,}1 écart-type. Que vaut l’ESS limite si l’on veut corriger ce biais par repondération ?

8. ★★ Divergence de Kullback–Leibler entre gaussiennes. Démontrer la formule \mathrm{KL}(\mathcal{N}(m_1, s_1^2)\,\|\,\mathcal{N}(m_2, s_2^2)) = \log\frac{s_2}{s_1} + \frac{s_1^2 + (m_1 - m_2)^2}{2s_2^2} - \frac12, puis calculer les deux divergences entre \mathcal{N}(0, 1) et \mathcal{N}(0, 4).

9. ★★ Maximum d’entropie. (a) Sur [0, \infty), trouver la densité d’entropie maximale de moyenne m. (b) Sur \mathbb{R}, trouver la densité d’entropie maximale de moyenne \mu et de variance v. (c) Sur un ensemble fini d’états \{1, \dots, K\} d’énergies E_k, trouver la loi d’entropie maximale d’énergie moyenne fixée.

10. ★★★ I-projections explicites. (a) Soit q = \mathcal{N}(0, 1). Calculer l’I-projection sous la contrainte \mathbb{E}[X^2] = v (v > 0), le multiplicateur \lambda^\star et la divergence minimale. (b) Montrer que \mathcal{N}(m, s^2) inclinée par e^{\lambda x} est proportionnelle à \mathcal{N}(m + \lambda s^2, s^2), avec la constante e^{\lambda m + \lambda^2s^2/2}. En déduire l’inclinaison d’un mélange de gaussiennes, et la moyenne inclinée de q = \frac12\mathcal{N}(-2, s^2) + \frac12\mathcal{N}(2, s^2). (c) Avec s = 0{,}5, montrer que \mathrm{ESS}/N \to \frac12e^{-\lambda^2/4} quand \lambda est grand, et expliquer l’écart entre multiplicateur estimé et exact pour c = 3{,}5 dans le code du chapitre 4.

11. ★★ Calculs de distances de Wasserstein. (a) Calculer W_1 et W_2 entre \mu = \frac12\delta_0 + \frac12\delta_1 et \nu = \frac12\delta_0 + \frac12\delta_3, et comparer au coût du couplage indépendant. (b) Calculer W_1 et W_2 entre \delta_0 et \frac12\delta_{-1} + \frac12\delta_1. (c) Pour une translation \nu = \mu(\cdot - a) en dimension 1, montrer que W_p(\mu, \nu) = |a| pour tout p.

12. ★★ Distance d’énergie. (a) Calculer \mathcal{E}(\delta_0, \delta_a) directement puis par la proposition 5.12. (b) Pour \mu = \mathcal{N}(0, 1) et \nu = \mathcal{N}(m, 1), montrer que \mathcal{E} \approx m^2/\sqrt\pi quand m \to 0. Comparer à W_1 et à la KL.

13. ★★ Calculs d’Itô. (a) Calculer \mathrm{d}(e^{W_t}) et en déduire \mathbb{E}[e^{W_t}]. (b) Pour le mouvement brownien géométrique \mathrm{d}X = aX\,\mathrm{d}t + bX\,\mathrm{d}W, calculer \mathrm{d}(\log X_t) et la solution. Comparer le taux de croissance de \mathbb{E}[X_t] à celui d’une trajectoire typique. (c) Vérifier que p(x, t) = \mathcal{N}(m(t), v(t)) résout l’équation de Fokker–Planck d’Ornstein–Uhlenbeck si et seulement si m' = -\theta m et v' = -2\theta v + \sigma^2.

14. ★★★ Double puits et identité de stationnarité (code). Soit \mathrm{d}X = -V'(X)\,\mathrm{d}t + \sigma\,\mathrm{d}W avec V(x) = \frac{x^4}{4} - \frac{x^2}{2} et \sigma = 0{,}8. (a) Donner la densité stationnaire et le rapport p_\infty(1)/p_\infty(0). (b) En appliquant la forme faible \mathbb{E}_{p_\infty}[\mathcal{L}\varphi] = 0 à \varphi(x) = x^2, établir une identité entre \mathbb{E}[X^4] et \mathbb{E}[X^2] qui ne demande pas de connaître p_\infty. (c) Vérifier (a) et (b) par Euler–Maruyama et par quadrature.

import numpy as np
from scipy.integrate import quad
rng = np.random.default_rng(8)

# dX = -V'(X) dt + sigma dW avec V(x) = x^4/4 - x^2/2 ; p_inf ∝ exp(-2 V / sigma^2)
sigma, dt, T, N = 0.8, 2e-3, 10.0, 10_000
V = lambda x: x ** 4 / 4 - x ** 2 / 2
X = rng.standard_normal(N)
for _ in range(int(T / dt)):
    X += -(X ** 3 - X) * dt + sigma * np.sqrt(dt) * rng.standard_normal(N)

Z = quad(lambda x: np.exp(-2 * V(x) / sigma ** 2), -np.inf, np.inf)[0]
mom = lambda k: quad(lambda x: x ** k * np.exp(-2 * V(x) / sigma ** 2), -np.inf, np.inf)[0] / Z
print(f"E[X^2] : EM = {np.mean(X ** 2):.4f} ; p_inf = {mom(2):.4f}")
print(f"E[X^4] : EM = {np.mean(X ** 4):.4f} ; p_inf = {mom(4):.4f}")
print(f"identité E[X^4] - E[X^2] = sigma^2/2 = {sigma ** 2 / 2:.4f} : EM = {np.mean(X ** 4 - X ** 2):.4f} ; p_inf = {mom(4) - mom(2):.4f}")
print(f"p_inf(1) / p_inf(0) = {np.exp(-2 * (V(1) - V(0)) / sigma ** 2):.4f} = exp(1/(2 sigma^2))")
# sortie :
# E[X^2] : EM = 0.8455 ; p_inf = 0.8447
# E[X^4] : EM = 1.1655 ; p_inf = 1.1647
# identité E[X^4] - E[X^2] = sigma^2/2 = 0.3200 : EM = 0.3200 ; p_inf = 0.3200
# p_inf(1) / p_inf(0) = 2.1842 = exp(1/(2 sigma^2))

15. ★★★ Flow matching entre deux gaussiennes. Soit x_0 \sim \mathcal{N}(0, 1) et x_1 \sim \mathcal{N}(m, s^2) indépendants, x_t = (1 - t)x_0 + t\,x_1 et \sigma_t^2 = (1 - t)^2 + t^2s^2. (a) Donner la loi p_t. (b) Calculer la vitesse marginale v(x, t) = \mathbb{E}[x_1 - x_0 \mid x_t = x]. (c) Vérifier que x(t) = tm + \sigma_tx_0 résout \dot x = v(x, t), et identifier l’application x_0 \mapsto x(1). (d) Montrer que la valeur minimale de la perte de flow matching au temps t vaut s^2/\sigma_t^2 ; que faut-il en conclure sur la courbe d’apprentissage ?

Lien avec le stage

Le stage MilaDy demande des modèles génératifs conditionnels (diffusion, flow matching) de champs MHD paramétriques, qui respectent la divergence nulle, le bilan d’énergie et les conditions aux limites, et qui soient évalués par des spectres, des hélicités, la divergence et des indicateurs topologiques. Chaque chapitre de ce module intervient à un endroit précis.

• Ce qu’est un modèle génératif (chapitres 2 et 6). La loi générée est une mesure image (T_\theta)_\#\mathcal{N}(0, I) ; en flow matching, T_\theta est le flot d’une vitesse v_\theta entraînée par la régression de la proposition 6.15. Le module Modèles génératifs : diffusion et flow matching part de là.
• Expérience E1 (chapitres 2, 3 et 6). Le flow matching conditionnel y apprend la loi invariante du modèle à trois modes de Gissinger (2012). Trois tests de ce module y sont employés : la symétrie S_\#q = q avec S(Q, D, V) = (-Q, -D, V) (exemple 2.3), les identités de stationnarité \mathbb{E}_q[\mathcal{L}\varphi] = 0 (proposition 6.11) et la fréquence des renversements par la formule de Rice (remarque 6.12). Les données de référence sont prélevées le long de trajectoires, donc corrélées dans le temps : la remarque 3.6 rappelle qu’elles valent moins de tirages indépendants que leur nombre ne le laisse croire, ce qu’il faut garder en tête en lisant les écarts entre loi générée et référence. Les résultats sont dans le rapport, section E1.
• Expérience E2 (chapitres 4 et 5). Des champs MHD 2D forcés sont générés conditionnellement à la diffusivité \eta, et cinq mécanismes de contrainte de divergence sont comparés (libre, pénalité, projection finale, flot projeté, potentiels). Le volet « contraintes en loi » compare une I-projection (théorème 4.11) réalisée par repondération sur trois identités (énergie, hélicité croisée, \langle a^2\rangle), dont l’ESS mesure le coût (propositions 3.11 et 4.14), au recalage par échantillon et à une perte de stationnarité à l’entraînement (V5), seule à améliorer la loi. Les métriques suivent la section 5.6. Voir le rapport, section E2.
• Expérience E3. Le taux de croissance de la dynamo cinématique ABC en fonction de Rm est un problème aux valeurs propres déterministe, et ce module n’y sert pas directement. Il servirait dans un prolongement que E3 ne réalise pas : générer des champs de vitesse aléatoires autour de l’écoulement ABC et comparer les lois des taux de croissance obtenus. Voir le rapport, section E3.
• Formulation et taxonomie des contraintes. La distinction entre contraintes par échantillon (divergence nulle de chaque champ) et contraintes linéaires en la loi (bilan d’énergie moyen, identités de stationnarité, remarque 4.15) structure la formulation du rapport et sa taxonomie des contraintes.
• Sujet de thèse. Le sujet de thèse prolonge le stage. Sa partie théorique peut s’appuyer sur l’I-projection (chapitre 4) et sur la forme faible de l’équation de Fokker–Planck (chapitre 6), et son plan expérimental sur les métriques en loi du chapitre 5.

Les modules suivants du parcours reprennent ces outils : Apprentissage profond pour la physique pour l’entraînement des réseaux, et Modèles génératifs sous contraintes physiques pour la combinaison des contraintes par échantillon et en loi.

Pour aller plus loin

• G. Peyré et M. Cuturi, Computational Optimal Transport (2019). La référence pratique du transport optimal numérique : programmes linéaires, régularisation entropique (algorithme de Sinkhorn), cas 1D et gaussien, transport tranché. Les chapitres 2 et 3 prolongent directement le chapitre 5.
• F. Santambrogio, Optimal Transport for Applied Mathematicians (2015) ; C. Villani, Optimal Transport: Old and New (2009). Le premier est l’exposé le plus lisible de la théorie (Monge, Kantorovich, Brenier, formulation de Benamou–Brenier, flots de gradient) ; le second est l’ouvrage de référence, à consulter pour les démonstrations admises ici.
• A. B. Owen, Monte Carlo theory, methods and examples (2013). Livre en accès libre, très clair sur l’échantillonnage d’importance, les poids autonormalisés et l’ESS ; à compléter par Elvira, Martino et Robert (2022) pour une discussion critique de l’ESS.
• G. A. Pavliotis, Stochastic Processes and Applications (2014). Mouvement brownien, EDS, équation de Fokker–Planck, processus d’Ornstein–Uhlenbeck et convergence vers l’équilibre, avec un point de vue de physicien-mathématicien ; Higham (2001) pour une introduction pratique à Euler–Maruyama.
• G. Papamakarios et al., « Normalizing Flows for Probabilistic Modeling and Inference » (2021), puis Y. Lipman et al., « Flow Matching for Generative Modeling » (2023). La revue fait le pont entre la formule du jacobien et les architectures modernes ; l’article de flow matching montre comment la vitesse marginale de la proposition 6.15 s’apprend en pratique. C’est la transition naturelle vers le module p7.
• T. M. Cover et J. A. Thomas, Elements of Information Theory (2006) ; I. Csiszár (1975). Entropie, divergence de Kullback–Leibler, maximum d’entropie (chapitre 12 de Cover et Thomas) ; l’article de Csiszár est la source géométrique de l’I-projection et de l’identité de Pythagore.

Sources

Références vérifiées (DOI, arXiv ou page de l’éditeur consultés).

• Albergo, M. S. et Vanden-Eijnden, E. (2023). Building Normalizing Flows with Stochastic Interpolants. ICLR 2023. arXiv:2209.15571. https://arxiv.org/abs/2209.15571
• Benamou, J.-D. et Brenier, Y. (2000). A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem. Numerische Mathematik 84, 375–393. https://doi.org/10.1007/s002110050002
• Bonneel, N., Rabin, J., Peyré, G. et Pfister, H. (2015). Sliced and Radon Wasserstein Barycenters of Measures. Journal of Mathematical Imaging and Vision 51(1), 22–45. https://doi.org/10.1007/s10851-014-0506-3
• Chen, R. T. Q., Rubanova, Y., Bettencourt, J. et Duvenaud, D. (2018). Neural Ordinary Differential Equations. NeurIPS 2018. arXiv:1806.07366. https://arxiv.org/abs/1806.07366
• Cover, T. M. et Thomas, J. A. (2006). Elements of Information Theory, 2e éd. Wiley. https://doi.org/10.1002/047174882X
• Csiszár, I. (1975). I-Divergence Geometry of Probability Distributions and Minimization Problems. The Annals of Probability 3(1), 146–158. https://doi.org/10.1214/aop/1176996454
• Dowson, D. C. et Landau, B. V. (1982). The Fréchet distance between multivariate normal distributions. Journal of Multivariate Analysis 12(3), 450–455. https://doi.org/10.1016/0047-259X(82)90077-X
• Durrett, R. (2019). Probability: Theory and Examples, 5e éd. Cambridge University Press. https://doi.org/10.1017/9781108591034
• Elvira, V., Martino, L. et Robert, C. P. (2022). Rethinking the Effective Sample Size. International Statistical Review 90(3), 525–550. https://doi.org/10.1111/insr.12500
• Gissinger, C. (2012). A new deterministic model for chaotic reversals. The European Physical Journal B 85(4), 137. https://doi.org/10.1140/epjb/e2012-20799-5
• Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B. et Smola, A. (2012). A Kernel Two-Sample Test. Journal of Machine Learning Research 13(25), 723–773. https://jmlr.org/papers/v13/gretton12a.html
• Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B. et Hochreiter, S. (2017). GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium. NeurIPS 2017. arXiv:1706.08500. https://arxiv.org/abs/1706.08500
• Higham, D. J. (2001). An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations. SIAM Review 43(3), 525–546. https://doi.org/10.1137/S0036144500378302
• Jaynes, E. T. (1957). Information Theory and Statistical Mechanics. Physical Review 106(4), 620–630. https://doi.org/10.1103/PhysRev.106.620
• Jordan, R., Kinderlehrer, D. et Otto, F. (1998). The Variational Formulation of the Fokker–Planck Equation. SIAM Journal on Mathematical Analysis 29(1), 1–17. https://doi.org/10.1137/S0036141096303359
• Khalifa, M., Elsahar, H. et Dymetman, M. (2021). A Distributional Approach to Controlled Text Generation. ICLR 2021. arXiv:2012.11635. https://arxiv.org/abs/2012.11635
• Kloeden, P. E. et Platen, E. (1992). Numerical Solution of Stochastic Differential Equations. Springer. https://doi.org/10.1007/978-3-662-12616-5
• Kullback, S. et Leibler, R. A. (1951). On Information and Sufficiency. The Annals of Mathematical Statistics 22(1), 79–86. https://doi.org/10.1214/aoms/1177729694
• Lienen, M., Lüdke, D., Hansen-Palmus, J. et Günnemann, S. (2024). From Zero to Turbulence: Generative Modeling for 3D Flow Simulation. ICLR 2024. arXiv:2306.01776. https://arxiv.org/abs/2306.01776
• Lipman, Y., Chen, R. T. Q., Ben-Hamu, H., Nickel, M. et Le, M. (2023). Flow Matching for Generative Modeling. ICLR 2023. arXiv:2210.02747. https://arxiv.org/abs/2210.02747
• Liu, X., Gong, C. et Liu, Q. (2023). Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow. ICLR 2023. arXiv:2209.03003. https://arxiv.org/abs/2209.03003
• Mun, J. (2026). Multifidelity Computer Model Emulation Via Diffusion Model Steering and Targeted Maximum Likelihood. Préprint arXiv:2608.30096. https://arxiv.org/abs/2608.30096
• Øksendal, B. (2003). Stochastic Differential Equations: An Introduction with Applications, 6e éd. Springer (Universitext). https://doi.org/10.1007/978-3-642-14394-6
• Owen, A. B. (2013). Monte Carlo theory, methods and examples. Livre en ligne. https://artowen.su.domains/mc/
• Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S. et Lakshminarayanan, B. (2021). Normalizing Flows for Probabilistic Modeling and Inference. Journal of Machine Learning Research 22(57), 1–64. arXiv:1912.02762. https://arxiv.org/abs/1912.02762
• Pavliotis, G. A. (2014). Stochastic Processes and Applications. Springer (Texts in Applied Mathematics). https://doi.org/10.1007/978-1-4939-1323-7
• Peyré, G. et Cuturi, M. (2019). Computational Optimal Transport: With Applications to Data Science. Foundations and Trends in Machine Learning 11(5–6), 355–607. https://doi.org/10.1561/2200000073
• Rabin, J., Peyré, G., Delon, J. et Bernot, M. (2012). Wasserstein Barycenter and Its Application to Texture Mixing. Scale Space and Variational Methods in Computer Vision (SSVM 2011), Lecture Notes in Computer Science, 435–446. https://doi.org/10.1007/978-3-642-24785-9_37
• Rezende, D. J. et Mohamed, S. (2015). Variational Inference with Normalizing Flows. ICML 2015. arXiv:1505.05770. https://arxiv.org/abs/1505.05770
• Rice, S. O. (1944). Mathematical Analysis of Random Noise. Bell System Technical Journal 23(3), 282–332. https://doi.org/10.1002/j.1538-7305.1944.tb00874.x
• Rice, S. O. (1945). Mathematical Analysis of Random Noise (suite, parties III et IV). Bell System Technical Journal 24(1), 46–156. https://doi.org/10.1002/j.1538-7305.1945.tb00453.x
• Robert, C. P. et Casella, G. (2004). Monte Carlo Statistical Methods, 2e éd. Springer. https://doi.org/10.1007/978-1-4757-4145-2
• Santambrogio, F. (2015). Optimal Transport for Applied Mathematicians. Birkhäuser. https://doi.org/10.1007/978-3-319-20828-2
• Sejdinovic, D., Sriperumbudur, B., Gretton, A. et Fukumizu, K. (2013). Equivalence of distance-based and RKHS-based statistics in hypothesis testing. The Annals of Statistics 41(5), 2263–2291. https://doi.org/10.1214/13-AOS1140
• Singhal, R., Horvitz, Z., Teehan, R., Ren, M., Yu, Z., McKeown, K. et Ranganath, R. (2025). A General Framework for Inference-time Scaling and Steering of Diffusion Models. ICML 2025, PMLR 267, 55810–55827. https://proceedings.mlr.press/v267/singhal25b.html
• Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S. et Poole, B. (2021). Score-Based Generative Modeling through Stochastic Differential Equations. ICLR 2021. arXiv:2011.13456. https://arxiv.org/abs/2011.13456
• Székely, G. J. et Rizzo, M. L. (2013). Energy statistics: A class of statistics based on distances. Journal of Statistical Planning and Inference 143(8), 1249–1272. https://doi.org/10.1016/j.jspi.2013.03.018
• Tauberschmidt, J., Fellenz, S., Vollmer, S. J. et Duncan, A. B. (2026). Physics-Constrained Fine-Tuning of Flow-Matching Models for Generation and Inverse Problems. ICLR 2026. arXiv:2508.09156. https://arxiv.org/abs/2508.09156
• Tong, A., Fatras, K., Malkin, N., Huguet, G., Zhang, Y., Rector-Brooks, J., Wolf, G. et Bengio, Y. (2024). Improving and generalizing flow-based generative models with minibatch optimal transport. Transactions on Machine Learning Research. arXiv:2302.00482. https://arxiv.org/abs/2302.00482
• Uhlenbeck, G. E. et Ornstein, L. S. (1930). On the Theory of the Brownian Motion. Physical Review 36(5), 823–841. https://doi.org/10.1103/PhysRev.36.823
• Villani, C. (2009). Optimal Transport: Old and New. Springer (Grundlehren der mathematischen Wissenschaften). https://doi.org/10.1007/978-3-540-71050-9
• Wu, L., Trippe, B. L., Naesseth, C. A., Blei, D. M. et Cunningham, J. P. (2023). Practical and Asymptotically Exact Conditional Sampling in Diffusion Models. NeurIPS 2023. arXiv:2306.17775. https://arxiv.org/abs/2306.17775