On Algorithms: Metropolis, Heat Bath, Molecular Dynamics and Hybrid Monte Carlo

For simplicity, we will start with the matrix harmonic oscillator. The action is given explicitly by
$$\begin{eqnarray} S_{\rm Bo}=\int_{0}^{\beta} dt {\rm Tr}\big[\frac{1}{2}\dot{\Phi}^2+\frac{1}{2}B\Phi^2\big]. \end{eqnarray}$$
In statistical mechanics (quantum field theory) we calculate the partition function (Euclidean path integral) given by
$$\begin{eqnarray} Z=\sum_ne^{-\beta E_n}=\int {\cal D}\Phi\exp\big(-\frac{S_{\rm Bo}[\Phi]}{\hbar}\big). \end{eqnarray}$$
This corresponds to the canonical ensemble where the temperature $\beta\equiv 1/\hbar$ is kept fixed. Here we usually employ the Metropolis algorithm, which is the most general of all Monte Carlo algorithms, to calculate this partition function. This is given explicitly by the probability distribution
$$\begin{equation}P(\Phi_n\longrightarrow \Phi_n^{\prime})={\rm min}(1,\exp(-\Delta S_{\rm Bo}))~,~\Delta S_{\rm Bo}=S_{\rm Bo}[\Phi^{\prime}]-S_{\rm Bo}[\Phi].\label{metropolis}\end{equation}$$
In more detail, the Metropolis algorithm consists of the following steps:

• We start from some initial configuration ${\Phi}_0$.


• We propose a new configuration $\tilde{\Phi}$.


• We compute the variation of the action $\Delta S_{\rm Bo}=S_{\rm Bo}[\Phi]-S_{\rm Bo}[\Phi_0]$.


• We accept the new proposal with a probability $P$ given by ($\ref{metropolis}$).


• We repeat starting from step $2$ until we reach thermalization.


• We measure thermalized configurations which are used to calculate expectation values of physical observables.

The meaning of equation ($\ref{metropolis}$) is as follows. After we compute the variation of the action in step $3$ we check its sign. If ${\Delta} S_{\rm Bo}<0$ then we should accept the proposal since it had resulted in a decrease of the action and thus getting us closer to the minimum. Otherwise, if ${\Delta} S_{\rm Bo}>0$ we accept the proposal only with a probability given by the Boltzmann distribution $\exp(-{\Delta} S_{\rm Bo})$. This is the part which is simulating quantum mechanics and it is implemented numerically via the Von Neumann method. In other words, choose a uniform random number $r$ between $0$ and $1$ and compare it to the Boltzmann distribution $P=\exp(-{\Delta} S)$. If $r<P$ then accept the proposal otherwise reject the proposal.

In many circumstances such as non-local theories it is found that the micro-canonical ensemble is much more favorable than the canonical ensemble. In this case it is the energy that is kept fixed and the preferred algorithm in this case is the hybrid Monte Carlo algorithm which synthesizes together the Metropolis algorithm, the molecular dynamics algorithm and the heat bath algorithm. This alternative formulation relies on the partition function
$$\begin{eqnarray}   Z=\int {\cal D}\Phi{\cal D}P \exp\bigg(-\frac{1}{2}\sum_{n=1}^{\Lambda}{\rm Tr}P^2_n-\frac{S_{\rm Bo}[\Phi]}{\hbar}\bigg). \end{eqnarray}$$
The Euclidean action $S_{\rm Bo}$ acts therefore as a potential term while the new fields $P_n$ act as conjugate momenta associated with the generalized coordinates $\Phi_n$, i.e. $\sum_n{\rm Tr}P_n^2/2$ is a kinetic energy. The Hamiltonian is then given by
$$\begin{eqnarray}   H_{\rm Bo}=\frac{1}{2}\sum_{n=1}^{\Lambda}{\rm Tr}P^2_n+\frac{S_{\rm Bo}[\Phi]}{\hbar}. \end{eqnarray}$$
This Hamiltonian defines a canonical evolution in a fictitious time denoted by $\tau$. Indeed, Hamilton equations of motion read explicitly
$$\begin{eqnarray}   \frac{d({P}_n)_{ab}}{d\tau}=-\frac{\partial H_{\rm Bo}}{\partial (\Phi_n)_{ba}}\Leftrightarrow -\frac{d({P}_n)_{ab}}{d\tau}=\frac{\partial S_{\rm Bo}}{\partial (\Phi_n)_{ba}}\equiv(F_n)_{ab}\nonumber\\    \frac{d({\Phi}_n)_{ab}}{d\tau}=\frac{\partial H_{\rm Bo}}{\partial (P_n)_{ba}}\Leftrightarrow  \frac{d({\Phi}_n)_{ab}}{d\tau}= ({P}_n)_{ab}. \end{eqnarray}$$
The force $F_n$ is given explicitly by

