Progress in Lattice Field Theory Algorithms

A. D. Kennedy의 논문 "Progress in Lattice Field Theory Algorithms"는 최근 lattice field theory 알고리즘 개발에 대한 한 해례를 정리한다. 논문의 주요 내용은 다음과 같다.

1. **Multicanonical Algorithm**: 이 알고리즘은 lattice field theory에서 사용되는 새로운 방법으로, 원래의 문제를 분해하고 각 부분 문제를 다른 확률 분포를 사용하여 해결하는 방식이다. 이 방법을 통해 원래의 문제에 대한 통계적 오류를 줄일 수 있다.
2. **Importance Sampling**: Importance sampling은 Monte Carlo 추정에서 유용한 기법으로, 특정 함수에 대한 기대값을 계산할 때 분포를 변경하여 오류를 줄일 수 있다. 논문에서는 optimal importance sampling 방법을 소개한다.
3. **Canonical Distribution**: Canonical distribution은 field theory와 통계 물리학에서 가장 자주 사용하는 확률 분포로, 엔트로피 함수 S를 최대화하는 문제이다. 논문에서는 canonical distribution에 대한 구체적인 예를 제시한다.
4. **Multicanonical Distribution**: Multicanonical distribution은 canonical distribution과 operator의 연산을 함께 고려하여 오류를 줄이는 방법이다. 논문에서는 multicanonical distribution을 사용한 사례와 이에 대한 개선 방안을 제시한다.

논문의 결과는 다음과 같다:

1. **Autocorrelation Time**: Autocorrelation time은 Monte Carlo 알고리즘에서 반복 시간의 평균을 계산할 때 발생하는 오류를 측정하는 지표이다. 논문에서는 multicanonical distribution을 사용하여 autocorrelation time을 줄일 수 있는 방법을 제시한다.
2. **Surface Energy**: Surface energy는 lattice field theory에서 중요한 물리학적 개념으로, 원래의 문제를 분해하고 각 부분 문제를 다른 확률 분포를 사용하여 해결하는 방식이다.

논문은 lattice field theory 알고리즘 개발에 대한 최근 연구결과를 정리하고, multicanonical algorithm을 포함한 다양한 방법을 소개한다. 논문의 결과는 lattice field theory 분야의 발전에 기여할 것으로 보인다.

Progress in Lattice Field Theory Algorithms

arXiv:hep-lat/9212017v1 14 Dec 1992Progress in Lattice Field Theory AlgorithmsA. D. Kennedy∗Supercomputer Computations Research InstituteFlorida State UniversityTallahassee FL 32306–4052, U.S.A.October 30, 1992SCRI preprint FSU–SCRI–92–160CRS4 preprint CRS4–PARCOMP–92–1hep-lat/9212017AbstractI present a summary of recent algorithmic developments for latticefield theories.

In particular I give a pedagogical introduction to thenew Multicanonical algorithm, and discuss the relation between theHybrid Overrelaxation and Hybrid Monte Carlo algorithms.I alsoattempt to clarify the rˆole of the dynamical critical exponent z andits connection with “computational cost.”∗This research was supported by the Florida State University Supercomputer Compu-tations Research Institute which is partially funded by the U.S. Department of Energythrough contract #DE–FC05–85ER250000. This research has been carried out with thepartial support of the Sardinia Regional Authorities.

I would like to thank the Universityof Edinburgh and CRS4 (Cagliari) for their hospitality while this review was being written.1

1INTRODUCTIONIn this review I shall concentrate on two areas in which there has been con-siderable work and significant progress during the last year. The first is theMulticanonical algorithm which, although the underlying ideas have beenknown for some time, has been successfully applied to studying many prop-erties of first-order transitions in lattice field theories.

I shall attempt togive a somewhat pedagogical introduction to this method. The second isthe Hybrid Overrelaxation algorithm; this method has been used for someyears, but recent analysis of its performance for the Gaussian model and ofits relationship with the Hybrid Monte Carlo algorithm are illuminating.Unfortunately there has been little activity or progress on what I considerthe two outstanding algorithmic challenges facing lattice field theory today:dynamical fermions and complex actions.

For the former the algorithms seemto work quite well, albeit slowly, whereas for the latter the situation is muchworse.There has been a lot of work on multigrid and cluster methods. Both havebeen shown to work well in some simple models (usually in two dimensions),but their efficacy has not yet been demonstrated for four dimensional non-abelian gauge theories.

Since these methods were discussed in considerabledetail in previous lattice conferences I shall only survey the progress verybriefly here.I shall also try to clarify the situation regarding what is known aboutcritical slowing down for the Hybrid Monte Carlo algorithm.2MULTICANONICAL METHODSAlthough the application of the multicanonical algorithm to lattice field the-ory is new [1, 2, 3], similar techniques have been used for some time [4, 5, 6].For a recent review of developments in this subject see [7].2.1Marginal DistributionsA statistical mechanical system or a quantum field theory may be describednot only in terms of their detailed microscopic states, but also in terms ofa set of macroscopic parameters (e.g., energy, density, magnetization) which2

suffice to specify all of its thermodynamic properties. We shall denote thespace of microstates by M, which has “volume” Z ≡RM dφ with respect tothe natural measure dφ, and we shall write ⟨·⟩f to indicate an average overM with respect to the measure ∝dφ f(φ).

