PB BIBLIO 1952 : L’algorithme de Metropolis

Le projet Manhattan demandant une quantité énorme de calculs, Metropolis se concentra sur la construction de MANIAC (mathematical analyzer numerical integrator and calculator), le premier ordinateur existant au laboratoire de Los Alamos, dont l’article fondateur [MRRTT] décrit l’une des premières utilisations pratiques.

Dynamique moléculaire et physique statistique

De quoi s’agit-il ? Le but est de décrire les propriétés statistiques d’un système composé d’un grand nombre \(N\) de particules en interaction. L’état à un instant donné d’un tel système est décrit par les \(6N\) coordonnées des positions et des vitesses des particules, que l’on peut regrouper dans un vecteur \(x\) de dimension \(6N\). Supposons, pour simplifier, que l’on a discrétisé les positions et les vitesses, de sorte que \(x\) prend ses valeurs dans un ensemble grand mais fini \(\mathcal{X}\) : son cardinal est \(|\mathcal{X}| = y^N\), où \(y\) est le nombre de positions et de vitesses autorisées.

Simuler les trajectoires de toutes les particules est trop coûteux en temps de calcul. Une idée cruciale, due à Edward Teller, est d’utiliser le fait qu’à l’équilibre thermodynamique, les particules suivent une distribution de probabilité, la distribution de Boltzmann–Gibbs. La probabilité de trouver le système dans la configuration \(x \in \mathcal{X}\) est donnée par
\[
p(x) = \frac{1}{Z} \text{e}^{-E(x)/(k_{\text{B}}T)}
\]
où \(E(x)\) est l’énergie de la configuration \(x\), \(T\) est la température, \(k_{\text{B}}\) est la [constante de Boltzmann->https://fr.wikipedia.org/wiki/Constante_de_Boltzmann], et \(Z\) est la constante de normalisation
\[
Z = \sum_{x\in\mathcal{X}} \text{e}^{-E(x)/(k_{\text{B}}T)}
\]
appelée fonction de partition.

On peut alors calculer la valeur moyenne de toute fonction \(F\) de la configuration (comme par exemple son énergie) par la formule
\[
\mathbb{E}(F) = \sum_{x\in\mathcal{X}} F(x) p(x)
= \frac{1}{Z} \sum_{x\in\mathcal{X}} F(x) \text{e}^{-E(x)/(k_{\text{B}}T)}
\]
Calculer cette somme de \(y^N\) termes reste trop coûteux numériquement. Toutefois, en 1946, John von Neumann et Stanislaw Ulam avaient inventé les méthodes de Monte-Carlo. Leur approche consiste à ne pas calculer tous les termes de la somme, mais seulement un échantillon relativement petit de valeurs représentatives, choisies avec probabilités proportionnelles à \(p(x)\).

Monte-Carlo par chaînes de Markov

Il reste toutefois le problème d’échantillonner correctement la mesure de probabilité \(p\), d’autant plus qu’on ne veut pas avoir à calculer précisément la fonction de partition \(Z\). Le problème est que dans la plupart des applications, \(p\) est très loin d’être uniforme : beaucoup de configurations ont une probabilité extrêmement faible, et les choisir rendrait l’algorithme inefficace. Comment faire alors pour échantillonner correctement les configurations ?

L’idée de génie introduite dans l’article [MRRTT] est la suivante. Au lieu de choisir les configurations \(x\) uniformément au hasard, on construit une trajectoire aléatoire \((X_1,X_2,\dots)\) dans l’espace \(\mathcal{X}\) qui tend à visiter les configurations les plus probables. D’une certaine manière, il s’agit d’une approximation de la dynamique de départ des particules, mais bien plus simple à simuler. Concrètement, à l’étape \(n\) de la simulation, on modifie légèrement la configuration \(X_n\) actuelle, en déplaçant une particule ou en modifiant sa vitesse au hasard. Notons \(X’_n\) cette nouvelle configuration. On calcule alors l’énergie \(E(X_n’)\), et on distingue deux cas :
-* si l’énergie de la nouvelle configuration est plus basse que celle de \(X_n\), alors on pose \(X_{n+1} = X’_n\) ;
-* si l’énergie est plus élevée, alors on pose \(X_{n+1} = X’_n\) avec probabilité \(p_n = \text{e}^{-[E(X_n’)-E(X_n)]/(k_{\text{B}}T)}\), et \(X_{n+1} = X_n\) avec probabilité \(1-p_n\).

Mathématiquement parlant, la suite des \(X_n\) forme une chaîne de Markov. On vérifie qu’elle admet bien la distribution de Boltzmann–Gibbs comme mesure invariante ²De plus, la chaîne de Markov est réversible par rapport à la mesure invariante \(p\). Cela signifie que si \(P(x\to x’)\) désigne la probabilité de transition de l’état \(x\) vers l’état \(x’\), alors on a la relation d’équilibre détaillé \(p(x) P(x\to x’) = p(x’) P(x’\to x)\)., et que la loi de \(X_n\) converge vers cette mesure. De plus, le nombre d’itérations nécessaires à calculer la moyenne de \(F\) avec une précision donnée est essentiellement indépendante du nombre de particules (et varie comme l’inverse de la précision requise).

Cet algorithme, ainsi que sa généralisation à des distributions autres que celle de Boltzmann–Gibbs, l’algorithme de Metropolis–Hastings, ont connu un succès gigantesque en dynamique moléculaire. Voici un exemple d’application au modèle d’Ising pour les matériaux ferromagnétiques :

Une description plus détaillée s’en trouve dans cet article sur Images des Mathématiques.

Si Metropolis a bien travaillé avec Ulam sur l’implémentation numérique de méthodes Monte-Carlo, selon le compte-rendu dans [GU] est la suivante. Au lieu de choisir les , le gros du travail conduisant à l’algorithme décrit ici a été fait par Marshall et Arianna Rosenbluth. L’idée décisive d’utiliser une approche de physique statistique et due à Edward Teller, et une première version de code numérique avait été écrite par Augusta Teller. Il est donc quelque peu injuste que seul le nom de Metropolis soit resté associé à cet algorithme fondateur.

Bibliographie

[Gu] J. E. Gubernatis, (2005), Marshall Rosenbluth and the Metropolis algorithm. Phys. Plasmas 12: 057303.

[MRRTT] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, (1953), Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21: 1087-1953.

[Ro] M. N. Rosenbluth, (2003), Genesis of the Monte-Carlo Algorithm for Statistical Mechanics. AIP Conference Proceedings 690, 22.