Instabilities and Non-Reversibility of Molecular Dynamics Trajectories

TL;DR

This paper investigates the causes of instability and irreversibility in molecular dynamics trajectories used in Hybrid Monte Carlo algorithms, analyzing numerical and chaotic factors, and providing evidence that these issues are manageable in quantum field theory simulations.

Contribution

It identifies two main sources of trajectory instability—numerical inaccuracies and intrinsic chaos—and analyzes their effects, offering insights into their impact on Hybrid Monte Carlo methods.

Findings

Instabilities due to numerical integration are finite volume effects.

Chaotic behavior's Lyapunov exponent decreases towards the continuum limit.

Empirical data from lattice QCD supports theoretical analysis.

Abstract

The theoretical justification of the Hybrid Monte Carlo algorithm depends upon the molecular dynamics trajectories within it being exactly reversible. If computations were carried out with exact arithmetic then it would be easy to ensure such reversibility, but the use of approximate floating point arithmetic inevitably introduces violations of reversibility. In the absence of evidence to the contrary, we are usually prepared to accept that such rounding errors can be made small enough to be innocuous, but in certain circumstances they are exponentially amplified and lead to blatantly erroneous results. We show that there are two types of instability of the molecular dynamics trajectories which lead to this behavior, instabilities due to insufficiently accurate numerical integration of Hamilton's equations, and intrinsic chaos in the underlying continuous fictitious time equations of motion themselves. We analyze the former for free field theory, and show that it is essentially a finite volume effect. For the latter we propose a hypothesis as to how the Liapunov exponent describing the chaotic behavior of the fictitious time equations of motion for an asymptotically free quantum field theory behaves as the system is taken to its continuum limit, and explain why this means that instabilities in molecular dynamics trajectories are not a significant problem for Hybrid Monte Carlo computations. We present data for pure $SU(3)$ gauge theory and for QCD with dynamical fermions on small lattices to illustrate and confirm some of our results.

Equations

$$H=\frac{1}{2}\sum _{p\in {\text{bbxii}Z}_V}\left({\pi }_p^2+{\omega }_p^2{\varphi }_p^2\right).$$

$$U(\delta\tau)=\left(\begin{array}[]{cc}1-{\scriptstyle{1\over 2}}\delta\tau^{2}&\delta\tau\\ -\delta\tau+{\scriptstyle{1\over 4}}\delta\tau^{3}&1-{\scriptstyle{1\over 2}}\delta\tau^{2}\end{array}\right).$$

$$U(\delta\tau)=\left(\begin{array}[]{cc}\cos[\kappa(\delta\tau)\delta\tau]&{\displaystyle\sin[\kappa(\delta\tau)\delta\tau]\over\displaystyle\rho(\delta\tau)}\\ -\rho(\delta\tau)\sin[\kappa(\delta\tau)\delta\tau]&{\cos[\kappa(\delta\tau)\delta\tau]}\end{array}\right),$$

$$\kappa (\delta \tau )=\frac{{\cos }^{-1}(1-\frac{1}{2}\delta {\tau }^2)}{\delta \tau },\rho (\delta \tau )=\sqrt{1-\frac{1}{4}\delta {\tau }^2}.$$

$$U(\tau)=\left(\begin{array}[]{cc}\cos[\kappa(\delta\tau)\tau]&{\displaystyle\sin[\kappa(\delta\tau)\tau]\over\displaystyle\rho(\delta\tau)}\\ -\rho(\delta\tau)\sin[\kappa(\delta\tau)\tau]&{\cos[\kappa(\delta\tau)\tau]}\end{array}\right),$$

$${\delta H}={1\over 2}\sum_{p\in{\bbxii Z}_{V}}\left(\begin{array}[]{c}\phi_{p}\\ \pi_{p}\end{array}\right)^{T}\Bigl{[}U_{p}^{T}U_{p}-{\bbxii I}\Bigr{]}\left(\begin{array}[]{c}\phi_{p}\\ \pi_{p}\end{array}\right).$$

$$P_{\delta H}(\xi )$$

$$\nu =\pm \frac{1}{\delta \tau }\mathrm{Re}\ln \left[\frac{1}{2}\delta {\tau }^2-1\pm \delta \tau \sqrt{\frac{1}{4}\delta {\tau }^2-1}\right],$$

$$P_{\delta H}(\xi )={\int }_{-\infty }^{\infty }\frac{d\eta }{2\pi }\frac{e^{i\eta \xi }}{\sqrt{1+T^{2}i\eta (1-i\eta )}}=e^{\xi /2}{\int }_{i/2}^{\infty +i/2}\frac{d{\eta }^{\prime }}{\pi }\frac{\cos ({\eta }^{\prime }\xi )}{\sqrt{1+T^2(\frac{1}{4}+{\eta }^{\prime 2})}}.$$

$$P_{\delta H}(\xi )=\frac{e^{\xi /2}}{\pi T}{\int }_0^{\infty }\frac{dx}{\sqrt{1+x^2}}\cos \left(\frac{x\xi }{2}\sqrt{1+\frac{4}{T^2}}\right),$$

$$P_{\delta H}(\xi )=\frac{e^{\xi /2}}{\pi T}{K}_0\left(\frac{|\xi |}{2}\sqrt{1+\frac{4}{T^2}}\right).$$

$$P_{\hbox{\rm\tiny acc}}\equiv\Big{\langle}\min(1,e^{-{\delta H}})\Big{\rangle}=\int_{-\infty}^{\infty}d\xi\,\min(1,e^{-\xi})P_{\delta H}(\xi).$$