Similarly, we shall call the spaceof macrostates O, with volume Z′ ≡RO dE and averages over this space withrespect to the measure ∝dE g(E) will be written as ⟨⟨·⟩⟩g. If Ωis some ther-modynamic observable which depends only on some macroscopic parametersE then the situation may be summarized by the following equationME−→OΩ−→R.

(1)Microscopic and macroscopic averages are related in a trivial way:1⟨Ω◦E⟩= 1ZZM dφ Ω(E(φ))(2)= 1ZZM dφZO dE δ(E −E(φ))Ω(E)(3)=ZO dE ρ(E)Ω(E) = ⟨⟨Ω⟩⟩ρ(4)where we have introduced the density of statesρ(E) ≡1ZZM dφ δ(E −E(φ)). (5)The probability distribution generated by ρ(E) is the marginal distributionon O, meaning that it is obtained from the full distribution over the muchlarger space M by averaging over all the other variables.2.2Importance SamplingImportance sampling is a fundamental technique to reduce the statisticalerrors in a Monte Carlo computation.Suppose we wish to measure theexpectation value of some quantity f : M →R: to this end we introduce anestimatorEµ"fµ#= 1TTXt=1f(φt)µ(φt)(6)1We use the notation Ω◦E to indicate functional composition: (Ω◦E)(φ) = ΩE(φ).3

where samples {φt} are chosen from the distribution µ(φt). It is easy to seethat the expectation value ⟨Eµ [f/µ]⟩µ = ⟨f⟩is independent of the choice ofdistribution µ; in fact the central limit theorem tells us the stronger resultthatEµ"fµ#= ⟨f⟩+ OsVµ [f/µ]T,(7)where the variance of a single sample isVµ"fµ#=*fµ −*fµ+µ2+µ(8)=*f 2µ+−⟨f⟩2 .

(9)Observe that we can reduce the error not only by increasing T but also bychoosing µ so as to make Vµ [f/µ] smaller.We may find the optimal importance sampling by a simple variationalcalculationδδµVµ"fµ#+λ ⟨µ⟩µ=µopt==−* f 2µ2opt++ λ = 0.The solution of this equation gives the optimal probability distributionµopt = |f|⟨|f|⟩;(10)this gives a variance per sample ofVµopt" fµopt#= ⟨|f|⟩2 −⟨f⟩2 ,(11)which vanishes ifff has the same sign ∀φ. While we cannot often achieve thisideal situation, generating a probability distribution which approximates theoperator being measured can greatly reduce the amount of computer timenecessary to get a reliable measurement of the operator’s expectation value.4

2.3Canonical DistributionIn field theory as well as for statistical systems we are most often interestedfinding the expectation values of operators with respect to the canonical dis-tribution. For statistical systems this distribution is obtained by maximizingthe entropy S ≡−⟨log P⟩P with respect to the probability distribution Psubject to the constraints that the ensemble average is at a given point in O:⟨1⟩P = 1, ⟨E⟩P = E ∈O.δδP {S + λ ⟨1⟩P −β ⟨E⟩P}P =Pc = 0(12)whose solution isPc(φ) = e−βE(φ)Z(β) ,(13)with the partition function Z and free energy F given byZ(β) =ZM dφ e−βE(φ) = e−βF (β),(14)andE = −∂∂β log Z(β).

(15)In terms of macroscopic averages the canonical distribution of Eq. (13) is⟨Ω◦E⟩e−βE = ⟨⟨Ω⟩⟩ρ(E)e−βE .

(16)Taking configurations from the canonical distribution Pc gives good impor-tance sampling if Ω(E) ≈1 for all E in the subspace of O of interest.It is important to note that we do not know Z(β) a priori, so using anyother importance sampling requires the computation of ratios of estimators.A ratio of unbiased estimators is not an unbiased estimator for the ratio, socare must be taken to avoid systematic errors; as we shall see, this ought tobe done for multicanonical computations.2.4Multicanonical DistributionIf Ω(E) has exponentially large peaks then important contributions maycome from regions where e−βEρ(E) is small.This situation is illustrated5

Figure 1: A situation in which the multicanonical method is useful. Thesolid line is the canonical distribution and the dashed line is the operator ofinterest.6

schematically in figure 1. A compromise to get good importance samplingsimultaneously for Ωe−βE and e−βE is to generate configurations with a uni-form density on macrostate space O.

This is the multicanonical distributionPMC(φ) ∝1ρ(E(φ)). (17)It leads to ⟨Ω◦E⟩1ρ◦E = ⟨⟨Ω⟩⟩and⟨Ω◦E⟩e−βE=Dρ ◦E e−βE Ω◦EE1ρ◦E⟨ρ ◦E Ω◦E⟩1ρ◦E(18)=DDρe−βEΩEE⟨ρe−βE⟩.

(19)2.5Algorithms to Generate theMulticanonical DistributionIf ˜ρ(E) ≈ρ(E) then we can readily generate an approximate multicanonicaldistribution˜PMC(φ) ∝1˜ρ(E(φ))(20)for which⟨Ω◦E⟩e−βE=D˜ρ ◦E e−βE Ω◦EE1˜ρ◦E⟨˜ρ ◦E Ω◦E⟩1˜ρ◦E(21)=DD˜ρe−βEΩEEρ/˜ρ⟨˜ρe−βE⟩ρ/˜ρ. (22)This is cheap if all the macroscopic parameters E are local operators onφ ∈M.

