Monte Carlo Hamiltonien

En physique computationnelle et en statistiques, l' algorithme de Monte Carlo hamiltonien (également connu sous le nom de Monte Carlo hybride), est une méthode de Monte Carlo par chaîne de Markov dont l'objectif est d'obtenir une séquence d'échantillons aléatoires qui convergent selon une distribution de probabilité cible, typiquement difficile à échantillonner directement. Cette séquence peut notamment être utilisée pour estimer des intégrales par rapport à une distribution cible (calcul d'espérances mathématiques).

Le Monte Carlo hamiltonien correspond à une instance de l'algorithme de Metropolis – Hastings où les déplacements proposés dans l'espace d'états sont issus d'un processus gouverné par une dynamique hamiltonienne[1],[2]et simulée à l'aide d'un intégrateur numérique réversible et préservant le volume (généralement la méthode saute-mouton).

Comparé à l'utilisation d'une distribution de proposition de marche aléatoire gaussienne dans l'algorithme Metropolis – Hastings, la méthode de Monte Carlo Hamiltonien réduit la corrélation entre les états échantillonnés successivement en proposant des déplacements vers des états distants qui maintiennent une forte probabilité d'acceptation en raison des propriétés de conservation d'énergie approximatives de la dynamique hamiltonienne simulée avec de l'utilisation d'un intégrateur symplectique . La corrélation réduite signifie que moins d'échantillons de chaîne de Markov sont nécessaires pour approcher les intégrales par rapport à la distribution de probabilité cible pour une erreur de Monte Carlo donnée. L'algorithme a été proposé à l'origine par Simon Duane, Anthony Kennedy, Brian Pendleton et Duncan Roweth en 1987 [3] pour des calculs en chromodynamique quantique sur un réseau .

Algorithme

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

d x i d t = ∂ H ∂ p i {\displaystyle {\frac {{\text{d}}x_{i}}{{\text{d}}t}}={\frac {\partial H}{\partial p_{i}}}}

et

d p i d t = − ∂ H ∂ x i {\displaystyle {\dfrac {{\text{d}}p_{i}}{{\text{d}}t}}=-{\dfrac {\partial H}{\partial x_{i}}}}

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

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

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

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


Supposons qu'après n {\displaystyle n} étapes, la chaîne soit dans l'état X n = x n {\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}} et introduisons x n ( 0 ) = x n {\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}} . L'algorithme consiste à proposer un nouvel état X n + 1 ∗ = x n ( L Δ t ) {\displaystyle \mathbf {X} _{n+1}^{*}=\mathbf {x} _{n}(L\Delta t)} , où x n {\displaystyle \mathbf {x} _{n}} est une approximation numérique de la dynamique hamiltonienne par la méthode saute-mouton avec Δ t {\displaystyle \Delta t} comme pas de discrétisation et où L {\displaystyle L} est un entier positif décrivant le nombre de pas à simuler à chaque étape. Le nouvel état proposé X n + 1 ∗ = x n ( L Δ t ) {\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 x n ( 0 ) = x n {\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}} et soit p n ( 0 ) {\displaystyle \mathbf {p} _{n}(0)} tiré de la loi gaussienne N ( 0 , M ) {\displaystyle {\text{N}}\left(\mathbf {0} ,M\right)} de moyenne 0 {\displaystyle \mathbf {0} } et de matrice de covariance M {\displaystyle M} . La méthode saute-mouton consiste à faire évoluer les vecteurs position x n {\displaystyle \mathbf {x} _{n}} et impulsion p n {\displaystyle \mathbf {p} _{n}} après le temps Δ t {\displaystyle \Delta t} de la façon suivante.

p n ( t + Δ t 2 ) = p n ( t ) − Δ t 2 ∇ U ( x ) | x = x n ( t ) {\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)}}
x n ( t + Δ t ) = x n ( t ) + Δ t M − 1 p n ( t + Δ t 2 ) {\displaystyle \mathbf {x} _{n}(t+\Delta t)=\mathbf {x} _{n}(t)+\Delta tM^{-1}\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}
p n ( t + Δ t ) = p n ( t + Δ t 2 ) − Δ t 2 ∇ U ( x ) | x = x n ( t + Δ t ) {\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 à x n ( 0 ) {\displaystyle \mathbf {x} _{n}(0)} et p n ( 0 ) {\displaystyle \mathbf {p} _{n}(0)} à L {\displaystyle L} reprises afin d'obtenir x n ( L Δ t ) {\displaystyle \mathbf {x} _{n}(L\Delta t)} et p n ( L Δ t ) {\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 X n = x n {\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}} à X n + 1 {\displaystyle \mathbf {X} _{n+1}} est donnée par

X n + 1 | X n = x n = { x n ( L Δ t ) with probability  α ( x n ( 0 ) , x n ( L Δ t ) ) x n ( 0 ) otherwise {\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ù

α ( x n ( 0 ) n , x n ( L Δ t ) ) = min ( 1 , exp ⁡ [ − H ( x n ( L Δ t ) , p n ( L Δ t ) ) ] exp ⁡ [ − H ( x n ( 0 ) , p n ( 0 ) ) ] ) {\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 X n + 1 , X n + 2 , X n + 3 , … {\displaystyle \mathbf {X} _{n+1},\mathbf {X} _{n+2},\mathbf {X} _{n+3},\ldots } .

Voir aussi

Articles connexes

Bibliographie

  • (en) Betancourt, « A Conceptual Introduction to Hamiltonian Monte Carlo », arXiv,‎ (Bibcode 2017arXiv170102434B, arXiv 1701.02434)
  • Liu, Jun S. (2004). Stratégies de Monte Carlo en informatique scientifique . Série Springer dans les statistiques, Springer. 189-203. (ISBN 978-0-387-76369-9).

Liens externes

  • Betancourt, « Efficient Bayesian inference with Hamiltonian Monte Carlo », MLSS Iceland 2014, sur YouTube
  • Hamiltonian Monte Carlo à partir de zéro
  • Optimisation et méthodes de Monte Carlo

Notes et références

  1. ↑ Buchholz, Alexander, « High dimensional Bayesian computation, Section 1.2 échantillonnage Monte Carlo », these.fr,‎
  2. ↑ (en) Radford Neal, « MCMC Using Hamiltonian Dynamics », dans Handbook of Markov Chain Monte Carlo, vol. 20116022, Chapman and Hall/CRC, (ISBN 978-1-4200-7941-8, DOI 10.1201/b10905-6, lire en ligne)
  3. ↑ Duane, Kennedy, Anthony D., Pendleton, Brian J. et Roweth, Duncan, « Hybrid Monte Carlo », Physics Letters B, vol. 195, no 2,‎ , p. 216–222 (DOI 10.1016/0370-2693(87)91197-X, Bibcode 1987PhLB..195..216D)
  • icône décorative Portail de la physique