$$P_{\mathrm{acc}}=\frac{2}{\pi T}{\int }_{-\infty }^{\infty }\frac{dx}{\sqrt{1+x^2}\left[1+x^2\left(1+\frac{4}{T^2}\right)\right]};$$

$$P_{\mathrm{acc}}=\frac{\sin \theta }{\pi }{\int }_0^{\infty }\frac{dy}{y^2+2y\cos \theta +1}=\frac{\theta }{\pi }=\frac{2}{\pi }{\tan }^{-1}\left(\frac{2}{T}\right),$$

$$P_{\mathrm{acc}}$$

$${1\over\prod_{q=1}^{N}\sqrt{1+T_{q}^{2}({\scriptstyle{1\over 4}}+\eta^{\prime 2})}}={({\scriptstyle{1\over 4}}+\eta^{\prime 2})^{-N/2}\over\prod_{q}T_{q}}\biggl{(}1-{\scriptstyle{1\over 2}}\sum_{q}{1\over T_{q}^{2}}({\scriptstyle{1\over 4}}+\eta^{\prime 2})^{-1}+\cdots\biggr{)}.$$

$$\int_{0}^{\infty}d\eta^{\prime}\,e^{-\eta^{\prime 2}/2\sigma^{2}}({\scriptstyle{1\over 4}}+\eta^{\prime 2})^{-k/2}=2^{k-2}\sqrt{\pi}U\Bigl{(}{1\over 2},{3-k\over 2};{1\over 8\sigma^{2}}\Bigr{)}$$

$$\nu \propto e^{-\beta /12{\beta }_0}\mathrm{with}{\beta }_0=\frac{11-\frac{2}{3}{n}_f}{16{\pi }^2},$$

Full text

1 Introduction

The goal of this paper is to study the effect of errors introduced by working with finite-precision arithmetic on Hybrid Monte Carlo (HMC) and related algorithms.

The HMC algorithm [1, 2, 3] allows us to generate an ensemble of field configurations which are selected from some probability distribution. The algorithm is “exact,” in the sense that it is a Markov process which converges to the desired distribution with no systematic errors provided that the fictitious time molecular dynamics (MD) trajectories within it are exactly reversible and area preserving, that the computation of the action for the Metropolis accept/reject step is exact, and that we have a supply of perfectly random numbers. It does not require that the MD integration or any conjugate gradient (CG) “inversion” of the fermionic kernel be carried out exactly111The initial CG vector must be chosen in a time-symmetric way.. Indeed, if we were to carry out all computations using exact arithmetic (i.e., with no rounding errors whatsoever) then the leapfrog integration scheme provides a method of integrating Hamilton’s equations while maintaining exact reversibility and area preservation.

All numerical computations carried out using floating point arithmetic are subject to rounding errors, but unless these errors are amplified exponentially we do not normally consider them to be a serious problem. The reason for this is that if a sensible rounding mode222If we were to truncate all floating point numbers towards zero after each operation then the error after $N$ operations would be $\varepsilon N$, and the cost would grow twice as rapidly with $N$ — but still logarithmically. is used for the elementary floating point operations then the error after a sequence of $N$ operations, each of which gives rise to an error of magnitude $\varepsilon$, is $\varepsilon\sqrt{N}$. The cost of working with $d$ fractional digits grows approximately linearly333On a typical computer chip doubling the number of digits requires doubling the number of gates in the adder and multiplier and increasing the clock period to allow for the worst-case carry propagation. Thus asymptotically the number of gates, the clock period, and the amount of memory would increase linearly with the number of digits $d$. In practice it is better to reduce the carry propagation time to grow as $\ln d$ at the cost of increasing the number of gates as $d\ln d$. This leads to an asymptotic cost which grows as $d(\ln d)^{2}$. with $d$ and gives $\varepsilon=10^{-d}$. Therefore if we require an answer with an error smaller in magnitude than $\delta$ in magnitude the cost is proportional to $-\log_{10}(\delta/\sqrt{N})$, which grows only like $\ln N$. This is a very small correction444For real computers the number of digits available is often “quantized” in units of words, and increasing the precision from about $7$ digits (IEEE single precision) to about $14$ digits (IEEE double precision) may be large, so asymptotic estimates must be treated with caution. to the growth of the cost of the HMC algorithm as the volume and correlation length of the system is increased.

All the arguments we present as to the conditions under which reversibility is maintained apply to area preservation as well, but we shall concentrate on the former as the latter is less easy to determine empirically.

2 A Brief Survey of Algorithms

Let us very briefly summarize the various advantages and disadvantages of several variants of the HMC algorithm for dynamical fermion Quantum Chromodynamics (QCD) computations.

2.1 Hybrid Monte Carlo (HMC)

The algorithm is exact even if only a cheap approximate CG solution is carried out for computing the force along a MD trajectory, and such a CG solution is needed once per integration step. Nevertheless, an accurate CG solution is needed for the Metropolis accept/reject test once per trajectory.

If we believe that the behavior of free field theory carries over to interacting theories like QCD, then choosing trajectory lengths which grow linearly with the correlation length of the system $\tau\propto\xi$ minimizes the cost per independent configuration and give a dynamical critical exponent $z\approx 1$. The coefficient of proportionality depends upon which operator is being optimized, but it will typically be of order unity.

2.2 Second Order Langevin Monte Carlo (L2MC)