We can use a Metropolis algorithm with acceptance probabilityP(φ →φ′)=min 1,˜PMC(φ′)˜PMC(φ)! (23)=min 1, ˜ρ(E(φ))˜ρ(E(φ′))!

(24)=min1, e−∆log ˜ρ◦E. (25)7

Presumably it would be better to use Hybrid Monte Carlo or Hybrid Over-relaxation algorithms if they have z ≈1.For low-dimensional O we may represent the approximate spectral densitylog ˜ρ(E) by binning. Berg and Neuhaus used a piecewise linear interpolation(canonical within each bin)log ˜ρ(E) =NXj=1χj(E)(βjE + αj)(26)[χj(E) is the characteristic function of bin j: it has the value 1 if E isin the bin and 0 otherwise], whereas Marinari and Parisi use a stochasticsuperposition of equiprobable linear interpolations (each term canonical)˜ρ(E) =NXj=1Pj(βjE + αj)(27)[Pj is the probability of choosing term j from the “sum”].

The former isillustrated in figure 2.2.6ApplicationsThere have been many calculations making use of the multicanonical methodduring the last year, including [8, 9, 10, 11, 12, 13, 14, 15, 16]. We shall justconsider two very simple examples in which it is effective:• The surface energy at a first order phase transition, and• The autocorrelation time (tunnelling time) at a first order phase tran-sition.The surface energy may be defined byPmaxPmin= 2FsurfaceLD−1(L →∞),(28)where Pmin and Pmax are illustrated in figure 3.

In order to relate this to ourformal analysis the appropriate operators may be expressed asPmaxmin = limγ→±∞1γ logDDeγP EE;(29)8

Figure 2: A piecewise linear approximation to log ρ.9

Figure 3: Schematic surface energy calculation.10

these operators are shown (for some large but finite value of γ) with dashedlines in figure 3. The lack of overlap between the operator giving Pmin andthe canonical distribution is obviously an example of the situation shown infigure 1.

For this reason it is clear why this operator will be sampled muchbetter by the multicanonical algorithm, which will generate configurations ofenergy E uniformly distributed over the range of interest.The autocorrelation time (tunnelling time) for local Monte Carlo algo-rithms at a first order phase transition is notoriously long. The problem hereis that any local algorithm requires the system to pass through the minimumseparating the two phases (figure 4), so the tunnelling time is approximatelyproportional to Pmin.

Since the logarithm of the density of states is an ex-tensive quantitylog ρ(E) ∝LD(30)where LD is the lattice volume, it follows thatτA ∝Pmin ∝e−∆LD,(31)where ∆is the error in the approximation for log ρ(E) at any fixed volume.For a true multicanonical distribution ∆= 0, and for an approximate one itis at least much smaller than the canonical barrier, as shown in figure 4.2.7Open QuestionsWhile the multicanonical method works very well in practice, there are stillseveral interesting algorithmic questions about it which ought to be ad-dressed.• For a fixed ˜ρ(E) the autocorrelation time τA is still exponential in LD,but with a much smaller exponent.This follows simply from equa-tion (31).• We can measure ρ(E) during the course of a multicanonical computa-tion by counting the number of configurations landing in each bin, andit can be used to improve the approximation ˜ρ(E).• Just as for simulated annealing, however, no one has given an algorithmfor evolving ˜ρ(E). In order to do this one would have to address thefollowing issues:11

Figure 4: Another situation in which the multicanonical method is useful;the canonical distribution at a first-order phase transition.12

– How do we distinguish a statistically significant difference betweenρ and ˜ρ from fluctuations?– How do we ensure stability, and that ˜ρ →ρ?– The number of bins must presumably grow as LD to avoid an expo-nential autocorrelation time, yet there must be sufficient statisticsin each bin to give a significant estimate for ˜ρ.– Is there an algorithm giving a power-law volume dependence forτA?– Do we care? In practical simulations we might not mind havingan exponential algorithm provided that the exponent was smallenough.

We might even prefer it to a polynomial algorithm witha large coefficient.3HYBRID OVERRELAXATIONThe name “Hybrid Overrelaxation” appears to have been coined by Wolff[17] who recently analyzed the dynamical critical behaviour of this algorithmin the Gaussian model. The algorithm, however, has been used in a varietyof models over the past few years [18, 19, 20, 21, 22, 23, 24, 25, 26, 27].

Itbelongs to a class of overrelaxed algorithms introduced by Adler [28] withdynamical critical exponent z ≈1 for the Gaussian model [29, 30, 31].3.1Adler OverrelaxationConsider the Gaussian model defined by the free field actionS(φ) = 12XxDXµ=1(∂µφ(x))2 + m2φ(x)2. (32)A single-site Adler overrelaxation (AOR) update [28] with parameter ζ re-places φ(x) byφ′(x) = (1 −ζ)φ(x) + ζFω2 +qζ(2 −ζ)ωη,(33)13

where ω2 ≡2D + m2 is the square of the highest frequency in the spectrum,the “force” on φ(x) due to its neighbours is F ≡P|x−y|=1 φ(y), and η is aGaussian-distributed random number.3.2z = 1 for Gaussian ModelThe lattice may be updated using a checkerboard scheme, alternating theupdate of all even and all odd sites. This is just a consequence of the localityof the action.

