Table Of Content

1 Cost of Generalised HMC Algorithms for Free Field Theory A. D. Kennedy and Brian Pendleton a aDepartment of Physics and Astronomy, The University of Edinburgh, Edinburgh, EH9 3JZ, Scotland We study analytically the computational cost of the Generalised Hybrid Monte Carlo (GHMC) algorithm for 0 free field theory. We calculate the autocorrelation functions of operators quadratic in the fields, and optimise 0 the GHMC momentum mixing angle, the trajectory length, and the integration stepsize. We show that long 0 2 trajectories are optimal for GHMC, and that standard HMC is much more efficient than algorithms based on the Second Order Langevin (L2MC) or Kramers Equation. We show that contrary to naive expectations HMC n and L2MC have the same volume dependence, but their dynamical critical exponents are z = 1 and z = 3/2 a respectively. J 8 2 1 1. GENERALISED HMC ing the validity of the method, except for some v special values for which the algorithm ceases to 1 The work reported here extends results pre- be ergodic. We may adjust these parameters to 3 sented previously [1,2], to which the reader is minimise the costofaMonteCarlocomputation, 0 referred for details. We begin by recalling that 1 andthemaingoalofthisworkistocarryoutthis thegeneralisedHMC [2]algorithmisconstructed 0 optimisation procedure for free field theory. 0 fromtwo kinds of update for a set of fields φ and HorowitzpointedoutthattheL2MCalgorithm 0 their conjugate momenta π. has the advantage of having a higher acceptance / Molecular Dynamics Monte Carlo: This t rate than HMC for a given step size, but he did a consists of three parts: (1) MD: an approximate l not take in to account that it also requires a - integration of Hamilton’s equations on phase p higheracceptanceratetogetthesameautocorre- space which is exactly area-preserving and re- e lations because the trajectory is reversed at each :h avnerdsibUle(τ;)U(=τ)U:((φ,τπ))−7→1. (φ(2′,)π′A) wmhoemreendteutmU =flip1 MC rejection. It is not obvious a priori which of v − these effects dominates. F :π π. (3) MC: a Metropolis accept/reject Xi 7→− The parameters τ and ϑ may be chosen in- test. dependently from some distributions P (τ) and R r Partial Momentum Refreshment: This a PM(ϑ) for each trajectory. In the following we mixes the Gaussian-distributed momenta π with shall consider various choices for the momentum Gaussian noise ξ: refreshmentdistributionP ,butweshallalways M π′ cosϑ sinϑ π take a fixed value for ϑ. = F (cid:18) ξ′ (cid:19) (cid:18) sinϑ cosϑ (cid:19)· (cid:18) ξ (cid:19) − 2. FREE FIELD THEORY The HMC algorithm [3] is the special case where ϑ= π2. The L2MC/Kramers algorithm [4,5] cor- ConsidZerasystemofharmonicoscillators{φp} respondstochoosingarbitraryϑbutMDMCtra- for p V. The Hamiltonian on phase space is ∈ jectories of a single leapfrog integration step. H = 1 Z π2+ω2φ2 . This describes free 2 p∈ V p p fieldthePoryinm(cid:0)omentums(cid:1)paceifthefrequencies 1.1. Tunable Parameters ω are chosen as p The GHMC algorithm has three free parameters, the trajectory length τ, the momentum d πp mixing angle ϑ, and the integration step size δτ. ωp2 =m2+4 sin2 Lµ. (1) These may be chosen arbitrarily without affect- µX=1 2 3. AUTOCORRELATION FUNCTIONS obtainedbyminimisingthecost,thatisbychoos- C C ing τ¯ so as to satisfy d /dτ¯=d /dϑ=0. 3.1. Simple Markov Processes We wish to compare the performance of the Let(φ ,φ ,...,φ ) be a sequence of field con- 1 2 N HMC, L2MC and GHMC algorithms for one di- figurations generated by an equilibrated ergodic mensional free field theory. To do this we com- Markovprocess,andlet Ω(φ) denotetheexpec- h i parethecostofgeneratingastatisticallyindepen- tation value of some connected operator Ω. We dent measurement of the magnetic susceptibility maydefineanunbiasedestimator Ω¯ overthefinite M2, choosing the optimal values for the angle ϑ sequenceofconfigurationsbyΩ¯ 1 N Ω(φ ), c ≡ N t=1 t andtheaveragetrajectorylengthτ¯. Wecanmin- Ω(φt+ℓ)ΩP(φt) imise the cost with respect to ϑ without having Asusual,wedefineC (ℓ) asthe Ω ≡ (cid:10) Ω(φ)2 (cid:11) to specify the form of the refresh distribution. autocorrelation function for Ω.(cid:10)The (cid:11)variance of The next step is to minimise the cost with re- the estimator Ω¯ is spect to the average trajectory length ξ = ωτ¯. Strictly speaking we should note that the accep- Ω(φ)2 N Ω¯2 = 1+2A 1+O exp , tance probability P¯ is a function of τ¯, but to h i { Ω}(cid:10) N (cid:11)(cid:20) (cid:18) N (cid:19)(cid:21) acc a good approximation we may assume that P acc ∞ depends only upon the integration step size δτ where A C (ℓ) is the integrated auto- Ω ≡ ℓ=1 Ω except in the case of very short trajectories. correlation funPction for the operator Ω and N exp is the exponential autocorrelation time. This re- 4.1. Exponentially Distributed Trajectory sult tells us that on average 1+2A correlated Ω Lengths measurements are needed to reduce the variance Toproceedfurtherweneedtochooseaspecific bythesameamountasasingletrulyindependent form for the momentum refresh distribution. In measurement. this section we will present results for the case of exponentially distributed trajectory lengths, 3.2. Autocorrelation Functions for P (τ)=re−rτ where the parameter r is just the Quadratic Operators R inverse average trajectory length r =1/τ¯. Inordertocarryoutthesecalculationswemake The cost at the point (c cosϑ ,ξ ) is the simplifying assumption that the acceptance opt ≡ opt opt probabilitymin 1,e−δH for eachtrajectorymay C exp(δτ,ϑ ,ξ ) = bereplacedbyit(cid:0)svaluea(cid:1)veragedoverphasespace Mc2 opt opt P min 1,e−δH ; we neglect correlations acc ≡ (7P¯ 3P¯2 4)ξ 2c 3+1 c 2 between(cid:10)succ(cid:0)essive tr(cid:1)a(cid:11)jectories. Including such  +a(cc−2P¯ a+cc−1)c o+pt (2oPp¯t 1−)c op3t  correlations leads to seemingly intractable com- − acc opt acc− opt +(P¯2 5P¯ +4)c ξ 2 plications. It is not obvious that our assumption  +( P¯ +4a)ccξ−2+ac(c 4+3oPp¯t op)tc 2ξ 2  correspondsto anysystematic approximationex-  − acc opt − acc opt opt . cept, of course, that it is valid when δH =0. P¯ δτω(P¯ 1)c 3ξ Details of the calculation of P for leapfrog  +ξ acPc¯ δτωacc−ξ P¯optδτoωpct 2  integration are published elsewheraecc[1,2,6]. oPp¯t aδccτω( −1+opPt¯ ac)cc ξ opt  − acc − acc opt opt  4. COMPARISON OF COSTS Thissolutionisafunctionofδτ andP¯ whichare acc not independent variables, and using the results Ifwe makethe reasonableassumptionthat the for P¯ (τ,δτ) [1,2] we can compute the cost as a acc cost of the computation is proportional to the function of P¯ as shown in Figure 1. acc total fictitious (MD) time for which we have to C integrate Hamilton’s equations, then the cost 4.1.1. HMC per independent configuration is proportional to The Hybrid Monte Carlo algorithm corre- (1+2A )τ¯/δτ withτ¯denotingtheaveragelength sponds to setting ϑ = π/2, and we find that the Ω of a trajectory. The optimal trajectory length is optimal trajectory length in this case is ξ = opt 3 27600 wherethescalingvariablex=Vδτ6isverysmall. GHMC 27400 We may then express c and P¯ as power se- LH2MMCC Cost 2277020000 ries in x, keeping only tohpet first feawcc terms. From 107 26800 theserelationswefindthattheminimumcostfor 26600 106 3.5 3.6 3.7- log10(1-Pacc)3.8 3.9 4.0 L2MC is ost105 C104 C = 10 1/4V5/4m−3/2. 103 opt (cid:18) π (cid:19) 1.5 This result tells us that not only does the tuned θopt1.0 L2MC algorithm have a dynamical critical expo- 0.5 0.0 nent z = 3/2, but also it has a volume depen- denceofexactlythesameformasHMC[1,2]. We 0.4 τξ/opt0.2 mthaaynutnhdeenrsatiavnedVw7h/6ymth−i1sbbyehtahveiofuorlloowcciunrgssriamthpeler 0.0 argument. 0.0 0.2 0.4 0.6 0.8 1.0 P IfP¯ <1thenthesystemwillcarryoutaran- acc acc domwalkbackwardsandforwardsalongatrajec- tory because the momentum, and thus the direc- tion of travel, must be reversedupon a Metropo- Figure1. CostasafunctionofaverageMetropolis lisrejection. Asimplemindedanalysisisthatthe acceptance rate for the GHMC algorithm com- averagetime betweenrejections must be O(1/m) pared to HMC and L2MC for free field theory. in order to achieve z = 1. This time is approxi- The operator under consideration is the “mag- mately neticsusceptibility”,i.e.,theconnectedquadratic ∞ P¯ δτ 1 operator depending only on the lowest frequency P¯n (1 P¯ )nδτ = acc = . acc − acc 1 P¯ m mode. The corresponding parameters, the mo- nX=0 − acc mentum mixing angle ϑ and the average tra- opt For small δτ we have 1 P¯ = erf√kVδτ6 jectorylengthmeasuredasafractionofthecorre- √Vδτ6 wherek isaconst−anta,ccandhencewemu∝st lationlength τ /ξ are also shown,all as a func- opt scale δτ so as to keep Vδτ4/m2 fixed. Since the tion of the acceptance rate P¯ . The inset graph acc L2MC algorithm has a naive dynamical critical showstheregionwheretheacceptancerateisvery exponent z = 1, this means that the cost should closetounitywhichiswheretheL2MCalgorithm varyasC V(Vδτ4m−2)1/4/mδτ =V5/4m−3/2. has its minimum cost. ∝ REFERENCES 1/ 4 P¯ , corresponding to a cost 1. A. D. Kennedy and B. J. Pendleton, Nucl. − acc p Phys. B (Proc. Suppl.) 20 (1991) 118. C = 2 4−P¯acc. 2. A.D.Kennedy,R.G.Edwards,H.Mino,and opt pP¯accδτω B.J.Pendleton, Nucl.Phys.B(Proc.Suppl.) This is also shown in Figure 1. 47 (1995) 781. 3. S. Duane, A. D. Kennedy, B. J. Pendleton, 4.2. Fixed Length Trajectories andD. Roweth,Phys.Lett. 195B(1987)216. Forfixedlengthtrajectoriesweshallonlyanal- 4. A.M.Horowitz, Phys.Lett.B268(1991)247. yse the case of L2MC for which the trajectory 5. M. Beccaria and G. Curci, Phys. Rev. D49 length ξ = ωδτ. In this case the value of ϑ opt (1994) 2578. and the corresponding cost are also plotted in 6. A.D.KennedyandB.J.Pendleton,inprepa- Figure 1. From this figure it is clear that the ration. minimum cost occurs for P¯ very close to unity, acc