This algorithm was introduced by Horowitz [4] (and is sometimes also called Kramers algorithm). It is similar to HMC except that only single step (Langevin) trajectories are taken, and that instead of choosing the fictitious momenta from a Gaussian heatbath after each step the final momenta of a trajectory are reversed and mixed with a small amount of Gaussian noise to serve as the initial momenta for the next trajectory.

The advantages are that because the trajectories are short they are not subject to any large exponential amplification of rounding errors (which will be discussed in this paper), and that the the free field theory mean acceptance rate per trajectory is $\mathop{\rm erfc}(c\sqrt{V\delta\tau^{6}})$ whereas for HMC it is $\mathop{\rm erfc}(c^{\prime}\sqrt{V\delta\tau^{4}})$, where $V$ is the number of lattice sites, $\delta\tau$ is the leapfrog integration step size, and $c$ and $c^{\prime}$ are constants of order unity [5, 6, 7, 8].

The disadvantages are that an accurate CG solution is required for every integration step (because each such step has an accept/reject test following it), and that upon a rejection the trajectory reverses itself (because the momenta are flipped). This last observation means that in order to minimize autocorrelations the acceptance rate has to be very close to one with a concomitantly tiny $\delta\tau$.

A free field theory analysis [7, 8] indicates that L2MC and HMC have essentially the same scaling behavior in both the thermodynamic and continuum limits, so the choice between them must be made on the grounds of the cost of specific implementations and on the possible costs associated with instabilities occurring for long trajectories in HMC.

2.3 Generalized Hybrid Monte Carlo (GHMC)

The Generalized Hybrid Monte Carlo (GHMC) algorithm [9, 8] combines the HMC and L2MC methods giving an algorithm with three tunable parameters: the trajectory length $\tau$, the mixing angle between the old fictitious momenta and Gaussian noise $\theta$, and the integration step size $\delta\tau$. A free field theory analysis indicates that if $\delta\tau$ is chosen to given a good acceptance rate then there is a valley in the $(\tau,\theta)$ plane for which the dynamical critical exponent $z=1$ and furthermore for which the cost per independent configuration does not vary much.

As HMC and L2MC are both special cases of GHMC, the optimal choice of parameters may well not correspond to either case but in fact be somewhere in the middle of this valley.

3 Causes of Irreversibility

Empirical studies [10, 11] have shown that HMC computations for gauge theories are not exactly reversible. Given that they are performed using floating point arithmetic this is not at all surprising. However, if $\delta\tau$ or, for dynamical fermion calculations, the CG residual are taken to be too large then very large violations of reversibility are observed. This large violation of reversibility was observed even when the initial CG vector was chosen to be zero.

We recall that the HMC algorithm is exactly reversible unless the initial CG vector is chosen in a time asymmetric way555If the CG solution was found exactly then it would of course be independent of the initial vector. In practice only enough CG iterations are carried out to reduce the residual to below some preset value, and the Krylov space thus explored depends on the initial vector. If this initial vector depends on the solution for the previous step then these Krylov spaces will be different for forward and backward trajectories, and we will induce reversibility errors depending on the CG residual used. except for the effects of finite-precision arithmetic. We are thus presented with the question of why the rounding errors are amplified by a large factor.

Three possible causes of this amplification suggest themselves:

• •

A small error changes the Metropolis accept/reject choice;

• •

A small error changes the number of CG iterations;

• •

The MD trajectories exhibit an exponential sensitivity to initial conditions (in other words they are unstable or chaotic) [11, 12].

It is easy to see that although the first two mechanisms produce a large change they occur correspondingly infrequently, and therefore are not compatible with the nature of the irreversibility observed. We therefore conclude that the irreversibility observed is caused by instabilities in the MD trajectories.

One may wonder to what extend irreversibility introduces systematic errors into the results of an HMC calculation. If the irreversibility is due to some underlying chaos in the equations of motion perhaps this just has the same effect as a different choice of random numbers for the fictitious momentum heatbath? Indeed, it is hard to see significant effects on physical observables unless truly huge violations of reversibility are induced (in which case the trajectories often become so unstable as to cause numerical overflows). We did observe666See for example Figure 8 on page 8. that when the trajectories become irreversible there is often a large change in the fictitious kinetic energy as well (even though the total fictitious energy ${\delta H}$ stays small). This means that the distribution of fictitious momenta is clearly wrong, and thus the computation is erroneous even if the effect on physically interesting observables is hard to detect.

4 Causes of Instability

There are two distinct mechanisms which cause the MD equations to become unstable:

• •

The discrete integration procedure (leapfrog) diverges exponentially from solutions of the classical equations of motion. This instability should grow with the number of integration steps, and thus have a short characteristic time scale for large volumes.

• •

The classical equations of motion are themselves chaotic777Chaos in classical Yang-Mills dynamics (but not in fictitious time) has been investigated by several authors. See [13] and references cited therein.. These may be expected to have a time scale independent of the step size, and presumably determined by some characteristic length scale of the system.

In section 5 we analyze the instabilities of the leapfrog integration scheme for free field theory, and in section 6 we consider the chaotic nature of the underlying continuous time equations of motion themselves and propose a hypothesis as to how the time scale of the corresponding instabilities behaves as the system approaches the continuum limit.

5 Free Field Theory Analysis

The phenomenon of the leapfrog integration scheme diverging from the true trajectory occurs even for free field theory, and it is instructive to see what happens there [7, 8].