This is also the basis of Wolff’s [17] analysis of the dynamicalcritical behaviour of overrelaxation in the Gaussian model: the fields on theeven and odd sublattices can Fourier transformed, and equation (33) becomesblock diagonal in momentum space. Determining the autocorrelations of themethod thus reduce to studying the eigenvalues of a 2 × 2 matrix.Let us now consider two special cases of the Adler overrelaxation algo-rithm, corresponding to different values of ζ.

For ζ = 1 this is the heatbath(HB) algorithm, because φ′(x) does not depend upon φ(x). The exponentialautocorrelation time isτHB =Dm2 + p2(p2 →0),(34)which corresponds to z = 2.

If we adjust ζ so as to minimize the autocorre-lation time we find thatζ = 2 −2m√D + O(m2),τAOR ≈√D2m ;(35)which gives z = 1.3.3Hybrid OverrelaxationFor ζ = 2 the new field value, φ′(x), does not depend on the noise η, so theupdate is not ergodic. The Hybrid Overrelaxation (HOR) algorithm curesthis by alternating N AOR steps with a single HB step.

Minimizing theautocorrelation time requires increasing N ∝1/m, which leads to τHOR = 1.Just as pointed out by Weingarten and Mackenzie [32, 33] for HybridMonte Carlo (HMC), N should be varied randomly around this value toavoid accidental non-ergodicity; although this is probably only a problem forthe Gaussian model and not for interacting field theories.14

3.4Overrelaxation and HybridMonte CarloI now want to show that the HOR algorithm — and the AOR algorithmtoo — are closely related to the HMC algorithm. To this end lets us con-sider the HMC algorithm applied to the Gaussian model: we introduce theHamiltonianH(π, φ) = 12π2 + S(φ)(36)on “fictitious” phase space.

The corresponding equation of motion is¨φx = −ω2φx + F,(37)whose solution in terms of the Gaussian distibuted random initial momentumπx and the initial field value φx isφx(t) = φx cos ωt + 1 −cos ωtω2F + πxω sin ωt,(38)which is exactly the AOR update considered before if we identify ζ ≡1 −cos ωt. If we use this exact solution of the equations of motion to generatecandidate configurations for HMC then the acceptance rate will be unity.For interacting field theories HMC provides an exact algorithm for any valueof the overrelaxation parameter ζ simply by dropping the interaction termsin the action and solving the equations of motion exactly for the resultingGaussian model.Let us summarize the differences between this variant of the HMC algo-rithm and the usual one.

For “conventional” HMC:• All sites are updated at once;• Each trajectory is of length τ ∝1/m and is a sequence of many stepsof length δτ;• There is a global Metropolis accept/reject step.For “local” HMC [34, 35, 36]:• Even and odd updates are alternated;• Each trajectory is of length τ ≈π/ω = O(1);• There is a local Metropolis accept/reject step.15

3.5Leapfrog vs. Free Field GuidanceIt is important to realize that there are two separate issues; one is whethera global accept/reject step is used, the other is whether an exact solution tothe free-field equations of motion or an approximate (leapfrog) solution tothe true equations of motion is used.• Local Metropolis acceptance/rejection:– Useful only for local (bosonic) theories;– Acceptance rate does not depend on the lattice volume.• Free-field instead of leapfrog guidance:– Leapfrog has δτ errors so Pacc < 1 even for free field theory;– Free field guidance has errors of order λ for interacting theories,whereas leapfrog has errors of order λδτ 2, where λ is the couplingconstant of the interaction part of the action.4DYNAMICAL CRITICALEXPONENTSLet us begin this topic by recalling the definition of the dynamical criticalexponent z. It relates a “dynamical” property of the algorithm generatingfield configurations, the autocorrelation time τA (either the exponentionalautocorrelation time or an integrated autocorrelation time will do), and a“static” property of the underlying field theory, the correlation length ξ:τA ∝ξz(ξ →∞).

(39)One of the main “selling points” of algorithms is that they reduce (z ≈1)or even eliminate (z ≈0) critical slowing down from the value (z ≈2)characteristic of “random walk” methods. In general terms it is possible toreduce z from 2 to 1 by using just the right amount of randomness in thealgorithm (see sections 3.2 and 4.2), whereas further reduction of z seems[31] to require an algorithm with more specific “knowledge” of the dynamicsof the theory (section 5).16

Although the theoretical analysis of critical slowing down is usually lim-ited to the Gaussian model, this is relevant because an algorithm is oftenbetter just because it is not quite so dumb in dealing with free fields. Howclose real theories are to free fields is an empirical question: for some in-teracting field theories the values of z have been determined empirically[37, 38, 39, 40, 41].4.1Do we care too much about z?The characterization of algorithms by their dynamical critical exponent hasperhaps been somewhat over-emphasized.

Care is required not only becauseour computations are not always peformed at parameters where the asymp-totic behaviour of the algorithms has set in, but also because the cost of acomputation is not just given by the autocorrelation time.We want to study the continuum limit (critical behaviour) of some latticemodel in a large enough box (thermodynamic limit). In order that the sys-tematic errors are under control we need to match both the short and longdistance behaviour of the lattice regularization to some analytic form:• For asymptotically free theories we hope we can match the ξ/a →∞scaling to perturbation theory.• Non-perturbative effects fall offexponentially fast in UV asymptopia.• For non asymptotically free theories it is unclear how we can verifythat a is small enough.• We hope to match results using finite size scaling for L/ξ →∞.We only need to carry out numerical computations for an essentially finiterange of lattice spacings.