These equations are solved using a molecular dynamics algorithm. In particular, the so-called leap-frog algorithm, which preserves the phase space volume and reversibility and only break Hamiltonian conservation, gives us explicitly the equations  
$$\begin{eqnarray}   &&(P_n)_{ab}(\tau+\frac{\delta\tau}{2})=(P_n)_{ab}(\tau)-\frac{\delta\tau}{2}(F_n)_{ab}(\tau).\label{lf1} \end{eqnarray}$$
$$\begin{eqnarray}   &&(\Phi_n)_{ab}(\tau+\delta\tau)=(\Phi_n)_{ab}(\tau)+\delta\tau (P_n)_{ab}(\tau+\frac{\delta\tau}{2}).\label{lf2} \end{eqnarray}$$
$$\begin{eqnarray}   &&(P_n)_{ab}(\tau+\delta\tau)=(P_n)_{ab}(\tau+\frac{\delta\tau}{2})-\frac{\delta\tau}{2}(F_n)_{ab}(\tau+\delta\tau).\label{lf3}. \end{eqnarray}$$
We remark that the force $F_n$ at the instant $\tau$ is needed to advance all the configurations $\Phi_n$ from $\tau$ to $\tau+\delta \tau$. Thus, we need to calculate the force $F_n$ for all $n$ at instant $\tau$, then apply equations ($\ref{lf1}$) and ($\ref{lf2}$) to advance all $\Phi_n$ from $\tau$ to $\tau+\delta \tau$, after which we calculate again the force $F_n$ for all $n$ at instant $\tau+\delta \tau$, before we finally  apply equation ($\ref{lf3}$) to advance all $P_n$ from $\tau+\delta \tau/2$ to $\tau+\delta \tau$.

The leap-frog will need two extra parameters: the step $\delta \tau$ and the number of iterations which we call $L$. The total time of the motion is given by $T\equiv L\delta \tau$. The initial values of $\Phi_n(0)$ and $P_n(0)$ at time $0$ will be specified and then by applying the above leap-frog algorithm we will obtain the final values of $\Phi_n(T)$ and $P_n(T)$ at time $T$.

Two essential remarks can be stated now:

• The molecular dynamics involves in an obvious way a systematic error.


• The molecular dynamics probes only classical physics.

These two problems can be solved at once via the so-called hybrid Monte Carlo algorithm which is the most general of all Monte Carlo algorithms. This algorithm involves in an essential way the Metropolis algorithm. Indeed, the configuration $\Phi_n(T)$ obtained from the molecular dynamics algorithm is the solution which we will propose as a possible update $\Phi_n^{\prime}$ to the Metropolis algorithm ($\ref{metropolis}$). The probability distribution in this case becomes
$$\begin{eqnarray} P(\Phi_n\longrightarrow \Phi_n^{\prime}, P_n\longrightarrow P_n^{\prime})={\rm min}(1,\exp(-\Delta H_{\rm Bo}))~,~\Delta H_{\rm Bo}=H_{\rm Bo}[\Phi^{\prime},P^{\prime}]-H_{\rm Bo}[\Phi,P].\label{hybrid} \end{eqnarray}$$
As it turns out, the  conjugate momenta $P_n$ should be  updated using a heat bath algorithm in order to avoid ergodic problems, i.e. to be able to reach every point in phase space. This means in particular that $P_n$ should be updated directly from a Gaussian distribution. This is indeed possible since the path integral over the conjugate momenta $P_n$ is only Gaussian. In fact, this path integral can be given by a closed-form expression of the form
$$\begin{eqnarray}   Z&=&\int {\cal D}P \exp(-\frac{1}{2}\sum_{n=1}^{\Lambda}{\rm Tr}P^2_n)\nonumber\\   &=&\bigg(\int \prod_{i}dP_{ii}\prod_{i>j}dP_{ij}dP_{ij}^* \exp(-\frac{1}{2}\sum_{i=1}^NP_{ii}^2-\sum_{i>j}P_{ij}P_{ij}^*)\bigg)^{\Lambda}. \end{eqnarray}$$
All these integrals are Gaussian of the form
$$\begin{eqnarray} \frac{a}{\pi}\int dzdz^*\exp(-azz^*)=\int_0^1 dv_1 \int_0^1 dv_2. \end{eqnarray}$$
$$\begin{eqnarray} z=re^{i\theta}~,~v_1=\exp(-ar^2)~,~v_2=\frac{\theta}{2\pi}. \end{eqnarray}$$
These two equations show that $v_1$ and $v_2$ are uniform random numbers between $0$ and $1$ and thus $z$ is a complex random number given by the formula
$$\begin{eqnarray} z=\sqrt{-\frac{1}{a}\ln v_1}(\cos 2\pi v_2+i\sin 2\pi v_2). \end{eqnarray}$$
The components of the conjugate momenta $(P_n)_{ij}$ are given by Gaussian random numbers of this form.



The action $S=b{\rm Tr}P^2_n$ at each lattice point can also be rewritten as an eigenvalue problem. As it turns out, the eigenvalue distribution is given by the Wigner semi-circle law
$$\begin{eqnarray} \rho(x)=\frac{2}{\pi\delta^2}\sqrt{\delta^2-x^2}~,~\delta^2=\frac{2N}{b}. \end{eqnarray}$$
In summary, the {hybrid Monte Carlo algorithm} is an algorithm in which two crucial extra steps are added to the Metropolis algorithm:

• The step number $2$ in the Metropolis algorithm is implemented via the {\bf molecular dynamics algorithm}. In other words, the new configurations ${\Phi}_n^{\prime}$, $P_n^{\prime}$ are given by the solutions  ${\Phi}_n(T)$, $P_n(T)$ of the molecular dynamics problem with   ${\Phi}_n(0)$, $P_n(0)$ as the initial conditions.




• As it turns out the path integral over the conjugate momenta $P_n$  should be sampled using the so-called {\bf heat bath algorithm} in order to avoid the ergodic problem, i.e. to be able to reach every point in phase space.