Consider a system of harmonic oscillators $\{\phi_{p}\}$ for $p\in{\bbxii Z}_{V}$. The Hamiltonian on fictitious phase space is

[TABLE]

5.1 Single Mode Analysis

The Hamiltonian is diagonal, so we may temporarily consider the evolution of a single mode and set its frequency $\omega_{p}=1$. In fact, we will find that the essential features of the instability are displayed by a single mode.

The leapfrog discretization of Hamilton’s equations can be written as the following matrix acting on the phase space vector $(\phi,\pi)$

[TABLE]

The general area-preserving reversible linear map on $(\phi,\pi)$ phase space888If F\equiv\left(\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right) is the momentum-flip operation then $F^{-1}U(\delta\tau)F=U(-\delta\tau)$, from which it immediately follows that $U_{1,1}$ and $U_{2,2}$ are even functions of $\delta\tau$, and $U_{1,2}$ and $U_{2,1}$ are odd functions of $\delta\tau$. Furthermore $U(\delta\tau)U(-\delta\tau)={\bbxii I}$ implies that $U_{1,1}=U_{2,2}$ except for the trivial case where $U_{1,2}=U_{2,1}=0$. may be parameterized as

[TABLE]

where $\kappa$ and $\rho$ are even functions, but not necessarily real valued. For the lowest-order leapfrog algorithm $\kappa$ and $\rho$ are

[TABLE]

With this parameterization it is easy to see that the evolution operator for a trajectory of $\tau/\delta\tau$ leapfrog steps is

[TABLE]

and thus the change in fictitious energy over this trajectory is999The subscript $p$ indicates that all fictitious times are to be multiplied by the frequency $\omega_{p}$.

[TABLE]

It is clear that some kind of critical behavior will occur at $\delta\tau=2$.

Given this evolution matrix we can compute the probability distribution of ${\delta H}$,

[TABLE]

where $T\equiv\left|\left(\rho-{\displaystyle 1\over\displaystyle\rho}\right)\sin(\kappa\tau)\right|$. It is easy to see that $T$ rapidly approaches its asymptotic behavior of growing proportionally to $e^{\nu\tau}$ where

[TABLE]

(which falls as $2\ln(\delta\tau)/\delta\tau$ as $\delta\tau\to\infty$). For $\delta\tau\leq 2$ not only is the exponent $\nu$ zero, but $T$ is bounded by a constant. That the exponent decreases as $\delta\tau\to\infty$ merely reflects the fact that $T$ is growing exponentially in the number of MD steps but only algebraically in the step size $\delta\tau$. This is illustrated in Figure 1.

For a single mode the probability distribution of ${\delta H}$ may be evaluated by changing variable to $\eta^{\prime}\equiv\eta+i/2$

[TABLE]

The singularities in the integrand occur when $1+T^{2}({\scriptstyle{1\over 4}}+\eta^{\prime 2})=0$, for which $\eta^{\prime 2}<-{\scriptstyle{1\over 4}}$ and hence $|\mathop{\rm Im}\eta^{\prime}|>{\scriptstyle{1\over 2}}$. We may thus shift the integration contour to the real axis, and defining $x\equiv T\eta^{\prime}/\sqrt{1+{\scriptstyle{1\over 4}}T^{2}}$ we obtain

[TABLE]

which can be expressed in closed form in terms of a modified Bessel function

[TABLE]

The logarithm of the acceptance rate is shown in Figure 2.

The average Metropolis acceptance rate is

[TABLE]

For a single mode we may use eq. (3), and evaluating the integral over $\xi$ gives

[TABLE]

setting $x\equiv{1\over 2}\left(\sqrt{y}-{\displaystyle 1\over\displaystyle\sqrt{y}}\right)$ and ${2\over T}\equiv\tan{\theta\over 2}$ we easily obtain

[TABLE]

and this is shown in Figure 3 as a function of the step size and trajectory length. Observe that there is a “wall” at $\delta\tau=2$ which separates the region where the acceptance rate oscillates as a function of $\tau$ from the region where it plummets exponentially.

5.2 Multiple Mode Analysis

These results can be generalized to a field theory with many “stable” modes and some number of “unstable” ones. In the general case let us consider the acceptance rate $P_{\hbox{\rm\tiny acc}}$ defined in eq. (4) and use the expression of eq. (1) for the distribution of values of ${\delta H}$. As before we may make the change of variable $\eta^{\prime}\equiv\eta+i/2$ and observe that there are no singularities in the strip of the complex $\eta^{\prime}$ plane satisfying $|\mathop{\rm Im}\eta^{\prime}|<{\scriptstyle{1\over 2}}$. We thus obtain

[TABLE]

Let us now assume that there are a small number $N$ of “unstable” modes with $T_{q}\gg 1$ and a large number $V$ of “stable” modes with $T_{p}\ll 1$ such that101010If $\delta\tau\ll 1$ then $T_{p}=-{\scriptstyle{1\over 4}}(\omega_{p}\delta\tau)^{2}\sin(\omega_{p}\tau)+O(\delta\tau^{4})$, and thus $\sum_{p=1}^{V}T_{p}^{2}=kV\delta\tau^{4}+O(\delta\tau^{6})$ where the spectral average $k\equiv{1\over 16V}\sum_{p=1}^{V}\omega_{p}^{4}\sin^{2}(\omega_{p}\tau)$ is of order unity. In order to obtain a non-negligible acceptance rate we must take $\delta\tau$ to be of order $V^{-{\scriptstyle{1\over 4}}}$, in which case the quantity $\sigma$ will also be of order one. $\sum_{p=1}^{V}T_{p}^{2}\equiv\sigma^{-2}$ and $\sum_{p=1}^{V}T_{p}^{2n}=O(V^{1-n})$. We can thus compute an asymptotic expansion valid for large $V$,