Good asymptotic dynamical scaling behaviour ofan algorithm is useful only if• scaling sets in for the lattice sizes we are interested in,• the coefficient in equation (39) is small.4.2Critical Slowing Down and Finite17

Size ScalingOur definition of z (equation (39)) assumes L ≫ξ, which is usually thecase for lattice gauge theory computations. If we want to study finite sizescaling computations are sometimes done in the regime L ≪ξ, and thedynamical critical exponent is defined by τA ∝Lz [42].

It is not clear thatthe definitions are equivalent if some parameters have to be tuned to valuesdepending upon ξ, e.g., τ ∝ξ in HMC. For HMC in the Gaussian modelwith L ≫ξ, z = 2 if τ =constant, and z = 1 if τ ∝ξ [43, 44].

If ξ ≫L thenz = 2 for τ =constant, but I know of no analysis for any other cases.4.3Computational CostFor the global HMC algorithm z is not the whole story; the computationalcost explicitly depends on the lattice volume even at fixed ξ [45, 46]Tcomp ∝V 5/4 = L5D/4(40)because the integration step size δτ has to be reduced so as to keep theMetropolis acceptance rate constant. It must be remembered that Tcomp ∝Vis true even for local algorithms (one has to look at every site!).

For higher-order leapfrog schemes this cost can be reduced to Tcomp ∝V 1+ǫ (an interest-ing new way of understanding these higher-order algorithms is presented in[47]).Before leaving this subject I would like to mention a couple of relatedtopics. First, some interesting results on the scaling behaviour of HMC inthe presence of dynamical fermions are presented in [48, 49].

Second, aninteresting variant of HMC which adds a small amount of randomness to themomenta at each integration step (instead of a complete momentum refresh-ment after an entire trajectory) was found by Horowitz [50]. Unfortunately,the deterministic part of the momenta has to have its sign flipped after everystep in order to satisfy the detailed balance condition, and this means themodified algorithm does not perform any better (although a detailed proofhas not yet been presented).5MULTIGRID AND18

CLUSTER METHODSI mentioned previously that algorithms which have z < 1 appear to requiresome knowledge of the detailed dynamics of the model being simulated. Forsome spin models this is achieved in a non-trivial manner by cluster algo-rithms, but for lattice gauge theories like QCD all attempts in this directionare motivated by free field theory dynamics.

This is not a bad idea, becausesuch theories are asymptotically free. There has not been much activity inthe area of Fourier acceleration, but there has been a great deal of work onMultigrid methods.

Since such methods were reviewed in detail last year Ishall merely survey work in this area extremely briefly.Grabenstein and Pinn have presented a formalism for computing the ac-ceptance rate for Multigrid Monte Carlo [51, 52]. There have been numerousapplications and tests of multigrid methods, including those of Hasenbuschand Meyer [53, 54, 25] (strictly speaking these are really “unigrid” algo-rithms), Laursen and Vink [55], and Edwards et al.

[56].Multigrid methods are also under active study as a means of inverting thelattice Dirac operator more efficiently [57, 58, 59, 60, 61]. Another methodused for this purpose has been proposed by Vyas [62, 63].

Whether thesetechniques can significantly outperform the conjugate gradient algorithm isnot yet settled.Multigrid methods have also been used to expedite gauge fixing [64] (seealso the work of van der Sijs [65]).Cluster algorithms have also been under active study.An interestingway of understanding them has been suggested by Wolff[66]. Their area ofapplicability is still mainly for spin models of various kinds [67, 68] wherethey are extremely successful and are often the clear method of choice.

Nev-ertheless, their use is limited to models whose spin manifold possesses aninvolutive isometry whose fixed point manifold has codimension one (i.e., adiscrete quotient of a product of spheres), as was shown by Sokal and collab-orators [69]. B. Bunk [70] describes a cluster algorithm which is applicableto Z(2) gauge theory.The parallelism inherent in FFT, multigrid, and cluster algorithms ismore complex than the simple grid-like structure needed for the algorithmsdiscussed in previous sections (they require only nearest-neighbour commu-nication with infrequent global summations).

The implementation of clusteralgorithms on parallel and vector computers has been discussed in several19

papers [71, 72, 73]. In the future it will be interesting to characterize thesealgorithms in terms of the type of communications network required to im-plement them efficiently on a parallel computer (e.g., combining networks,infinite dimensional grids, etc.

).6OTHER ALGORITHMSThere are, of course, numerous interesting results which I have not had timeor space to discuss in detail. I do wish to mention at least a few of them here,so that the interested reader will be able to refer to the original literature.An interesting algorithm called “microcanonical cluster Monte Carlo” wasintroduced by Creutz [74].

This algorithm allows a discrete spin variable tointeract with a heat bath (provided by a family of “demons”), and for thedemon energies to be refreshed occasionally. In some ways this method is adiscrete analogue of the HMC method.Sloan, Kusnezov, and Bulgac [41] introduce a “chaotic molecular dynam-ics” algorithm, which couples chaotic degrees of freedom to a system to ensureergodicity (instead of, or perhaps as well as, refreshing the fictitious momentaof the Hybrid algorithm).

