1952 : L’algorithme de Metropolis

El 11 febrero 2020  - Escrito por  Nils Berglund Ver los comentarios

L’algorithme de Metropolis, qui devrait plutôt s’appeler algorithme de Metropolis-Rosenbluth-Rosenbluth-Teller-Teller, est un outil fondamental pour l’échantillonnage de distributions de probabilité par chaînes de Markov.

Il résulta d’une collaboration, dans les années 1940-50 à Los Alamos, entre Nicholas Metropolis, le couple Arianna Rosenbluth et Marshall Rosenbluth, et le couple Augusta Teller et Edward Teller. Leurs travaux culminèrent dans la publication, en 1953, de l’article Equation of State Calculations by Fast Computing Machines, qui allait révolutionner la dynamique moléculaire.

Les détails de l’histoire de cette collaboration sont longtemps restés obscurs, et sont aujourd’hui encore controversés. Des précisions sont apparues dans les articles [Ro] et [Gu], peu après que Marshall Rosenbluth eut assisté à une conférence célébrant les 50 ans de l’algorithme au Laboratoire national de Los Alamos dans la ville éponyme. Ce laboratoire fut créé en secret en 1943, dans le but de développer la bombe atomique, et les physiciens Edward Teller et Nicholas Metropolis, ainsi que l’informaticienne Augusta H. Teller, le rejoignirent la même année. Edward Teller recruta le physicien des plasmas Marshall Rosenbluth en 1950, qui y rencontrera sa première femme, la physicienne et informaticienne Arianna W. Rosenbluth.

PNG - 57.5 KB
Nicholas Metropolis
PNG - 95.5 KB
Marshall Rosenbluth
JPEG - 25.9 KB
Edward Teller et Augusta Teller

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, 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 [1], 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], 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.

Article édité par Laurent Bartholdi

Notas

[1De 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)$.

Comparte este artículo

Para citar este artículo:

Nils Berglund — «1952 : L’algorithme de Metropolis» — Images des Mathématiques, CNRS, 2020

Comentario sobre el artículo

Dejar un comentario

Foro sólo para inscritos

Para participar en este foro, debe registrarte previamente. Gracias por indicar a continuación el identificador personal que se le ha suministrado. Si no está inscrito/a, debe inscribirse.

Conexióninscribirse¿contraseña olvidada?

Registros

Este artículo es parte del registro «Les 80 ans du CNRS» consulte el registro
La traducción del sitio del francés al castellano se realiza gracias al apoyo de diversas instituciones de matemáticas de América Latina.