[TABLE]

where the “stable” modes correspond to subscript $p$ and the unstable ones to subscript $q$. If there are a few modes with $V^{-1/4}\ll\omega_{p}\delta\tau\ll 1$ then these may be taken into account by computing further terms in the asymptotic expansion above. On the other hand, we may carry out an expansion in powers of $1/T_{q}$ for the “unstable” modes by writing

[TABLE]

Using this expansion in the integrand above reduces both the large $T$ and small $T$ expansions to a sum of integrals of the form

[TABLE]

where $U$ is a confluent hypergeometric function.111111This is related to Whittaker’s function $W_{\lambda,\mu}(z)=z^{\mu+{\scriptstyle{1\over 2}}}e^{-z/2}U(\mu-\lambda+{\scriptstyle{1\over 2}},2\mu+1;z)$. It is easily verified that for $N=0$ the leading term in the asymptotic expansion is just $\mathop{\rm erfc}(1/\sqrt{8\sigma^{2}})$.

In Figure 4 we show the $P_{\hbox{\rm\tiny acc}}$ for the case where there are many “stable” modes with $\sigma=1$ as a function of the value of $\sigma$ for one additional “unstable” mode. It is clear that the expansions calculated above are very good as long as $T$ is not close to one.

In Figure 5 we plot $\ln P_{\hbox{\rm\tiny acc}}$ as a function of $T$ for various values of $N$, the number of “unstable” modes. Each additional unstable mode clearly reduces the acceptance rate by a large factor for these parameters, and this makes it reasonable to assume that at least the qualitative behavior of the instabilities in the system is apparent from the case where there is just a single unstable mode.

In Figure 6 we show a comparison of the dependence of the mean acceptance rate on $T$ for a single mode and for a system of many “stable” modes and one “unstable” one. This demonstrates that the behavior of a single mode dominates the instabilities in the case of free field theory.

We notice that in order to keep the acceptance rate constant as the lattice volume $V\to\infty$ we must decrease $\delta\tau$ so as to keep $V\delta\tau^{4}$ fixed. Thus as we approach the thermodynamic limit the instabilities will go away: in other words the leapfrog instability is a finite volume effect.

5.3 Interacting Field Theories

When interactions are present it is no longer meaningful to talk of independent modes, but for weak interactions like those for asymptotically free field theories at short distance it is still useful to consider the system as a set of weakly coupled modes. We expect that the onset of the integration instability will still be caused by the highest frequency mode, and will occur at $\omega_{\max}\delta\tau=2$. The forces acting on the highest frequency mode due to the other modes will fluctuate in some complicated way, and consequently we would expect that the “wall” at $\omega_{\max}\delta\tau=2$ gets smeared out.

In order to illustrate this effect we have made numerical studies of a simple model of a single harmonic oscillator whose frequency is randomly chosen from a Gaussian distribution with mean unity and standard deviation $\sigma$ before each MD step. The behavior of this model is just like that of the free field theory considered above, except that the “wall” at $\delta\tau=2$ gets spread out. Our numerical results are shown in Figure 7, where the analytic result for $\sigma=0$ is taken from eq. (2).

We conjecture that interacting asymptotically free field theories like non-Abelian gauge theories behave similarly, and this is supported by our numerical results.

We do not expect the introduction of dynamical fermions to produce any qualitatively different behavior. For QCD computations using the HMC algorithm the pseudofermions produce a force which we expect to be proportional to $1/m_{f}$, so the “highest frequency” of the system, which is responsible for the integration instability, will grow as $\omega_{\max}\propto 1/m_{f}$ too. In particular we expect $\omega_{\max}$ to become large as $m_{f}\to 0$, so the critical value of $\delta\tau$ will become small as $\kappa\to\kappa_{c}$ for Wilson fermions. An inaccurate CG solution will produce errors which are multiplied by the scale of this fermionic force, and will become important when the size of the errors in the CG solution become comparable with the size of the fluctuations in the forces due to the fields themselves.

It is interesting to observe that reversibility can be violated significantly even though the acceptance rate does not become small. For example we will see an example of this in Figure 10 on page 10 for $\delta\tau=0.005$ and $\tau=50$. This is because the leapfrog integration scheme tries to conserve energy even when it has wandered far away from the true continuous time trajectory. It is therefore prudent to verify that reversibility is satisfied to the precision required, even when the HMC algorithm has an average acceptance rate which is close to one.

6 Chaos in Continuous Time Evolution

Jansen and Liu [11] observed that there is a second cause of instability in the MD trajectories, which occurs even when the integration step size $\delta\tau\to 0$. This is because the underlying continuous fictitious time equations of motion are themselves chaotic. In the chaotic regime two nearby classical trajectories diverge exponentially with a characteristic exponent $\nu$ called the Liapunov exponent. Obviously such behavior cannot be studied in the context of free field theory.

The instability due to chaos cannot be removed by reducing $\delta\tau$, and is therefore not a finite size effect. If the system we are studying exhibits critical slowing down,121212This is probably not the case for most current lattice QCD computations. and we believe that the mechanism for reducing this to $z=1$ in free field theory is applicable to this system too, and we wish to use the HMC algorithm, then we must scale the trajectory length with the correlation length, $\tau\propto\xi$. This means that there is the potential for exponential amplification of rounding errors.