Neal [75] suggests a new procedure which may serveto increase the acceptance rate of the HMC algorithm. Finally, Bhanot andAdler [76] describe a parallel algorithm for updating spin models.References[1] Bernd A. Berg and Thomas Neuhaus.Multicanonical algorithms forfirst order phase transitions.

Phys. Lett., B267:249–253, 1991.

[2] Bernd A. Berg and Thomas Neuhaus. Multicanonical ensemble: a newapproach to simulate first order phase transitions.Phys.

Rev. Lett.,68:9–12, 1992.

[3] Enzo Marinari and Georgio Parisi. Simulated Tempering: A new MonteCarlo scheme.

Europhys. Lett., 19:451–458, 1992.

[4] G. M. Torrie and J. P. Valleau. J. Comp.

Phys., 23:187, 1977. [5] I. S. Graham and J. P. Valleau.

J. Chem. Phys., 94:7894, 1990.20

[6] J. P. Valleau. J. Comp.

Phys., 96:193, 1991. [7] Bernd Berg.

The multicanonical ensemble: a new approach to computersimulations. Technical report, FSU, August 1992.

FSU–HEP–92–0817. [8] Ulli Hansmann Bernd A. Berg and Thomas Neuhaus.

Recent results ofmultimagnetical simulations of the Ising model. Technical report, SCRI,July 1992.

FSU–SCRI–92C–87. [9] Bernd A. Berg and Tarik Celik.

A new approach to spin glass simula-tions. Technical report, SCRI, April 1992.

FSU–SCRI–92–58. [10] Bernd A. Berg.Multicanonical Potts model simulations.Technicalreport, Universit¨at Bielefeld, 1991.

BI–TP–91–21. [11] Bernd A. Berg, Tarik Celik, and Ulrich Hansmann.

Multicanonical studyof the 3-d Ising spin glass. Technical report, SCRI, September 1992.FSU–SCRI–92–121.

[12] Bernd A. Berg and Tarik Celik. Multicanonical spin glass simulations.Technical report, July 1992.

PRINT–92–0295. [13] Wolfhard Janke, Bernd A. Berg, and Mohammad Katoot.

Monte Carlocalculation of the surface free energy for the 2-d seven state Potts modeland an estimate for 4-d SU(3) gauge theory. Technical report, SCRI,February 1992.

FSU–SCRI–92–40. [14] Alain Billoire.

Numerical study of finite size scaling for first order phasetransitions. Technical report, Saclay, June 1992.

SPHT–92–097. [15] B. Grossmann, M. L. Laursen, T. Trappenberg, and U.-J.

Wiese. A mul-ticanonical algorithm and the surface free energy in SU(3) pure gaugetheory.

Technical report, HLRZ, 1992. HLRZ–92–31.

[16] B. Grossmann, M. L. Laursen, T. Trappenberg, and U.-J. Wiese.

Theconfined–deconfined interface tension and the spectrum of the transfermatrix. Technical report, HLRZ, 1992.

HLRZ–92–47. [17] Ulli Wolff.

Dynamics of Hybrid Overrelaxation in the Gaussian model.Phys. Lett., B288:166–170, 1992.21

[18] Michael Creutz. Monte Carlo algorithms for lattice gauge theory.

Phys.Rev., D36:515, 1987. [19] Frank R. Brown and Thomas J. Woch.Overrelaxed heat bath andMetropolis algorithms for accelerating pure gauge Monte Carlo calcula-tions.

Phys. Rev.

Lett., 58:2394, 1987. [20] Rajan Gupta, Jerry DeLapp, George G. Batrouni, Geoffrey C. Fox,Clive F. Baillie, and John Apostolakis.The phase transition in the2-D XY model.

Phys. Rev.

Lett., 61:1996, 1988. [21] Rajan Gupta, Gregory W. Kilcup, Apoorva Patel, Stephen R. Sharpe,and Philippe de Forcrand.

Comparison of update algorithms for puregauge SU(3). Mod.

Phys. Lett., A3:1367–1378, 1988.

[22] Karsten M. Decker and Philippe de Forcrand. Pure SU(2) lattice gaugetheory on 324 lattices.

In Nicola Cabibbo, Enzo Marinari, Giorgio Parisi,Roberto Petronzio, Luciano Maiani, Guido Martinelli, and Roberto Pet-torino, editors, Lattice ’89, volume B17 of Nuclear Physics (ProceedingsSupplements), pages 567–570, 1990. Proceedings of the 1989 Symposiumon Lattice Field Theory, Capri, Italy, 18–21 September 1989.

[23] John Apostolakis, Clive F. Baillie, and Geoffrey C. Fox. In Heller et al.

[78], page 678. Proceedings of the International Conference on LatticeField Theory, Tallahassee, Florida, USA, 8–12 October 1990.

[24] John Apostolakis, Clive F. Baillie, and Geoffrey C. Fox. Investigationof the two-dimensional O(3) model using the overrelaxation algorithm.Phys.

Rev., D43:2687–2693, 1991. [25] Martin Hasenbusch and Steffen Meyer.

