Chaînes de Markov
Une chaîne de Markov, c’est un jeton qui saute de case en case au hasard, avec une seule règle : pour décider où aller, il ne regarde que la case où il se trouve. Peu importe le chemin qui l’y a mené. Le futur ne dépend du passé que par le présent.
Cette amnésie paraît être une limitation sévère. C’est au contraire ce qui rend ces objets si maniables, et si répandus : la météo du lendemain, une partie de pile ou face, un jeu de cartes qu’on bat, un marcheur égaré, un moteur de recherche, ou les milliards de spins d’un aimant qu’on simule sur ordinateur. Autant d’exemples de la polyvalence des chaînes de Markov.
En 1913, Andreï Markov publie la première application de ses chaînes, sur un terrain inattendu : le roman en vers Eugène Onéguine de Pouchkine. Il dispose les 20 000 premières lettres du texte (sans espaces ni ponctuation) en tableaux de cent lettres, puis compte les voyelles, les consonnes et les paires de lettres consécutives. Sur les 19 999 paires, les voyelles représentent 43 % des lettres, mais une voyelle n'est suivie d'une autre voyelle que 1 104 fois sur 8 638 (13 %), alors qu'une consonne est suivie d'une voyelle 7 534 fois sur 11 361 (66 %). Les lettres ne sont pas indépendantes : chacune « se souvient » de la précédente. Markov essaie même une mémoire de deux lettres, qui rend encore mieux compte des fluctuations observées dans le texte.
La théorie, elle, datait de 1906 et était née d'une querelle. Le mathématicien Pavel Nekrassov soutenait que la loi des grands nombres (les moyennes se stabilisent) exigeait des événements indépendants, et il en tirait même des conclusions théologiques sur le libre arbitre. Markov construisit ses chaînes pour lui démontrer que des événements enchaînés obéissent eux aussi à la loi des grands nombres ; Onéguine en fournit, sept ans plus tard, l'illustration concrète.
Une promenade au hasard sur un graphe
Voici un modèle (très) simplifié de la météo : trois états, soleil, nuages et pluie. Chaque flèche porte la probabilité de passer d’un état à l’autre d’un jour au suivant. Après une journée ensoleillée, il fait encore beau six fois sur dix ; après une journée de pluie, il pleut encore quatre fois sur dix. Le jeton rouge parcourt le graphe en tirant au sort chaque nouvelle flèche.
Le jeton suit les flèches au hasard, avec les probabilités indiquées. En bas, la proportion des jours passés dans chaque état : elle fluctue d'abord, puis se cale sur les tirets, quelle que soit l'histoire du jeton. Accélérez pour voir la convergence.
Toute l’information tient dans un tableau, la matrice de transition, dont la case $(i,j)$ contient la probabilité d’aller de $i$ vers $j$ :
$\displaystyle P=\begin{pmatrix}0{,}6 & 0{,}3 & 0{,}1\\ 0{,}3 & 0{,}4 & 0{,}3\\ 0{,}2 & 0{,}4 & 0{,}4\end{pmatrix}$
Chaque ligne est une loi de probabilité (ses termes sont positifs et de somme $1$) : c’est le « dé » qu’on lance quand le jeton se trouve dans l’état correspondant.
Chaîne de Markov. Une suite d'états aléatoires $X_0, X_1, X_2,\dots$ à valeurs dans un ensemble fini $E$ est une chaîne de Markov de matrice de transition $P=(p_{ij})$ si, pour tout instant $t$ et tous états,
$\mathbb{P}(X_{t+1}=j \mid X_t=i,\ X_{t-1}=i_{t-1},\dots,X_0=i_0)=\mathbb{P}(X_{t+1}=j\mid X_t=i)=p_{ij}.$
La matrice $P$ est dite stochastique : ses coefficients sont positifs et chaque ligne a pour somme $1$.
Ce qui rend le formalisme si agréable, c’est que suivre la loi de probabilité revient à multiplier des matrices. Notons $\mu_t$ la loi de $X_t$, écrite comme un vecteur ligne (par exemple $\mu_0=(1,0,0)$ s’il fait beau le jour $0$). La formule des probabilités totales donne immédiatement :
$\displaystyle \mu_{t+1}=\mu_t\,P \qquad\text{et donc}\qquad \mu_t=\mu_0\,P^t.$
La case $(i,j)$ de $P^t$ est la probabilité d’aller de $i$ à $j$ en exactement $t$ pas. Toute la dynamique est contenue dans les puissances d’une matrice.
En 1948, dans l'article qui fonde la théorie de l'information, Claude Shannon reprend l'idée de Markov pour fabriquer du texte. On choisit chaque lettre au hasard en ne regardant que les $k$ précédentes, avec les fréquences observées dans un vrai texte. C'est une chaîne de Markov dont les états sont les suites de $k$ lettres. Voici la recette, appliquée aux pages de ce site.
1. Le texte de référence. On prend les pages du site et l'on n'en garde que la prose : titres, figures, tableaux, programmes et paragraphes contenant des formules sont écartés. Mis bout à bout, les paragraphes restants forment un seul long texte d'environ $390\,000$ caractères ($62\,000$ mots), écrit avec $123$ caractères différents (minuscules, majuscules, lettres accentuées, chiffres, ponctuation, espace).
2. La matrice de transition. L'état de la chaîne est la suite des $k$ derniers caractères écrits. Pour connaître la ligne de la matrice associée à un état, on relève toutes ses apparitions dans le texte et l'on compte le caractère qui le suit. Par exemple, avec $k=2$, l'état « qu » apparaît $4036$ fois ; il est suivi de « e » dans $61\,\%$ des cas, de « i » dans $14\,\%$, d'une apostrophe dans $12\,\%$ et de « a » dans $10\,\%$. Ces fréquences sont les probabilités de transition depuis « qu ».
3. Le tirage. On part du début d'une phrase prise au hasard dans le texte. Puis, à chaque pas, on tire le caractère suivant selon la ligne de l'état courant, on l'écrit, et l'état se décale d'un cran (on oublie le caractère le plus ancien). On s'arrête au premier point final après $150$ caractères.
En pratique, on ne construit jamais la matrice, qui serait énorme et presque vide. Il suffit de remarquer que tirer un caractère « selon les fréquences observées » revient à choisir au hasard l'une des apparitions de l'état dans le texte, et à recopier le caractère qui la suit. Shannon procédait ainsi à la main : il ouvrait un livre au hasard, lisait jusqu'à rencontrer les dernières lettres écrites, et notait la suivante. Le bouton ci-dessous fait exactement cela dans votre navigateur. Pour une mémoire comptée en mots, la recette est la même en remplaçant « caractère » par « mot » (les mots étant simplement les morceaux de texte séparés par des espaces, ponctuation comprise).
Deux lettres de mémoire :
Loet un touvermen les la mitégé ter l'un grant Want ormit a blendande enté de qu'onsen, parate le ille l'averation. Sarration, ent plus la maitreme ra bay façon.
pas forcés : 1 %
Six lettres de mémoire :
L'idée durant plus grand des grandeur de matière optique de la particle science. En juin 1696, il a une symétrique. La rétrocausalité. Le théorie est orthogonality Fléau de mémoire n'est de masse.
pas forcés : 45 %
Deux mots de mémoire :
On va d'abord chercher à mesurer tous les trois propriétés réflexivité, symétrie et transitivité. Mais on se heurte à la syntaxe. Pour le rendre plus scientifique, on peut prouver une formule à droite est une ligne droite infinie.
pas forcés : 50 %
Le français émerge à mesure que la mémoire s'allonge, sans que la machine « comprenne » quoi que ce soit. Mais il y a un revers, que mesure le pourcentage de pas forcés affiché sous chaque texte : la proportion de pas où l'état courant n'avait, dans tout le texte de référence, qu'une seule suite possible. Avec deux lettres, la chaîne invente vraiment : on ne rencontre que $1514$ états différents, et les plus courants ont chacun de nombreuses suites possibles. Avec six lettres, on compte $95\,637$ états différents, dont $77\,\%$ n'ont qu'une suite ; avec deux mots, $86\,\%$. La chaîne passe alors une bonne part de son temps à recopier des fragments du site, et ne bifurque que de loin en loin. Allonger encore la mémoire ne donnerait que des copies conformes : le nombre d'états possibles croît exponentiellement avec $k$, et aucun texte n'est assez long pour les faire tous apparaître plusieurs fois.
Les modèles de langage actuels prolongent cette idée (prédire le mot suivant à partir de ceux qui précèdent) en contournant précisément cet obstacle : au lieu de consulter un tableau de fréquences, ils calculent les probabilités de la suite avec un réseau de neurones, capable de traiter de façon semblable des contextes qui se ressemblent sans avoir jamais été vus tels quels.
Le point fixe de la danse des probabilités
Dans l’animation, l’histogramme se stabilise. Ce n’est pas un hasard : la loi $\mu_t$ elle-même converge. Avec trois états, une loi de probabilité $(a,b,c)$ est un point d’un triangle, le simplexe (chaque sommet est une certitude, le centre est l’incertitude maximale, voir la page Simplexes et probabilités). Multiplier par $P$ déplace ce point. Itérons depuis différents départs :
Chaque trajectoire suit la loi du temps qu'il fait, jour après jour, en partant d'une certitude (un sommet) ou d'une loi quelconque (en rouge), tirée au hasard à chaque lancement. Toutes convergent vers le même point $\pi$ : la chaîne oublie d'où elle est partie. Cliquez dans le triangle pour lancer une nouvelle trajectoire depuis le point de votre choix.
Le point d’arrivée $\pi$ ne bouge plus quand on le multiplie par $P$ : c’est un point fixe, la loi stationnaire (ou invariante). Ici, $\pi=\left(\frac{24}{61},\frac{22}{61},\frac{15}{61}\right)\approx(0{,}39\;;\;0{,}36\;;\;0{,}25)$ : c’est la hauteur des tirets de l’histogramme. Sur le long terme, il fait beau 39 % du temps, quel que soit le temps qu’il faisait au départ.
Comment trouver $\pi$ sans attendre la convergence ? En écrivant précisément qu’elle ne bouge pas : $\pi=\pi P$. Lue état par état, cette égalité est un bilan de flux :
$\displaystyle \pi(j)=\sum_{i}\pi(i)\,p_{ij}\qquad\text{pour chaque état } j.$
La probabilité d’être en $j$ le lendemain est la somme des probabilités d’y arriver depuis chacun des états ; à l’équilibre, ce qui entre dans chaque état compense exactement ce qui en sort. Il suffit d’y ajouter la condition $\sum_j \pi(j)=1$ pour obtenir un système linéaire, que l’on résout comme n’importe quel autre.
le calcul pour la météo
Notons $\pi=(\pi_{\text{s}},\pi_{\text{n}},\pi_{\text{p}})$ les probabilités du soleil, des nuages et de la pluie. Chaque équation se lit sur une colonne de $P$ : pour avoir du soleil demain, il faut du soleil aujourd’hui puis rester (probabilité $0{,}6$), ou des nuages puis une éclaircie ($0{,}3$), ou de la pluie puis une éclaircie ($0{,}2$). D’où le système :
$\begin{cases}\pi_{\text{s}}=0{,}6\,\pi_{\text{s}}+0{,}3\,\pi_{\text{n}}+0{,}2\,\pi_{\text{p}}\\ \pi_{\text{n}}=0{,}3\,\pi_{\text{s}}+0{,}4\,\pi_{\text{n}}+0{,}4\,\pi_{\text{p}}\\ \pi_{\text{p}}=0{,}1\,\pi_{\text{s}}+0{,}3\,\pi_{\text{n}}+0{,}4\,\pi_{\text{p}}\\ \pi_{\text{s}}+\pi_{\text{n}}+\pi_{\text{p}}=1\end{cases}$
Les trois premières équations ne sont pas indépendantes : en les additionnant, on obtient une identité, puisque chaque ligne de $P$ a pour somme $1$. Elles ne déterminent $\pi$ qu’à un facteur près, et c’est la quatrième, la normalisation, qui fixe ce facteur. En multipliant par $10$ et en regroupant, les trois premières deviennent :
$\begin{cases}4\,\pi_{\text{s}}=3\,\pi_{\text{n}}+2\,\pi_{\text{p}}\\ 6\,\pi_{\text{n}}=3\,\pi_{\text{s}}+4\,\pi_{\text{p}}\\ 6\,\pi_{\text{p}}=\pi_{\text{s}}+3\,\pi_{\text{n}}\end{cases}$
La troisième donne $\pi_{\text{s}}=6\,\pi_{\text{p}}-3\,\pi_{\text{n}}$. Reportée dans la première : $24\,\pi_{\text{p}}-12\,\pi_{\text{n}}=3\,\pi_{\text{n}}+2\,\pi_{\text{p}}$, soit $\pi_{\text{n}}=\frac{22}{15}\,\pi_{\text{p}}$, puis $\pi_{\text{s}}=\frac{24}{15}\,\pi_{\text{p}}$ (la deuxième équation est alors automatiquement vérifiée). La normalisation donne enfin $\frac{24+22+15}{15}\,\pi_{\text{p}}=1$, d’où
$\pi=\left(\dfrac{24}{61},\ \dfrac{22}{61},\ \dfrac{15}{61}\right).$
Dans le langage de l’algèbre linéaire, $\pi$ est un vecteur propre à gauche de $P$ pour la valeur propre $1$ (autrement dit un vecteur propre de la transposée $P^{\mathsf T}$). Cette valeur propre existe toujours : le vecteur colonne $(1,1,1)$ vérifie $P\,(1,1,1)^{\mathsf T}=(1,1,1)^{\mathsf T}$, justement parce que chaque ligne a pour somme $1$. Le théorème de Perron-Frobenius garantit de plus que, pour une chaîne irréductible, ce vecteur propre peut être choisi à coefficients strictement positifs, et qu’il est unique une fois normalisé.
Pour une très grande chaîne (des milliards de pages web, par exemple), on ne résout pas le système : on part d’une loi quelconque et on la multiplie par $P$ encore et encore, exactement comme dans le triangle ci-dessus. C’est la méthode de la puissance, et le théorème de convergence assure qu’elle aboutit.
Théorème de convergence
Si la chaîne est irréductible (on peut aller de tout état à tout autre) et apériodique (elle n’est pas condamnée à revenir sur ses pas selon un rythme fixe), alors il existe une unique loi $\pi$ telle que $\pi P=\pi$, et pour toute loi de départ $\mu_0$ :
$\mu_0P^t\xrightarrow[t\to\infty]{}\pi.$
De plus, la proportion du temps passé dans chaque état $x$ tend vers $\pi(x)$ (loi des grands nombres de Markov).
Chacune des deux hypothèses exclut une façon de ne pas converger.
À gauche, le jeton alterne à coup sûr : partie de A, sa loi vaut $(1,0)$, $(0,1)$, $(1,0)$… et ne converge jamais (pourtant $\left(\frac12,\frac12\right)$ est bien stationnaire). À droite, une fois dans l'un des deux mondes, on n'en sort plus : la limite dépend du point de départ.
pourquoi ça converge : le cas de la grenouille
Le livre de Levin et Peres ouvre sur une grenouille qui vit sur deux nénuphars, E et O. Chaque matin, elle lance la pièce posée sur son nénuphar et saute sur l’autre si elle obtient pile : avec probabilité $p$ depuis E, $q$ depuis O. La matrice est $P=\begin{pmatrix}1-p&p\\ q&1-q\end{pmatrix}$ et la loi stationnaire vaut $\pi=\left(\frac{q}{p+q},\frac{p}{p+q}\right)$.
Notons $\Delta_t=\mu_t(\text{E})-\frac{q}{p+q}$ l’écart à l’équilibre. Comme $\mu_{t+1}(\text{E})=(1-p)\mu_t(\text{E})+q\big(1-\mu_t(\text{E})\big)$, un calcul d’une ligne donne $\Delta_{t+1}=(1-p-q)\Delta_t$, donc :
$\Delta_t=(1-p-q)^t\,\Delta_0.$
L’écart fond comme une suite géométrique de raison $1-p-q$, qui est la seconde valeur propre de $P$. Il ne fond pas quand $|1-p-q|=1$ : soit $p=q=0$ (la grenouille ne bouge jamais, la chaîne est réductible), soit $p=q=1$ (elle saute tous les jours, la chaîne est périodique). On retrouve exactement les deux obstacles.
Pour la météo, c’est pareil avec trois états : les valeurs propres de $P$ sont $1$, $\frac{2+\sqrt3}{10}\approx0{,}37$ et $\frac{2-\sqrt3}{10}\approx0{,}03$. La première correspond à $\pi$, les deux autres mesurent la vitesse à laquelle le reste s’efface : l’écart à l’équilibre est divisé par près de $3$ chaque jour, d’où la convergence en quelques pas visible dans le triangle.
En 1998, Sergey Brin et Larry Page cherchent à classer les pages du web. Leur idée : imaginer un internaute qui clique au hasard sur l'un des liens de la page où il se trouve, et qui, de temps en temps (15 % du temps), saute sur une page quelconque. C'est une chaîne de Markov dont les états sont les pages. Le score d'une page, son PageRank, est la proportion du temps que l'internaute y passe à long terme : sa loi stationnaire.
Six pages et leurs liens ; la taille de chaque disque est proportionnelle à son PageRank. Une page est importante si des pages importantes pointent vers elle : la définition est circulaire, et c'est précisément l'équation $\pi=\pi P$ qui la rend cohérente.
Le saut aléatoire n'est pas qu'un détail : il rend la chaîne irréductible et apériodique, ce qui garantit l'existence, l'unicité et la convergence du classement, et il évite que l'internaute reste piégé dans un recoin du web sans liens sortants.
Converger, oui, mais en combien de temps ? Battre un jeu de cartes, c'est faire avancer une chaîne de Markov dont les états sont les $52!\approx8\times10^{67}$ ordres possibles du jeu ; sa loi stationnaire est la loi uniforme, le mélange parfait. Pour le battement à l'américaine (on coupe le jeu en deux et on entrelace les deux paquets), Dave Bayer et Persi Diaconis ont calculé en 1992 la distance exacte au mélange parfait après chaque battement.
Un jeu de 52 cartes, colorées selon leur position de départ, après 0 à 12 battements. À droite, la distance exacte au mélange parfait (la plus grande erreur possible sur la probabilité d'un événement). Elle reste collée à $1$ pendant cinq battements, puis s'effondre en trois ou quatre coups.
Ce décrochage brutal porte un nom, le phénomène de seuil (cutoff) : rien ne se passe, puis tout se passe d'un coup. D'où le célèbre « sept battements suffisent », à nuancer : après sept battements, la distance vaut encore $0{,}33$. Peter Doyle a inventé une réussite qu'on gagne avec probabilité $\frac12$ avec un jeu parfaitement mélangé, mais avec probabilité $0{,}80$ après sept battements d'un jeu initialement trié. Pour descendre à 1 % d'erreur, comparable à l'avantage d'un casino, il faut plutôt onze ou douze battements.
Le jeu de Penney
Deux joueurs choisissent chacun une séquence de trois lancers de pièce, par exemple FFP (face, face, pile) et PFF. On lance la pièce jusqu’à ce que l’une des deux séquences apparaisse : son propriétaire gagne. Toutes les séquences de trois lancers sont équiprobables, alors le jeu semble équitable. Il ne l’est pas du tout.
Le tableau complet des probabilités de victoire est saisissant :
Probabilité que la séquence du joueur 2 (en ligne) sorte avant celle du joueur 1 (en colonne), calculée exactement. Les cases encadrées sont la meilleure réponse du joueur 2 : quelle que soit la séquence adverse, il existe une réponse qui gagne au moins deux fois sur trois.
Il n’y a donc pas de meilleure séquence : celui qui choisit en second gagne. Les meilleures réponses forment même un cycle, comme à pierre-feuille-ciseaux : FFP est battue par PFF, elle-même battue par PPF, battue par FPP, qui est battue par… FFP.
Les meilleures réponses, avec leur probabilité de victoire. La recette, due à John Conway : prendre le contraire du deuxième lancer adverse, puis y accoler les deux premiers lancers adverses (contre FFP : P suivi de FF, soit PFF). Sous chaque séquence, le nombre moyen de lancers nécessaires pour la voir sortir.
Essayez de battre l’ordinateur, qui applique la recette de Conway :
Choisissez votre séquence, puis lancez des parties. L'ordinateur répond toujours avec la recette de Conway ; sur 100 parties, la barre ne ment pas.
Une course entre deux chaînes de Markov
Pour comprendre, suivons la partie comme un jeton sur un graphe. L’état est la fin de la suite des lancers qui peut encore servir à l’un des joueurs. Voici la course FFP contre PFF :
Chaque flèche a une probabilité $\frac12$. Les états FFP et PFF sont absorbants : une fois atteints, la partie est finie.
Le dessin dit tout. Dès qu’un P sort, le jeton entre dans la zone rose, et il n’en ressortira que par PFF : tout F qui suit un P prépare PFF, et FF ne peut plus apparaître sans être précédé d’un P. La seule chance de FFP est donc que les deux premiers lancers soient FF, soit une chance sur quatre. PFF gagne trois fois sur quatre. Le secret de la réponse de Conway est là : se placer juste devant la séquence adverse, pour qu’elle serve de tremplin.
Le paradoxe des temps d’attente
On pourrait croire qu’une séquence gagne parce qu’elle sort plus vite en moyenne. Il n’en est rien. Il faut en moyenne 14 lancers pour voir FFF et seulement 8 pour FFP ; pourtant, face à face, elles ont exactement une chance sur deux chacune (dès que FF sort, le lancer suivant tranche). Inversement, FFP et PFF demandent toutes deux 8 lancers en moyenne, mais PFF gagne trois fois sur quatre. Le temps d’attente moyen ne dit rien de l’issue d’une course : c’est la façon dont les deux séquences se recouvrent qui décide.
D’où viennent ces nombres ? Une séquence qui « se recouvre » elle-même est lente à venir : un échec à la dernière lettre de FFF ne laisse aucun acquis, alors qu’un échec de FFP (on tire F au lieu de P) laisse FF, prêt à servir.
calculer le temps d'attente avec un casino (l'argument des parieurs)
Voici un argument d’une élégance rare, présenté dans le livre de Levin et Peres. Installons un casino équitable où l’on parie sur chaque lancer : une mise de $k$ euros sur le bon côté rapporte $2k$, sinon elle est perdue.
Avant chaque lancer, un nouveau parieur arrive avec 1 euro et suit la séquence attendue, disons FPF : il mise tout sur F ; s’il gagne, il remise ses 2 euros sur P ; s’il gagne encore, il remise ses 4 euros sur F. Au premier échec, il repart les poches vides.
Arrêtons tout quand FPF sort, au lancer $\tau$. À ce moment, $\tau$ parieurs ont misé chacun 1 euro. Qui repart avec de l’argent ? Celui arrivé trois lancers plus tôt, qui a vu sortir FPF en entier : il emporte $8$ euros. Celui arrivé au dernier lancer a misé sur F et gagné : $2$ euros. Celui de l’avant-dernier lancer avait misé sur F et vu sortir P : rien. Les gains totaux valent donc $8+2=10$, quoi qu’il arrive.
Le casino étant équitable, les joueurs ne peuvent rien gagner en moyenne : $\mathbb{E}(\text{gains}-\text{mises})=0$, soit $10-\mathbb{E}(\tau)=0$. D’où $\mathbb{E}(\tau)=10$.
En général, le temps d’attente moyen d’une séquence est la somme des $2^k$ pour toutes les longueurs $k$ telles que les $k$ dernières lettres de la séquence coïncident avec ses $k$ premières : $8$ pour FFP (aucun recouvrement), $8+2=10$ pour FPF, $8+4+2=14$ pour FFF. Avec deux familles de parieurs, une par joueur, le même argument donne la formule de Conway pour les probabilités de victoire.
la méthode générale : la matrice fondamentale d'une chaîne absorbante
Rangeons d’abord les états transitoires (départ, F, P, FF…), puis les états absorbants (les deux séquences gagnantes). La matrice de transition prend la forme par blocs
$P=\begin{pmatrix}Q & R\\ 0 & I\end{pmatrix},$
où $Q$ contient les transitions entre états transitoires et $R$ celles qui mènent à un état absorbant. Par récurrence,
$P^n=\begin{pmatrix}Q^n & (I+Q+\dots+Q^{n-1})R\\ 0 & I\end{pmatrix}.$
Comme le jeton finit forcément absorbé, $Q^n\to0$. La série $I+Q+Q^2+\dots$ converge alors vers la matrice fondamentale :
$N=(I-Q)^{-1}.$
Son coefficient $N_{ij}$ est le nombre moyen de passages par l’état $j$ en partant de $i$ (chaque terme $(Q^n)_{ij}$ est la probabilité d’y être au temps $n$). On en tire tout :
- la somme de la ligne $i$ de $N$ est la durée moyenne de la partie en partant de $i$ ;
- la matrice $B=NR$ donne les probabilités d'absorption : $ B_{ i k } $ est la probabilité de finir dans l'état gagnant $k$ en partant de $i$.
Revenir chez soi
Lâchons maintenant un marcheur sur un réseau infini : à chaque pas, il choisit au hasard l’une des directions possibles (gauche ou droite sur une droite, quatre directions dans le plan, six dans l’espace). C’est une chaîne de Markov, avec une infinité d’états. Question simple : reviendra-t-il un jour à son point de départ ?
Proportion de 20 000 marcheurs simulés déjà revenus à l'origine après un nombre de pas donné. Sur une droite et dans le plan, la courbe monte vers 100 % (lentement dans le plan). Dans l'espace, elle plafonne vers 34 %.
C’est le théorème de Pólya (1921) : en dimension 1 et 2, le marcheur revient presque sûrement à son point de départ, et même une infinité de fois (la marche est dite récurrente) ; en dimension 3 et au-delà, il a une probabilité strictement positive (environ 66 % en dimension 3) de ne jamais revenir (la marche est transiente). Shizuo Kakutani l’a résumé d’une formule : un homme ivre finit toujours par retrouver sa maison, un oiseau ivre peut se perdre à jamais.
pourquoi la dimension change tout
L’argument se déroule en quatre temps.
1. Revenir à coup sûr, c’est revenir une infinité de fois. Notons $r$ la probabilité que le marcheur revienne au moins une fois à son point de départ. S’il y revient, il se retrouve exactement dans la situation initiale, sans mémoire de ce qui précède (c’est la propriété de Markov) : tout recommence, et il reviendra encore avec la même probabilité $r$.
- Si $r=1$, il revient, puis revient encore, et ainsi de suite : une infinité de visites.
- Si $r<1$, chaque retour est comme un tirage à pile ou face qu'il finit par perdre. La probabilité de revenir au moins $k$ fois vaut $r^k$, et le nombre moyen de retours vaut $r+r^2+r^3+\dots=\frac{r}{1-r}$, qui est fini.
La question « revient-il à coup sûr ? » équivaut donc à la question « le nombre moyen de visites à l’origine est-il infini ? ».
2. Ce nombre moyen est une somme. Notons $p_t$ la probabilité que le marcheur soit à l’origine après $t$ pas. Le nombre de visites s’obtient en comptant $1$ à chaque instant où il s’y trouve ; sa moyenne est donc simplement la somme de ces probabilités :
$\displaystyle \text{nombre moyen de visites}=p_1+p_2+p_3+\dots=\sum_{t\geqslant1}p_t .$
Tout se ramène à une question d’analyse : cette série diverge-t-elle ?
3. La probabilité d’être à l’origine décroît comme $t^{-d/2}$. Après $t$ pas, le marcheur ne s’est éloigné que d’environ $\sqrt t$ (c’est la loi de la diffusion : les pas se compensent en partie). En dimension $d$, il se trouve donc quelque part dans une région de rayon $\sqrt t$, qui contient environ $(\sqrt t)^d=t^{d/2}$ sites, tous à peu près aussi probables. La probabilité d’être sur l’un d’eux en particulier, l’origine, vaut donc environ
$\displaystyle p_t\approx\frac{1}{t^{d/2}} .$
4. Il reste à sommer. Une série $\sum\frac{1}{t^{\alpha}}$ converge si et seulement si $\alpha>1$ (critère de Riemann). Ici $\alpha=\frac d2$ : la frontière tombe exactement sur $d=2$.
| dimension | $p_t$ environ | la somme | verdict |
|---|---|---|---|
| $1$ (droite) | $\frac{1}{\sqrt t}$ | diverge vite | retour certain, et fréquent |
| $2$ (plan) | $\frac{1}{t}$ | diverge de justesse, comme $\ln t$ | retour certain, mais très tardif |
| $3$ (espace) | $\frac{1}{t\sqrt t}$ | converge | retour incertain |
La dimension $2$ est le cas limite : la série harmonique diverge, mais avec une lenteur extrême, ce qui explique que la courbe bleue monte si paresseusement vers $100\,\%$.
L’image à retenir. En $t$ pas, le marcheur visite au plus $t$ sites, alors que la région où il se promène en contient $t^{d/2}$. Sur une droite, cela fait $\sqrt t$ sites pour $t$ pas : il repasse partout, de nombreuses fois. Dans le plan, $t$ sites pour $t$ pas : c’est tout juste suffisant. Dans l’espace, $t\sqrt t$ sites pour $t$ pas : il y a trop de place, l’essentiel de la région n’est jamais visité, et l’origine a de bonnes chances d’en faire partie.
Revenir est certain sur une droite, mais combien de temps faut-il attendre ? Pour une chaîne finie, la réponse tient en une formule d’une simplicité désarmante, due à Mark Kac :
Temps de retour
Pour une chaîne irréductible de loi stationnaire $\pi$, le temps moyen pour revenir à un état $x$ en partant de $x$ vaut
$\mathbb{E}_x\left(\tau_x^+\right)=\dfrac{1}{\pi(x)}.$
L’intuition est immédiate : si le jeton passe une fraction $\pi(x)$ de son temps en $x$, ses visites sont espacées en moyenne de $1/\pi(x)$ pas.
Cette formule a résolu l’un des grands paradoxes de la physique. En 1907, Paul et Tatiana Ehrenfest proposent un modèle jouet de gaz : $N$ particules réparties entre deux compartiments ; à chaque instant, une particule tirée au hasard change de compartiment. Toutes au départ à gauche, elles se répartissent vite moitié-moitié, comme un gaz qui se détend.
L'urne d'Ehrenfest avec 100 particules. L'équilibre est atteint en quelques centaines de pas, puis le système fluctue autour de 50, selon la loi stationnaire (binomiale) dessinée à droite.
Or Poincaré avait démontré qu’un système mécanique isolé revient toujours, aussi près qu’on veut, de son état initial. Comment concilier ce retour obligatoire avec l’irréversibilité observée, le gaz qui ne se recomprime jamais tout seul dans un coin ? Le modèle d’Ehrenfest donne la réponse, chiffrée. Le retour à l’état « tout à gauche » est certain, mais la loi stationnaire lui donne la probabilité $\pi=\frac{1}{2^N}$ : d’après Kac, il faut attendre en moyenne $2^N$ pas. Pour 100 particules, c’est $2^{100}\approx10^{30}$ pas : à un million de pas par seconde, des millions de fois l’âge de l’Univers. Pour un vrai gaz, avec $N\approx10^{23}$, le nombre défie toute image. Les deux points de vue sont justes : le retour est certain, mais il n’arrive jamais à notre échelle. La flèche du temps est une affaire de nombres gigantesques, pas une impossibilité (voir la page Entropie).
Échantillonner l’impossible : le Monte-Carlo par chaînes de Markov
Jusqu’ici, la chaîne était donnée et l’on cherchait sa loi stationnaire. Renversons la question : on se donne une loi $\pi$, et l’on fabrique une chaîne dont elle est la loi stationnaire. Il suffira alors de laisser tourner la chaîne assez longtemps pour obtenir des tirages selon $\pi$. C’est le Monte-Carlo par chaînes de Markov (MCMC), sans doute l’algorithme le plus utilisé de la physique numérique.
Pourquoi un tel détour ? Parce qu’en physique statistique, la loi à échantillonner est celle de Boltzmann : à la température $T$, une configuration $x$ d’énergie $E(x)$ a la probabilité
$\displaystyle \pi(x)=\frac{\mathrm{e}^{-E(x)/kT}}{Z},\qquad Z=\sum_y \mathrm{e}^{-E(y)/kT}.$
Le numérateur se calcule sans peine, mais la constante $Z$ (la fonction de partition) est une somme sur toutes les configurations : $2^{10\,000}$ termes pour un petit aimant de $100\times100$ spins. Impossible à calculer, et donc impossible de tirer directement selon $\pi$.
Ce que l’on cherche vraiment : des moyennes
En pratique, on ne veut pas connaître $\pi$ configuration par configuration (que ferait-on d’une liste de $2^{10\,000}$ nombres ?), mais des moyennes : l’énergie moyenne, l’aimantation moyenne, leurs fluctuations. Pour une grandeur $A$,
$\displaystyle \langle A\rangle=\sum_x A(x)\,\pi(x).$
Première idée, le Monte-Carlo naïf : tirer des configurations au hasard, uniformément, et faire la moyenne de $A$ en pondérant chacune par son facteur $\mathrm{e}^{-E/kT}$. C’est un échec complet. Presque toutes les configurations tirées ainsi sont du pur désordre, dont le poids de Boltzmann est infime à basse température ; les configurations qui comptent (presque ordonnées) sont une proportion si dérisoire de l’ensemble qu’on ne les tirerait jamais, même en un temps égal à l’âge de l’Univers. On passerait son temps à additionner des zéros.
Il faut donc tirer directement les configurations importantes, avec la fréquence $\pi(x)$ (c’est l’échantillonnage préférentiel). La moyenne devient alors une simple moyenne arithmétique. Et c’est là qu’intervient la chaîne de Markov : si $\pi$ est sa loi stationnaire, le théorème ergodique (la loi des grands nombres des chaînes de Markov, celle-là même que Markov cherchait en 1906) assure que la moyenne le long d’une seule trajectoire converge vers la moyenne selon $\pi$ :
$\displaystyle \frac{1}{n}\sum_{t=1}^{n}A(X_t)\;\xrightarrow[n\to\infty]{}\;\langle A\rangle .$
La chaîne ne visite qu’une fraction ridicule des configurations, mais elle visite les bonnes, dans les bonnes proportions. Reste à la construire.
La règle de Metropolis
L’astuce tient en une règle : depuis la configuration $x$, on propose une petite modification $y$ (retourner un spin, déplacer une particule), puis on l’accepte
- toujours si elle fait baisser l'énergie ;
- avec la probabilité $\mathrm{e}^{-\Delta E/kT}$ si elle l'augmente de $\Delta E$ ;
et sinon on reste sur place. Le rapport $\frac{\pi(y)}{\pi(x)}=\mathrm{e}^{-\Delta E/kT}$ ne fait intervenir que la différence d’énergie : la constante $Z$ a disparu.
Trois précisions, qui sont autant de pièges classiques.
- Un refus compte. Si la modification est refusée, la chaîne reste en $x$, et cet état est compté une fois de plus dans les moyennes. C'est précisément ce « surplace » qui donne aux états de basse énergie leur poids supplémentaire. Jeter les refus fausserait tout.
- La proposition doit être symétrique : proposer $y$ depuis $x$ aussi souvent que $x$ depuis $y$. (Sinon, il faut corriger la probabilité d'acceptation : c'est la variante de Hastings.)
- Les propositions doivent permettre d'aller partout (irréductibilité) : en retournant les spins un par un, on peut relier deux configurations quelconques.
Pourquoi de petites modifications ? Parce qu’une grande modification tirée au hasard ferait presque toujours grimper l’énergie de beaucoup, et serait presque toujours refusée. Les petits pas sont acceptés souvent ; le prix à payer est que deux configurations successives se ressemblent énormément.
En haut, la règle sur un paysage d'énergie à deux vallées : la bille descend toujours, mais remonte parfois, d'autant plus facilement que la température est élevée. En bas, après 400 000 pas, l'histogramme des positions visitées colle à la loi de Boltzmann.
pourquoi la loi stationnaire est bien celle de Boltzmann : l'équilibre détaillé
Une loi $\pi$ est stationnaire dès que le « trafic » entre deux états est le même dans les deux sens :
$\pi(x)\,P(x,y)=\pi(y)\,P(y,x)\quad\text{pour tous } x, y.$
C’est l’équilibre détaillé (la chaîne est dite réversible). En sommant sur $x$, on obtient $\sum_x\pi(x)P(x,y)=\pi(y)\sum_x P(y,x)=\pi(y)$, c’est-à-dire $\pi P=\pi$.
Vérifions-le pour Metropolis, avec une proposition symétrique (on propose $y$ depuis $x$ aussi souvent que $x$ depuis $y$). Supposons $E(y)>E(x)$ : alors $P(x,y)=\text{(proposition)}\times\mathrm{e}^{-(E(y)-E(x))/kT}$ et $P(y,x)=\text{(proposition)}\times1$. D’où
$\pi(x)P(x,y)=\dfrac{\mathrm{e}^{-E(x)/kT}}{Z}\,\mathrm{e}^{-(E(y)-E(x))/kT}\times\text{(prop.)}=\dfrac{\mathrm{e}^{-E(y)/kT}}{Z}\times\text{(prop.)}=\pi(y)P(y,x).$
L’équilibre détaillé est satisfait : la loi de Boltzmann est stationnaire, et si la chaîne est irréductible et apériodique, elle converge vers elle. C’est une version microscopique du principe d’équilibre des réactions chimiques : à l’équilibre, chaque processus est exactement compensé par son inverse.
En pratique, une simulation se déroule en deux temps.
- La thermalisation. La configuration de départ est arbitraire, donc atypique. On laisse d'abord tourner la chaîne sans rien mesurer, le temps qu'elle « oublie » son point de départ et rejoigne la loi stationnaire. Combien de temps ? C'est exactement la question du temps de mélange, celle des sept battements du jeu de cartes.
- La mesure. On accumule ensuite les valeurs de $A$. Mais les configurations successives sont très corrélées : mille mesures consécutives ne valent pas mille tirages indépendants. Si la chaîne met $\tau$ pas à se décorréler, $n$ mesures n'en valent qu'environ $\frac{n}{2\tau}$, et l'incertitude sur la moyenne décroît comme $\sqrt{2\tau/n}$ au lieu de $1/\sqrt n$.
Tout l’art du MCMC consiste à garder ces deux durées raisonnables. Nous allons voir sur un exemple qu’elles peuvent exploser.
Le modèle d’Ising et sa transition de phase
Le terrain de jeu favori de cet algorithme est le modèle d’Ising, le plus simple des aimants. Sur une grille, chaque site porte un spin qui pointe vers le haut ($+1$) ou vers le bas ($-1$). Deux voisins alignés abaissent l’énergie :
$\displaystyle E=-J\sum_{\text{voisins } i,j}\sigma_i\,\sigma_j .$
Deux tendances s’affrontent. L’énergie pousse les spins à s’aligner tous dans le même sens (un aimant) ; l’agitation thermique les pousse au désordre. Qui l’emporte ?
Pour appliquer Metropolis, la modification proposée est la plus simple possible : retourner un seul spin, choisi au hasard. Tout le charme du modèle est que la variation d’énergie ne dépend que de ses quatre voisins. Les seuls termes de $E$ où figure $\sigma_i$ sont $-J\sigma_i\times(\text{somme des quatre voisins})$, et ils changent de signe quand $\sigma_i$ se retourne :
$\displaystyle \Delta E=2J\,\sigma_i\,(\sigma_{\text{haut}}+\sigma_{\text{bas}}+\sigma_{\text{gauche}}+\sigma_{\text{droite}}).$
Il n’y a que cinq cas possibles, selon le nombre de voisins avec lesquels le spin est d’accord (on prend $J=k=1$) :
| voisins alignés avec le spin | $\Delta E$ si on le retourne | probabilité d’acceptation |
|---|---|---|
| $0$ ou $1$ (il est en minorité) | $-8$ ou $-4$ | $1$ |
| $2$ (égalité) | $0$ | $1$ |
| $3$ | $+4$ | $\mathrm{e}^{-4/T}$ |
| $4$ (il est en plein accord) | $+8$ | $\mathrm{e}^{-8/T}$ |
À $T=1{,}5$, ces deux dernières probabilités valent $7\,\%$ et $0{,}5\,\%$ : un spin entouré de ses semblables ne bouge presque jamais. À $T=4$, elles valent $37\,\%$ et $14\,\%$ : tout s’agite.
Le programme qui anime la figure ci-dessous tient alors en une dizaine d’instructions. En voici le cœur (en JavaScript, à peine simplifié) :
// s : grille de N×N cases valant +1 ou -1
// T : température
function balayage() {
// un « balayage » = N² tentatives
for (let k = 0; k < N * N; k++) {
// 1. un site au hasard
const i = hasard(N), j = hasard(N);
// 2. somme des 4 voisins
// (le « % N » referme la grille)
const voisins = s[(i + 1) % N][j]
+ s[(i + N - 1) % N][j]
+ s[i][(j + 1) % N]
+ s[i][(j + N - 1) % N];
// 3. coût du retournement
const dE = 2 * s[i][j] * voisins;
// 4. règle de Metropolis
if (dE <= 0 || Math.random() < Math.exp(-dE / T))
s[i][j] = -s[i][j]; // accepté
// sinon : refusé, on ne change rien
}
}
Quelques commentaires sur ce code.
- Aucune physique n'est « programmée » : ni domaines, ni température critique, ni transition. Il n'y a que la règle locale. Tout ce que l'on voit émerger à l'écran est une conséquence collective de ces quelques lignes.
- La grille est refermée sur elle-même : grâce au reste de la division par $N$, le voisin de droite de la dernière colonne est la première colonne, et de même verticalement. La grille est donc un tore, sans bord. Ce détail aura son importance.
- L'unité de temps est le balayage : $N^2=25\,600$ tentatives, soit en moyenne une par spin. La figure en effectue deux entre deux images, donc une centaine par seconde.
- Le test
Math.random() < pest la façon standard de réaliser un évènement de probabilité $p$. (Le vrai programme calcule à l'avance $\mathrm{e}^{-4/T}$ et $\mathrm{e}^{-8/T}$, les deux seules valeurs utiles, au lieu d'appeler l'exponentielle des millions de fois par seconde.) - La mesure : à chaque image, on calcule l'aimantation par spin $m=\frac{1}{N^2}\sum_i\sigma_i$, et le point rouge affiche $|m|$, légèrement lissé sur les dernières images.
Faites l’expérience : la grille de $160\times160$ spins ci-dessous évolue selon cette règle, en direct.
Chaque pixel est un spin (bleu vers le haut, clair vers le bas). Le curseur règle la température (avec $J=k=1$). Le point rouge est l'aimantation mesurée, comparée à la valeur exacte. Partez du désordre à basse température : des domaines se forment et grossissent lentement (le bouton « Accélérer » fait défiler le temps dix fois plus vite). Placez-vous à $T_c$ : des amas de toutes tailles apparaissent. Au-dessus, le désordre l'emporte.
Il y a une température critique, $T_c=\frac{2J}{k\ln(1+\sqrt2)}\approx2{,}27\,J/k$. En dessous, une aimantation spontanée apparaît : une majorité de spins s’alignent d’eux-mêmes, sans champ extérieur. Au-dessus, elle est nulle. C’est une véritable transition de phase, celle d’un aimant qui perd son aimantation quand on le chauffe au-delà de son point de Curie (voir la page Transitions de phase). Kramers et Wannier avaient localisé $T_c$ dès 1941 par un argument de symétrie, et Lars Onsager a résolu le modèle exactement en 1944, un tour de force ; mais pour presque tous les autres modèles, seule la simulation par MCMC donne accès aux propriétés de la transition.
Quand l’aimantation reste bloquée : domaines et bandes
Faites l’essai : réglez $T=1{,}5$, bien en dessous de $T_c$, et cliquez sur « Tout désordonner ». La courbe exacte annonce $|m|\approx0{,}99$. Or, une fois sur trois environ, le point rouge se fige à une valeur quelconque ($0{,}2$, $0{,}4$…) et n’en bouge plus. Ce n’est pas un bug.
Il faut bien distinguer deux choses. La courbe d’Onsager décrit l’équilibre, c’est-à-dire la loi stationnaire, la limite $t\to\infty$ de la chaîne. Ce que montre l’écran, c’est la dynamique : le chemin que suit la chaîne pour s’y rendre depuis un état de départ très atypique. Et ce chemin passe par deux étapes.
Une trempe à $T=1{,}5$ depuis le désordre complet, sur la grille de $160\times160$. Le système a choisi « haut » dans une région et « bas » dans une autre ; les deux bandes font le tour de la grille et l'aimantation reste bloquée vers $0{,}4$ au lieu de $0{,}99$.
Première étape : la croissance des domaines. Juste après la trempe, chaque spin s’aligne sur ses voisins immédiats. Mais rien ne dit à une région de choisir « haut » plutôt que « bas » : les deux choix sont équivalents, et des régions éloignées choisissent indépendamment. D’où une mosaïque de petits domaines. Ensuite, seules les parois entre domaines évoluent, et leur moteur est la courbure : un spin situé sur une bosse de la paroi a plus de voisins adverses que de voisins amis, il bascule volontiers, et la bosse se résorbe. Les petits domaines, très courbés, fondent au profit des gros. La taille typique des domaines croît comme la racine carrée du temps, $\ell(t)\propto\sqrt t$ (loi de Lifshitz, Allen et Cahn) : pour doubler la taille des domaines, il faut quatre fois plus de balayages. Pendant toute cette phase, les domaines « haut » et « bas » occupent à peu près la même surface, et $|m|$ reste faible.
Seconde étape : la loterie finale. Quand les domaines atteignent la taille de la grille, deux issues sont possibles.
- Un domaine finit par englober l'autre, dont il ne reste qu'une île ; celle-ci, courbée, rétrécit et disparaît. L'aimantation rejoint la courbe exacte.
- Ou bien un domaine parvient à faire le tour de la grille (rappelez-vous : c'est un tore). On obtient deux bandes, séparées par deux parois droites. Or une paroi droite n'a plus de courbure, donc plus de moteur : chacun de ses spins a trois voisins amis contre un adverse, et le retourner coûte $\Delta E=+4$. Rien ne pousse la paroi dans un sens plutôt que dans l'autre. Le système est piégé, et la valeur de $|m|$ ne reflète que la largeur des bandes, fixée par le hasard de la trempe.
Ce piège est bien connu des physiciens. À température nulle, où l’on n’accepte jamais de hausse d’énergie, les bandes sont figées pour toujours, et l’on a montré (Spirin, Krapivsky et Redner, 2001) qu’environ une trempe sur trois finit ainsi sur un tore carré. Cette proportion se calcule même exactement, grâce à un lien inattendu avec la théorie de la percolation. À température non nulle, les bandes ne sont pas éternelles : les parois fluctuent, errent comme des marcheurs ivres, et finissent par se rencontrer et s’annihiler. Mais cette marche au hasard est d’une lenteur désespérante : sur seize essais numériques à $T=1{,}5$, cinq étaient en bandes après $4000$ balayages, et trois l’étaient encore après $20\,000$, soit près de trois minutes de l’animation en vitesse normale. Pour une grille de côté $L$, le temps de survie des bandes croît comme $L^3$, quand la croissance des domaines n’en demande que $L^2$.
Deux expériences pour s’en convaincre.
- À $T=1{,}5$, cliquez sur « Tout aligner » : le point rouge se pose sur la courbe en un instant, et y reste. L'état d'équilibre est bien celui qu'annonce Onsager ; c'est seulement le chemin depuis le désordre qui est semé d'embûches.
- Une fois les bandes installées, réchauffez vers $T\approx2$ : les parois deviennent rugueuses et mobiles, et les bandes finissent par disparaître bien plus vite. (Remarquez au passage que l'aimantation d'équilibre n'est pas $1$ : à $T=2$, elle vaut $0{,}91$, car l'agitation thermique retourne en permanence quelques spins.)
Pour qui utilise le MCMC, la leçon est sévère. Imaginez que vous ne connaissiez pas la solution d’Onsager : une simulation bloquée dans un état en bandes a toutes les apparences d’une simulation qui a convergé. Les mesures sont stables, les fluctuations faibles, et la valeur obtenue est fausse. Rien, dans une chaîne de Markov, ne signale qu’elle n’a pas encore atteint sa loi stationnaire. D’où les précautions d’usage : lancer plusieurs chaînes depuis des états de départ très différents (ici, l’ordre et le désordre) et vérifier qu’elles donnent la même chose ; refroidir progressivement plutôt que brutalement (c’est l’idée du recuit simulé).
Les mathématiciens ont donné à cette difficulté une forme précise. En dessous de $T_c$, le temps de mélange de cette dynamique sur une grille de côté $L$ croît exponentiellement avec $L$ (pour passer d’une majorité de « haut » à une majorité de « bas », il faut créer une paroi qui traverse toute la grille, ce qui coûte une énergie proportionnelle à $L$). Au-dessus de $T_c$, il est au contraire très court, de l’ordre de $\ln L$ balayages. Et exactement à $T_c$, il croît comme une puissance de $L$ : c’est le ralentissement critique. Les amas de toutes tailles que l’on voit à l’écran évoluent d’autant plus lentement qu’ils sont grands, et le temps de décorrélation $\tau$ croît à peu près comme $L^{2{,}2}$ : précisément là où la physique est la plus intéressante, l’algorithme est le plus lent. On retrouve la question des sept battements : la loi stationnaire existe, mais la chaîne peut mettre un temps démesuré à l’atteindre. Le temps de mélange est le talon d’Achille du MCMC, et l’objet même du livre de Levin et Peres.
Le remède, quand il existe, consiste à changer de chaîne tout en gardant la même loi stationnaire. Plutôt que de retourner les spins un par un, les algorithmes d’amas (Swendsen et Wang en 1987, Wolff en 1989) retournent d’un coup des paquets entiers de spins alignés, choisis de façon à respecter l’équilibre détaillé : le ralentissement critique disparaît presque complètement. Même loi $\pi$, autre matrice de transition, et un temps de mélange sans commune mesure : toute la théorie de cette page est là.
L'algorithme est publié en 1953 par Nicholas Metropolis, Arianna et Marshall Rosenbluth, Augusta et Edward Teller, sous le titre Equation of State Calculations by Fast Computing Machines. Il s'agit de calculer les propriétés d'un fluide de disques durs sur MANIAC, l'un des premiers ordinateurs, construit à Los Alamos pour les calculs des armes nucléaires. Selon les témoignages ultérieurs, l'essentiel du travail revient aux Rosenbluth : Arianna Rosenbluth écrivit le programme, Metropolis dirigeait la machine et signa en premier. En 1970, le statisticien W. K. Hastings généralise la méthode à des propositions non symétriques : c'est l'algorithme de Metropolis-Hastings, aujourd'hui au cœur des statistiques bayésiennes, de la cosmologie à la génétique.
Levin, D. A. et Peres, Y. Markov Chains and Mixing Times, 2e édition, American Mathematical Society (2017). 🌐 Petruszewycz, M. « Chaînes de Markov discrètes dans le domaine linguistique : l’article de 1913 de A. A. Markov », Publications des séminaires de mathématiques et informatique de Rennes, fascicule 2 (1981). 🌐 Hayes, B. « First Links in the Markov Chain », American Scientist 101 (2013). 🌐