Even if there is such an exponential amplification of rounding errors, reversibility can be maintained to any desired precision by a linear increase in the number of digits used for floating point arithmetic. This is no longer a negligible addition to the computational complexity of the problem (i.e., it is not just a logarithmic correction), but it is still a small correction to the power dependence of the computational cost on the correlation length and volume.

This problem can be avoided using the GHMC algorithm, but not only is this unnecessary until chaotic amplification of rounding errors becomes important, but also as we shall see it might never be necessary at all.

The reason for this is that the Liapunov exponent (averaged over the equilibrium distribution) is not constant as a function of $\beta$. If the chaotic dynamics is not only a property of the underlying continuous fictitious time evolution, but is also a property of the underlying (space-time) continuum field theory, then the Liapunov exponent would be constant when measured in “physical” units, that is $\nu\xi$ would be constant as $\xi\to\infty$. In other words this hypothesis would say that $\nu\propto 1/\xi$, or for QCD with $n_{f}$ flavors of fermions

[TABLE]

where the last relation follows from perturbation theory for $\beta$ large enough that the system exhibits asymptotic scaling. In section 7.2 we shall present numerical evidence that this hypothesis may indeed be true. If this is the case, then tuning the HMC algorithm by varying the trajectory length proportionally to the correlation length does not lead to any change in the amplification of rounding errors as we change $\xi$.

7 Monte Carlo Results

We have carried out extensive numerical studies of reversibility errors for (quenched) $SU(3)$ gauge theory and for full QCD with two flavors of dynamical Wilson fermions. Measurements were made principally on $4^{4}$ lattices, but additional data on $8^{4}$ lattices was used for an analysis of finite size effects.

7.1 Numerical Determination of Reversibility Errors

In order to measure the deviation from reversibility of MD trajectories we chose some initial point in fictitious phase space $(U_{i},\pi_{i})$ and followed an MD trajectory of length $\tau$ to the final point $(U_{f},\pi_{f})$; we then reversed the fictitious momenta and followed the backward trajectory of length $\tau$ from $(U_{f},-\pi_{f})$ to $(U^{\prime}_{i},\pi^{\prime}_{i})$. At the end of this backward trajectory we measured the amount by which it failed to return to its starting point by computing two quantities, the change in energy $\Delta\delta H\equiv H(U_{i}^{\prime},\pi_{i}^{\prime})-H(U_{i},\pi_{i})$, and the norm of the change in the gauge field131313The configuration norm $||\Delta\delta U||$ could have been made into a phase space norm by adding $||\Delta\delta\pi||$ to it, but we do not believe that this would make any significant difference in practice (and it would require saving the initial momenta $\pi_{i}$, which are not needed in the usual HMC algorithm but are needed for the L2MC and GHMC algorithms). $||\Delta\delta U||\equiv||U_{i}^{\prime}-U_{i}||$, where ||U||^{2}\equiv\sum_{x,\mu}\sum_{a,b}\Bigl{|}[U_{\mu}(x)]_{ab}\Bigr{|}^{2}, $a$ and $b$ being $SU(3)$ color indices. Observe that both quantities are extensive in the lattice volume, that $\Delta\delta H$ is gauge invariant whereas $||\Delta\delta U||$ is not, and that $||\Delta\delta U||$ cannot be small unless $\Delta\delta H$ is too. These quantities can also be considered as the difference between the change in energy or gauge field over the forward and the backward trajectories, which serves to explain our notation. Both of these quantities vanish identically for reversible trajectories, of course.

7.1.1 $SU(3)$ Gauge Theory

For $SU(3)$ gauge theory an example of a long trajectory is shown in Figure 8. This shows how the energy ${\delta H}$ changes along a long trajectory. The change in energy over the forward trajectory was ${\delta H}=1.2$ (corresponding to an acceptance probability of 30%), which differed from the change in energy over the backward trajectory by $\Delta\delta H=0.36$. This clearly shows that the leapfrog integration scheme tries very hard to conserve energy even though it has wandered far from the true path: indeed, the backward trajectory ends up a distance $||\Delta\delta U||=20$ away from the starting configuration.

This example demonstrates that there can be severe violations of reversibility even though the acceptance rate is quite reasonable, but it should be noted that a trajectory length of $\tau=40$ is unreasonably long for a system whose correlation length $\xi=O(1)$, especially on a $4^{4}$ lattice!

By carrying out such forward and backward trajectories for a variety of step sizes and trajectory lengths we produced Figures 9 and 10.

In the former we plot the logarithm of $|\Delta\delta H|$, $|{\delta H}|$, and $||\Delta\delta U||$ as a function of the integration step size. The differences between measurements made on different equilibrated configurations were very small. The “kinks” in some of the curves are just a consequence of rounding $\tau$ to the nearest integer multiple of $\delta\tau$. The top and bottom graphs clearly show the integration instability “wall” at $\delta\tau\approx 0.6$, which has spread out just as we would expect. The middle graph shows that by the time one has reached the “wall” ${\delta H}=O(10^{3})$, so the integration instabilities are of no practical importance for this system. If one looks carefully at the left end of the bottom graph one sees that $||\Delta\delta U||$ increases a little as $\delta\tau\to 0$. This is because the number of integration steps $n=\tau/\delta\tau$ is increasing, and therefore rounding errors are being amplified by the usual $\sqrt{n}$ factor.