Testing accelerated algorithmsin the lattice CP 3 model. Phys.

Rev., D45:4376–4380, 1992. [26] Ulli Wolff.

Scaling topological charge in the CP 3 spin model. Phys.Lett., B284:94–98, 1992.

[27] UKQCD Collaboration. Towards the continuum limit of SU(2) latticegauge theory.

Phys. Lett., B275:424–428, 1992.22

[28] Stephen L. Adler. An Overrelaxation method for the Monte Carlo eval-uation of the partition function for multiquadratic actions.

Phys. Rev.,D23:2901, 1981.

[29] Stephen L. Adler. Overrelaxation algorithms for lattice field theories.Phys.

Rev., D37:458, 1988. [30] Herbert Neuberger.Adler’s overrelaxation algorithm for Goldstonebosons.

Phys. Rev.

Lett., 59:1877, 1987. [31] Georgios Bathas and Herbert Neuberger.

A possible barrier at z ≈1 forlocal algorithms. Phys.

Rev., D45:3880–3883, 1992. [32] Don Weingarten.

Monte Carlo algorithms for QCD. In Kronfeld andMackenzie [79].

Proceedings of the 1988 Symposium on Lattice FieldTheory, Fermilab. [33] Paul B. Mackenzie.

An improved Hybrid Monte Carlo method. Phys.Lett., B226:369, 1989.

[34] Pietro Rossi. Study of local Hybrid Monte Carlo for SU(3) lattice gaugetheory.

Private Communication, 1992. [35] Steffen Meyer, Brian Pendleton, and Christine Davies.

Phase structureof the SU(2) sector of the standard model. In Z. Horvath, L. Palla, andA.

Patk´os, editors, Frontiers in Nonperturbative Field Theory, page 106.World Scientific, 1989. Eger, Hungary, 18–23 August 1988.

[36] Brian Pendleton. SU(2) Higgs models with dynamical fermions at zeroand finite temperature.

In Kronfeld and Mackenzie [79], page 82. Pro-ceedings of the 1988 Symposium on Lattice Field Theory, Fermilab.

[37] Khalil M. Bitar, Robert G. Edwards, Urs M. Heller, and A. D. Kennedy.On the dynamics of light quarks in QCD. Technical Report FSU–SCRI–92C–161, SCRI, October 1992.

[38] Sourendu Gupta. Dynamical critical properties of the Hybrid MonteCarlo algorithm.

Nucl. Phys., B370:741–761, 1992.23

[39] Sourendu Gupta. Dynamical properties of the Hybrid Monte Carlo algo-rithm.

In Fukugita et al. [77], page 617.

Proceedings of the InternationalSymposium on Lattice Field Theory, Tsukuba, Japan, 5–9 November1991. [40] Alessandro Vaccarino.

Comparison of exact Hybrid and Langevin algo-rithms for QCD with Kogut-Susskind fermions. Phys.

Rev., D44:3254–3257, 1991. [41] John Sloan, Dimitri Kusnezov, and Aurel Bulgac.

A test of the chaoticmolecular dynamics algorithm in the XY model. Technical report, OhioState University, 1992.

YCTP–N1–92. [42] Bernhard Mehlig, A. L. C. Ferreira, and Dieter W. Heermann.

Criticaldynamics of the Hybrid Monte Carlo algorithm. 1992.

[43] A. D. Kennedy and Brian J. Pendleton. Acceptances and autocorrela-tions in Hybrid Monte Carlo.

In Heller et al. [78], pages 118–121.

Talkpresented at “Lattice ’90,” Tallahassee. [44] A. D. Kennedy and Brian J. Pendleton.

Some exact results for HybridMonte Carlo. In preparation, 1992.

[45] Michael Creutz. Global Monte Carlo algorithms for many fermion sys-tems.

Phys. Rev., D38(4):1228–1238, August 1988.

[46] Rajan Gupta, Gregory W. Kilcup, and Stephen R. Sharpe. Tuning theHybrid Monte Carlo algorithm.

Phys. Rev., D38(4):1278–1287, August1988.

[47] J. C. Sexton and D. H. Weingarten.Hamiltonian evolution for theHybrid Monte Carlo algorithm. Nucl.

Phys., B380:665–678, March 1992. [48] Sourendu Gupta, Anders Irb¨ack, Frithjof Karsch, and Bengt Petersson.The acceptance probability in the Hybrid Monte Carlo method.

Phys.Lett., B242:437–443, 1990. [49] Sourendu Gupta.Acceptance and autocorrelation in Hybrid MonteCarlo.

Int. J. Mod.

Phys., C3:43–52, 1992.24

[50] Alan M. Horowitz. A generalized guided Monte Carlo algorithm.

Phys.Lett., B268:247–252, 1991. [51] Martin Grabenstein and Klaus Pinn.Acceptance rates in multigridMonte Carlo.

Phys. Rev., D45:4372, 1992.

[52] Martin Grabenstein and Klaus Pinn. Kinematics of multigrid MonteCarlo.

Technical report, DESY, 1992. DESY–92–094.

[53] Martin Hasenbusch and Steffen Meyer.Multigrid acceleration forasymptotically free theories. Phys.

Rev. Lett., 68:435–438, 1992.

[54] Martin Hasenbusch, Steffen Meyer, and G. Mack. Noncritical multigridMonte Carlo: O(3) nonlinear σ model.

