Monte Carlo Hamiltonien

Supposons que la distribution cible à échantillonner soit ${\displaystyle f(\mathbf {x} )}$ et qu'une chaîne d'échantillons ${\displaystyle \mathbf {X} _{0},\mathbf {X} _{1},\mathbf {X} _{2},\ldots }$ soit requise. Les équations de la mécanique hamiltonienne se lisent

${\displaystyle {\frac {{\text{d}}x_{i}}{{\text{d}}t}}={\frac {\partial H}{\partial p_{i}}}}$

et

${\displaystyle {\dfrac {{\text{d}}p_{i}}{{\text{d}}t}}=-{\dfrac {\partial H}{\partial x_{i}}}}$

où ${\displaystyle x_{i}}$ et ${\displaystyle p_{i}}$ sont les ${\displaystyle i}$ ème composante du vecteur position et impulsion respectivement et où le hamiltonien ${\displaystyle H}$ est de la forme

${\displaystyle H(\mathbf {x} ,\mathbf {p} )=U(\mathbf {x} )+{\dfrac {1}{2}}\mathbf {p} ^{\text{T}}M^{-1}\mathbf {p} }$

où ${\displaystyle U(\mathbf {x} )}$ est l' énergie potentielle et ${\displaystyle M}$ est une matrice symétrique et définie positive. Dans le but d'échantilloner d'une mesure cible ${\displaystyle f}$, l'énergie potentielle est typiquement donnée par

${\displaystyle U(\mathbf {x} )=-\ln f(\mathbf {x} )}$

Supposons qu'après ${\displaystyle n}$ étapes, la chaîne soit dans l'état ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$et introduisons ${\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}}$ . L'algorithme consiste à proposer un nouvel état ${\displaystyle \mathbf {X} _{n+1}^{*}=\mathbf {x} _{n}(L\Delta t)}$, où ${\displaystyle \mathbf {x} _{n}}$ est une approximation numérique de la dynamique hamiltonienne par la méthode saute-mouton avec ${\displaystyle \Delta t}$ comme pas de discrétisation et où ${\displaystyle L}$ est un entier positif décrivant le nombre de pas à simuler à chaque étape. Le nouvel état proposé ${\displaystyle \mathbf {X} _{n+1}^{*}=\mathbf {x} _{n}(L\Delta t)}$ est ensuite accepté ou rejeté selon la règle de l'algorithme de Metropolis – Hastings.

Plus formellement, soit ${\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}}$ et soit ${\displaystyle \mathbf {p} _{n}(0)}$ tiré de la loi gaussienne ${\displaystyle {\text{N}}\left(\mathbf {0} ,M\right)}$ de moyenne ${\displaystyle \mathbf {0} }$ et de matrice de covariance ${\displaystyle M}$. La méthode saute-mouton consiste à faire évoluer les vecteurs position ${\displaystyle \mathbf {x} _{n}}$ et impulsion ${\displaystyle \mathbf {p} _{n}}$ après le temps ${\displaystyle \Delta t}$ de la façon suivante.

${\displaystyle \mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)=\mathbf {p} _{n}(t)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t)}}$

${\displaystyle \mathbf {x} _{n}(t+\Delta t)=\mathbf {x} _{n}(t)+\Delta tM^{-1}\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}$

${\displaystyle \mathbf {p} _{n}(t+\Delta t)=\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t+\Delta t)}}$

Ces équations doivent être appliquées à ${\displaystyle \mathbf {x} _{n}(0)}$ et ${\displaystyle \mathbf {p} _{n}(0)}$ à ${\displaystyle L}$ reprises afin d'obtenir ${\displaystyle \mathbf {x} _{n}(L\Delta t)}$ et ${\displaystyle \mathbf {p} _{n}(L\Delta t)}$ .

Comme la méthode saute-mouton est une méthode numérique et ne résout pas exactement les équations de la mécanique hamiltonienne, une étape Metropolis – Hastings est utilisée en complément. La transition de ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$ à ${\displaystyle \mathbf {X} _{n+1}}$ est donnée par

${\displaystyle \mathbf {X} _{n+1}|\mathbf {X} _{n}=\mathbf {x} _{n}={\begin{cases}\mathbf {x} _{n}(L\Delta t)&{\text{with probability }}\alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)\\\mathbf {x} _{n}(0)&{\text{otherwise}}\end{cases}}}$

où

${\displaystyle \alpha \left(\mathbf {x} _{n}(0)_{n},\mathbf {x} _{n}(L\Delta t)\right)={\text{min}}\left(1,{\dfrac {\exp \left[-H(\mathbf {x} _{n}(L\Delta t),\mathbf {p} _{n}(L\Delta t))\right]}{\exp \left[-H(\mathbf {x} _{n}(0),\mathbf {p} _{n}(0))\right]}}\right)}$

Cette opération est ensuite répétée afin d'obtenir ${\displaystyle \mathbf {X} _{n+1},\mathbf {X} _{n+2},\mathbf {X} _{n+3},\ldots }$ .

• Portail de la physique