For extremely long trajectories reversibility is violated even for very small values of $\delta\tau$; in order to understand this it is instructive to look at Figure 10, where we plot the logarithm of ${\delta H}$ and $||\Delta\delta U||$ as a function of the trajectory length for two values of $\beta$ and three values of $\delta\tau$. The lower graph immediately shows us that even for very small step sizes the errors are amplified exponentially, with an exponent which is independent of $\delta\tau$ (i.e., the lines for the same $\beta$ value but different $\delta\tau$ are parallel). This is evidence for chaos in the continuous time equations of motion, as suggested by Jansen and Liu [11]. The lines for small step size on the lower graph eventually curve downwards for the largest $\tau$ values because $SU(3)$ is a compact manifold so there is a maximum distance two configurations can be from one another, and this bound is becoming saturated. The upper graph shows that ${\delta H}$ stays small for small $\delta\tau$ even though the rounding errors have been enormously amplified.

7.1.2 QCD with Wilson Fermions

We have carried out the corresponding analysis for full QCD with two flavors of Wilson fermions with configurations taken from two ensembles, one at $\beta=5.1$ and $\kappa=0.16$ which corresponds to quite heavy quarks, and one at $\beta=3.5$ and $\kappa=0.2225$ which is extremely close to $\kappa_{c}$. We took care that both ensembles were in the “confined phase,” if we use thermodynamic terminology to approximately describe the finite size effects on our symmetric lattices.

Figure 11 shows a typical trajectory for the light dynamical fermion system. For this trajectory ${\delta H}=3.9\times 10^{-2}$ corresponding to an acceptance probability of 96%, and $||\Delta\delta U||=1.4\times 10^{-2}$. This trajectory exhibits some typical behavior: the potential energy oscillates smoothly with a period of about $2$ MD time units; the fermionic energy is much more jittery; and the kinetic energy adjusts itself to cancel the other two components leaving a small residual ${\delta H}$ whose magnitude is about a percent of that of its constituent parts. The three constituent contributions to the energy are of similar magnitudes, which indicates that the fermions are indeed very light: with heavier fermions the fermionic energy is much smaller.

From many such trajectories on several equilibrated configurations we produced the data shown in Figures 12 and 13. For simplicity141414The data for the light fermions is similar in form but on a very different $\delta\tau$ scale. The results obtained from all the data will be shown in Figure 14 on page 14. we only show the data for the heavy fermion case in Figure 12. The results are very similar to the pure $SU(3)$ gauge theory results of Figures 9 and 10, especially if one notes that the trajectory lengths are somewhat longer in the present case.

Just as in that case we observe the smeared-out “wall” at which the leapfrog integration scheme becomes unstable, but that the acceptance rate is already completely negligible by the time the wall is reached. We also see that long trajectories with small step sizes exhibit chaotic behavior.

All our data show a clear exponential instability in $||\Delta\delta U||$ as a function of $\tau$, so we fitted the data over the range of $\tau$ for which this exponential behavior was obvious and extracted a characteristic exponent $\nu$ (which we shall henceforth call the Liapunov exponent, although this might not always be quite the correct terminology). In Figure 14 we show the Liapunov exponent as a function of the step size $\delta\tau$ for pure $SU(3)$ gauge theory, and QCD with heavy and light dynamical fermions.

The results show the same qualitative behavior as that of our toy model which was exhibited in Figure 7 on page 7, except that

• •

The “wall” where the integration instability sets in is at a different value of $\delta\tau$ (about $0.6$ for the pure gauge theory and heavy fermion cases, and about $0.1$ for the light fermion case). This probably just reflects the different highest frequencies of the systems.

• •

The exponent is not zero for small $\delta\tau$, but has some fixed non-vanishing value. This is just the evidence for chaotic continuous time dynamics discussed earlier.

The integration instability is easily avoided by choosing a sufficiently small step size, which is forced upon us anyhow as the lattice volume becomes large if we want to maintain a reasonable Metropolis acceptance rate. This, of course, is exactly what we mean by saying that the integration instability is a finite volume effect. Furthermore, it is clear that even a $4^{4}$ lattice is sufficiently large for us to be driven well away from the “wall” for all the parameters we have considered.

We shall therefore turn our attention to the chaotic instability which is present for any value of $\delta\tau$, and in the next section we shall investigate the dependence of $\nu$ upon the parameters of our lattice field theories.

7.2 Parameter Dependence of Liapunov Exponent

We investigated the dependence of the Liapunov exponent $\nu$ on the coupling constant $\beta$ for pure $SU(3)$ gauge theory.151515Somewhat similar results have been reported by Jansen and Liu [12], although their data does is not obviously consistent with ours, especially for large $\beta$. The results are shown in Figure 15 and 16. In the former we plot $\ln||\Delta\delta U||$ against $\tau$ for a variety of values of $\beta$. All of the curves exhibit a very clear exponential instability: for small $\beta$ the relation is almost completely linear on this semi-logarithmic plot, whereas for large $\beta$ there is some initial curvature before a linear region is reached. The slope (i.e., the Liapunov exponent) decreases very rapidly between $\beta=5.4$ and $6.0$, which is just below the finite-size analogue of the finite temperature deconfinement transition on a $8^{4}$ lattice.