In Heller et al. [78], page 110.Proceedings of the International Conference on Lattice Field Theory,Tallahassee, Florida, USA, 8–12 October 1990.

[55] M. L. Laursen and Jeroen C. Vink. Multigrid updating of the compactU(1) gauge fields in four dimensions.

Technical report, HLRZ, March1992. HLRZ–92–14.

[56] R. G. Edwards, S. J. Ferreira, J. Goodman, and A. D. Sokal. Multi-grid Monte Carlo.

III. Two-dimensional O(4)-symmetric nonlinear sigmamodel.

Nucl. Phys., B380:612–664, 1992.

[57] Thomas Kalkreuter. Improving multigrid and conventional relaxationalgorithms for propagators.

Technical report, DESY, 1992. DESY–92–108.

[58] Arjan Hulsebos, Jan Smit, and Jeroen C. Vink. Multigrid inversion oflattice fermion operators.

Nucl. Phys., B368:379–389, 1992.

[59] R. Ben-Av, A. Brandt, M. Harmatz, E. Katznelson, P.G. Lauwers,S.

Solomon, and K. Wolowesky. Fermion simulations using parallel trans-ported multigrid.

Phys. Lett., B253:185–192, 1991.

[60] P. G. Lauwers, R. Ben-Av, and S. Solomon. Inverting the Dirac ma-trix for SU(2) lattice gauge theory by means of parallel transportedmultigrid.

Nucl. Phys., B374:249–262, 1992.25

[61] P. G. Lauwers and T. Wittlich. Inversion of the fermion matrix in latticeQCD by means of Parallel-Transported Multigrid (PTMG).

Technicalreport, Universit¨at Bonn, September 1992. BONN–92–29.

[62] Vikram Vyas. A multigrid algorithm for calculating fermionic propaga-tors using Migdal–Kadanofftransformation.

Technical report, Univer-sit¨at Wuppertal, 1991. WUB–91–10.

[63] Vikram Vyas. An efficient algorithm for calculating the quark propa-gators using Migdal–Kadanofftransformation.

Technical Report WUB–92–30, Universit¨at Wuppertal, September 1992. [64] Arjan Hulsebos, M. L. Laursen, and Jan Smit.

SU(N) multigrid Landaugauge fixing. Phys.

Lett., B291:431–436, 1992. [65] Arjan J. van der Sijs.Gauge invariant extremization on the lattice.Technical report, Oxford University, 1992.

OUTP–92–16P. [66] Ulli Wolff.Cluster algorithms for nonlinear σ models.Int.

J. Mod.Phys., C3:213–219, 1992.

[67] Hans Gerd Evertz, Matrin Hasenbusch, Mihail Marcu, Klaus Pinn, andSorin Solomon.Cluster algorithms for surfaces.Int. J. Mod.

Phys.,C3:235, 1992. [68] Hans Gerd Evertz, Gideon Lana, and Mihai Marcu.

Cluster algorithmfor vertex models. Technical report, SCRI, August 1992.

[69] Sergio Caracciolo, Robert G. Edwards, Andrea Pelissetto, and Alan D.Sokal. Wolff-type embedding algorithms for general nonlinear σ-models.Technical report, SCRI, 1992.

FSU–SCRI–92–65. [70] B. Bunk.

An efficient algorithm for cluster updates in Z(2) lattice gaugetheories. Technical report, Universit¨at Bielefeld, 1992.

BI–TP–91–22. [71] Mike Flanigan and Pablo Tamayo.

A parallel cluster labeling method forMonte Carlo dynamics. Technical report, Thinking Machines, August1992.

92–0339. [72] Pietro Rossi and G. P. Tecchiolli.

Finding clusters in a parallel environ-ment. 1992.26

[73] Hans Gerd Evertz. Vectorized cluster search.

In Fukugita et al. [77],pages 620–622.

Proceedings of the International Symposium on LatticeField Theory, Tsukuba, Japan, 5–9 November 1991. [74] Michael Creutz.

Microcanonical cluster Monte Carlo. Phys.

Rev. Lett.,69:1002–1005, 1992.

[75] Radford M. Neal. An improved acceptance procedure for the HybridMonte Carlo algorithm.

Technical report, Toronto, July 1992. Print–92–0346.

[76] Gyan Bhanot and Stephen L. Adler. Parallel acceleration algorithm forspin models.

Technical report, IAS Princeton, 1991. IASSNS–HEP–91–90.

[77] M. Fukugita, Y. Iwasaki, M. Okawa, and A. Ukawa, editors. volume B26of Nuclear Physics (Proceedings Supplements), 1992.

Proceedings of theInternational Symposium on Lattice Field Theory, Tsukuba, Japan, 5–9November 1991. [78] Urs M. Heller, A. D. Kennedy, and Sergiu Sanielevici, editors.

volumeB20 of Nuclear Physics (Proceedings Supplements), 1991. Proceedingsof the International Conference on Lattice Field Theory, Tallahassee,Florida, USA, 8–12 October 1990.

[79] Andreas S. Kronfeld and Paul B. Mackenzie, editors.volume B9 ofNuclear Physics (Proceedings Supplements), 1989. Proceedings of the1988 Symposium on Lattice Field Theory, Fermilab.27


출처: arXiv:9212.017원문 보기