In Figure 16 we plot the measured values of the Liapunov exponent $\nu$ as a function of $\beta$ for both $4^{4}$ and $8^{4}$ lattices. For small $\beta$ the lattice theory is in the strong coupling regime and does not obey the asymptotic scaling behavior mandated by the renormalization group equations and the perturbative $\beta$-function. For large $\beta$ the system is in a tiny box, and is thus in the deconfined phase. Therefore we can at best only trust the data in some “scaling window” near $\beta=5.7$. In this region we have fitted the $8^{4}$ data to our suggested asymptotic scaling form of eq. (7) and we find the fit surprisingly good, especially considering that there is only one free parameter.

It is dangerous to rely too much on asymptotic scaling results obtained on such small lattices, so we will content ourselves with the statement that our data is consistent with our hypothesis that $\nu\xi$ is a constant. The value of $\nu$ is also sufficiently small that the amplification of rounding errors for trajectories of length $\tau=\xi$ is a factor of $O(1)$, and thus unimportant.

We do not have a clear understanding of why the value of $\nu$ at large $\beta$ is so large for $4^{4}$ lattices, other than the obvious suggestion that this is a finite size artefact. The data from $8^{4}$ lattices shows a small decrease in $\nu$ for large $\beta$, but this is not convincing evidence that the data will eventually fall onto our hypothetical asymptotic scaling curve.

We can gain a little more understanding of the mechanism by which $\nu$ could decrease as we approach the continuum limit by noting that $\nu$ depends on $\beta$ in two ways: there is an explicit $\beta$-dependence of the equations of motion, and there is an implicit $\beta$-dependence in the equilibrium ensemble of configurations over which $\nu$ is measured. This latter effect is illustrated in Figure 17, where we plot $\ln||\Delta\delta U||$ against $\tau$ for $\beta=5.1$ and $\kappa=0.16$ on a $4^{4}$ lattice for three different configurations, one hot, one cold, and one chosen at random from the equilibrium distribution.161616We verified that the results for the equilibrated configuration did not change significantly when we selected another such configuration. The hot configuration yields a larger Liapunov exponent than the equilibrium one, and the cold configuration is consistent with a power-law dependence of $||\Delta\delta U||$ on $\tau$. This result suggests that we should not expect to get the correct continuum result from lattices which are too small to be in the confined phase (even if we wanted to consider a finite temperature system in the deconfined phase we should at least ensure that the spatial volume is large enough).

In Figure 17 we also show two-parameter fits to the hot and equilibrium trajectories of the form $||\Delta\delta U||\propto e^{\nu\tau}$, and a one-parameter fit to the cold trajectory of the form $||\Delta\delta U||\propto\tau^{3/2}$. These fits are stable if we add more free parameters by fitting to a combination of an exponential and power dependence. The power-law behavior of the cold trajectory may be explained by a combination of a factor of $\sqrt{\tau}$ for the random walk evolution of rounding errors and a linear factor characteristic of the divergence of nearby trajectories for interacting dynamical systems in the stable (non-chaotic) region of phase space. .

7.3 Reversibility and Conjugate Gradient Accuracy

Our final topic is to investigate the effect of inaccurate CG solutions on instabilities and reversibility for dynamical fermion computations.

In Figure 18 we show the dependence of the Liapunov exponent $\nu$ and the Metropolis acceptance rate as a function of the logarithm of the CG residual for a time-symmetric (zero) initial CG vector. We define the CG residual for the approximate solution $x^{\prime}$ of the linear equations $Ax=b$ as $r\equiv||b-Ax^{\prime}||/||b||$. We also show the logarithm of the coefficient $k(r)$ for completeness.

The results clearly show that the choice of residual has no effect for $r<10^{-5}$, and that for $r>10^{-2}$ the acceptance rate is essentially zero. There would seem to be no benefit from finding the solution more accurately than is needed to give a reasonable acceptance rate.

The effect of improving the convergence of the CG algorithm by using a time-asymmetric initial vector is illustrated in Figure 19. For a residual $r=10^{-8}$ the solution vector is sufficiently independent of the starting vector that the choice of initial vector has no effect on the reversibility (or otherwise) of the trajectory. On the other hand, for $r=10^{-5}$ there is a large difference between the time-symmetric start and the intrinsically irreversible time-asymmetric starts.

It is clear that the reduction in the number of CG iterations required to reach a prescribed residual obtained by using an starting vector derived from previous nearby solutions [14] needs to be offset against the need to have a more accurate solution in order to avoid increased reversibility errors. To reach definitive conclusions more extensive studies would have to be done.

8 Conclusions

• •

Instabilities and the concomitant amplification of rounding errors caused by inaccurate numerical integration of the equations of motion should not be a problem if the integration step size is chosen suitably.

• •

The instabilities due to integration errors become unimportant as $V\to\infty$.

• •

We hypothesize that the Liapunov exponent falls as $\nu\propto\xi^{-1}$ for sufficiently large correlation lengths $\xi$, and therefore we can choose trajectories of length $\tau\propto\xi$ which may reduce critical slowing down without suffering from exponentially large amplification of rounding errors.

• •

The measured values of the Liapunov exponent are of $O(1)$, so the exponential amplification factor for trajectories of length $\tau\approx\xi$ are this only $e^{\nu\xi}=O(1)$.

• •

We have found examples with large reversibility errors but reasonable acceptance rate and small $|\Delta\delta H|$. It is always prudent to verify that $||\Delta\delta U||$ is small.

Acknowledgements

We would like to thank Karl Jansen for kindly sending us a draft of Ref. [12].

This research was supported by by the U.S. Department of Energy through Contract Nos. DE-FG05-92ER40742 and DE-FC05-85ER250000.