---

Multicanonical Multigrid Monte Carlo∗

arXiv:hep-lat/9305016v1 20 May 1993FUB-HEP 9/93May 4, 1993Multicanonical Multigrid Monte Carlo∗Wolfhard Janke1 and Tilman Sauer 21 Institut f¨ur Physik, Johannes Gutenberg-Universit¨at MainzStaudinger Weg 7, 6500 Mainz 1, Germany2 Institut f¨ur Theoretische Physik, Freie Universit¨at BerlinArnimallee 14, 1000 Berlin 33, GermanyAbstractTo further improve the performance of Monte Carlo simulations offirst-order phase transitions we propose to combine the multicanonicalapproach with multigrid techniques. We report tests of this proposi-tion for the d-dimensional Φ4 field theory in two different situations.First, we study quantum tunneling for d = 1 in the continuum limit,and second, we investigate first-order phase transitions for d = 2 in theinfinite volume limit.

Compared with standard multicanonical simu-lations we obtain improvement factors of several resp. of about oneorder of magnitude.∗Work supported in part by Deutsche Forschungsgemeinschaft under grant Kl256.

At first-order phase transitions standard Monte Carlo simulations in thecanonical ensemble exhibit a supercritical slowing down [1]. Here extremelylarge autocorrelation times are caused by strongly suppressed transitions be-tween coexisting phases which, on finite periodic lattices, can only proceedvia mixed phase configurations containing two interfaces.

Since the proba-bility of such configurations is suppressed by a factor exp(−2σLd−1), whereσ is the interface tension and Ld−1 the cross-section of the system, the au-tocorrelation times in the simulation grow exponentially with the size of thesystem, τ ∝exp(2σLd−1).A way to overcome this problem, known as umbrella [2] or multicanonical[3] sampling, is to simulate an auxiliary distribution in which the mixed phaseconfigurations have the same weight as the pure phases and canonical expec-tations are computed by reweighting [4]. Several tests for various models[5] have demonstrated that this method works well in practice and reducessupercritical slowing down to a power-like behavior with τ ∝V α = Ldα,where α ≈1.

While this is clearly an important step forward the remainingslowing down problem is still severe. In most cases it is even worse than forstandard (e.g., Metropolis or heat-bath) Monte Carlo simulations of criticalphenomena [6].For the latter applications several update algorithms have been developedwhich greatly reduce or even completely eliminate the critical slowing downproblem [7].

In addition to overrelaxation and cluster methods, an importantclass of such algorithms are multigrid techniques [8, 9]. Here the generalstrategy is to perform collective updates on different length scales by visitingvarious coarsened grids in a systematic, recursively defined way, generallyknown as V- or W-cycle [10].Because of their conceptual simplicity both the multicanonical reweight-ing approach and the multigrid update techniques are quite generally ap-plicable.

The purpose of this note is to show that the two approaches caneasily be combined and give a much better performance than each compo-nent alone. We report tests of this combination for the Φ4 lattice field theorywith negative mass term in two conceptually different situations.

We firstconsider the quantum mechanical tunneling problem in one dimension andstudy the performance of the new algorithm in the continuum limit. We thendiscuss field driven first-order phase transitions in the two-dimensional caseand investigate the behavior of the multicanonical multigrid algorithm in theinfinite volume limit.1

For Potts models an interesting different approach was proposed onlyrecently in Ref.[11]. Here the idea is to combine a multicanonical demonalgorithm with cluster update methods in a hybrid-like fashion.The basic idea of the multicanonical approach is to sample the mixedphase configurations with the same statistical weight as the configurations ofthe pure phases.

At a field driven first-order phase transition this can alwaysbe achieved by a suitably chosen reweighting factor w−1(m) ≡exp(−f(m)),where m =Pi Φi/V is the average field. In a temperature driven transition,m has simply to be replaced by the average energy.

Starting from an initialguess based on experience or on some analytical approximation, a few iter-ations are usually sufficient to adjust this factor. Once it is fixed, canonicalexpectation values ⟨O⟩can of any observable O can be computed from thebasic reweighting formula⟨O⟩can = ⟨wO⟩⟨w⟩,(1)where ⟨.

. .⟩without subscripts denote expectation values in the multicanon-ical distribution.

To update field values with a Metropolis algorithm in themulticanonical approach, we consider as usual local moves Φi →Φi + ∆Φiand compute the energy difference ∆E. The decision of whether such movesare accepted or not, however, is now based on the value of ∆E + f(m +∆Φi/V ) −f(m).The basic idea of multigrid techniques is to perform updates on dif-ferent length scales.Using the so-called linear interpolation scheme thisamounts, in the equivalent unigrid viewpoint, to proposing moves for blocksof 1, 2d, 4d, .

. .

, V = Ld = 2nd adjacent variables in conjunction, with thesequence of length scales 2k, k = 0, . .

. , n chosen in a specific, recursivelydefined order.

Particular successful sequences are the so-called V-cycle withk = 0, 1, . .

. , n −1, n, n −1, .

. .

, 1, 0 and the W-cycle whose graphical repre-sentation looks like the letter W (for n=3, e.g., this is 0, 1, 2, 3, 2, 3, 2, 1, 2,3, 2, 3, 2, 1, 0, and for large n the W looks more and more like a “fractal”).In a canonical simulation the update at level k thus consists in consideringa common move ∆Φ for all 2kd variables of one block, Φi −→Φi + ∆Φ,i ∈block.The modifications for a multicanonical multigrid simulation are quite triv-ial. Since at level k the proposed move would change the average field by2kd∆Φ/V , the decision of acceptance is now simply to be based on the value2

of ∆E + f(m + 2kd∆Φ/V ) −f(m), with ∆E computed as in the canonicalcase. For didactic reasons we have emphasized here the conceptually sim-pler unigrid viewpoint.

It should be stressed that in the recursive multigridformulation, which can be implemented more efficiently (similar to the fast-Fourier transformation FFT), the multicanonical modification is preciselythe same.We have tested the multicanonical multigrid algorithm for the scalar Φ4lattice field theory in d = 1 and d = 2 dimensions, defined by the partitionfunctionZ =LdYiZ ∞−∞dΦi/Aexp−ǫLdXi=1 12ǫ2(⃗∇Φi)2 −µ22 Φ2i + gΦ4i!,(2)with A =√2πǫ and µ2, g > 0. We always impose periodic boundary condi-tions.

For d = 1 we keep Lǫ = β fixed. Here the model describes the quantumstatistics of a particle tunneling back and forth in a double-well potential incontact with a heat-bath at temperature T = 1/β [12].

At fixed β the limitL →∞corresponds to the continuum limit.For d ≥2 we put ǫ = 1.Here reflection symmetry is spontaneously broken for all µ2 > µ2c(g) > 0 asL →∞, which is now the infinite volume limit. Consequently, if a termhPi Φi is added to the energy, the system exhibits a line of first-order phasetransitions driven by the field h.Even though in the one-dimensional case no spontaneous symmetry break-ing occurs, the numerical difficulties are quite similar.

This is due to thefact that for small quartic coupling g tunneling events are strongly sup-pressed by a factor ∼exp(−const/g), which plays a similar role as the factorexp(−2σLd−1) at a first-order phase transition. The important difference is,of course, that the suppression factor stays roughly constant in the contin-uum limit, while at a first-order phase transition it rapidly decreases in theinfinite volume limit.

Nevertheless, for small values of g analogous slowingdown problems in canonical simulations of the quantum problem are notori-ous, and a number of modified Monte Carlo schemes have been proposed inthe past [13, 14]. None of these techniques, however, is general enough to beeasily adapted to different potential shapes.To evaluate the performance of the multicanonical multigrid algorithm,we have recorded the time series for several observables and studied theirautocorrelation times.

In this note we shall concentrate on the average field3

m = Pi Φi/V , which reflects most directly the tunneling process.In previous investigations [5] emphasis was laid on the exponential au-tocorrelation time τ (0)mof m, i.e., directly on the multicanonical dynamics.While this nicely illustrated the absence of exponential slowing down, it isnot immediately clear how the remaining autocorrelations enter into the errorestimates for canonical expectation values computed according to (1). To beprecise, we are interested in the variance ǫ2 = σ2ˆO = ⟨ˆO2⟩−⟨ˆO⟩2 of the (weaklybiased) estimator for ⟨O⟩can, ˆO = PNm1w(mi)Oi/ PNm1w(mi) ≡wiOi/wi, ifNm (multicanonical) measurements are performed.To facilitate a directcomparison with canonical simulations, we hence define for multicanonicalsimulations an effective autocorrelation time τ effby the standard error for-mula for Nm correlated measurements,ǫ2 = σ2can2τ eff/Nm,(3)where σ2can = ⟨O2i ⟩can −⟨Oi⟩2can is the variance of the canonical distribution ofsingle measurements, which can be computed in a multicanonical simulationby using eq.

(1). The squared error ǫ2 can be estimated either by block-ing or better by jack-knife blocking [15] procedures, or by applying standarderror propagation to the variance of ˆO = wiOi/wi, which involves the (multi-canonical) variances and covariances of wiOi and wi, and the three associatedautocorrelation times τm;m ≡τm, τwm;wm ≡τwm, and τwm;m = τm;wm [16].

Bysymmetry, for O = m this simplifies toǫ2 = ⟨wimi; wimi⟩⟨wi⟩22τwmNm≡σ2muca2τwmNm,(4)where ⟨x; y⟩≡⟨xy⟩−⟨x⟩⟨y⟩and τx;y = 1/2 +Pk⟨x0; yk⟩/⟨x0; y0⟩is theintegrated autocorrelation time of multicanonical measurements. In this wayproperties of the multicanonical distribution (given by σ2muca) are disentangledfrom properties of the update algorithm (given by τwm).

Note that in τ eff=(σ2muca/σ2can)τwm, it is the autocorrelation time of w(m)m that enters and notthat of m, as previously investigated.Let us first discuss our results for the quantum mechanical case in d = 1,where we shall confine ourselves to the case µ2 = 1.0, g = 0.04 and β = 10.This choice of parameters may perhaps be better characterized by the firstfew energy eigenvalues (obtained by a numerical integration of the associatedSchr¨odinger equation), E0 = −0.913371, E1 = −0.892348, E2 = 0.029846,4

and E3 = 0.37813, or the probability ratio Pmin/Pmax = P(0)/Pmax ≈1.9 ×10−3, where P(m) is the probability distribution of the magnetization (or, inthe present interpretation, of the average path). In a canonical simulationit would consequently be about 500 times harder to sample configurationswith m = 0 than configurations contributing to the peaks of P(m).Toset up the multicanonical reweighting factor we started from a variationalapproximation [17] that works quite well for not too large β-values and isknown to provide locally a lower bound on P(m) [18].

Alternatively, as thedistribution depends only weakly on L, one could also use the distribution forsmall L (which is quite easy to generate by standard techniques) as input forthe other simulations. Since the m-values vary continuously, we introducedbins of size ∆m = 0.02 to store the weight factor w(m).

A single short runwas usually sufficient to improve the initial guess such that the multicanonicaldistribution P ′(m) had the desired flat shape, P ′(m) ≈const, between thetwo peaks.In this way we performed multicanonical Metropolis and multigrid simu-lations for L = 4, 8, 16, 32, 64, 128, and 256, and using the multigrid updatealso for L = 512. In the multigrid case we investigated both the V-cycle(using npre = npost = 1 pre- resp.

post-sweeps [8, 9, 10]) and the W-cycle(using npre = npost = 1 as well as npre = 1, npost = 0). In the log-log plotof Fig.

1 we show our results for τ effof the multicanonical Metropolis andW-cycle (without post-sweeps) update algorithm, and for comparison alsoprevious results [19] for the canonical counterparts. We see that canonicaland multicanonical simulations exhibit qualitatively the same behavior in thecontinuum limit L →∞.

For both distributions, the Metropolis update leadsto a power-law growth τ ∝Lz with z ≈2, while for the W-cycle update theautocorrelation times stay roughly constant [20]. Here it is the overall scalewhich is reduced in the multicanonical simulation.

For our choice of param-eters we obtain an improvement factor of about 60 (30) for the Metropolis(W-cycle) update. For smaller g, since this factor is essentially given by theinverse of the suppression factor exp(−const/g), we found the multicanoni-cal Metropolis update to be more favorable than the canonical W-cycle forreasonably large L. Eventually, however, there will always be a crossover atsome L. For g = 0.04, if we compare the Metropolis update and the W-cyclein the multicanonical simulation we obtain an improvement factor of about50, 000 for L = 512.In the continuum limit the distribution P(m) and the ratio Pmin/Pmax5

stay roughly constant. This implies that the multicanonical approach canonly give an improvement factor which is independent of L. In this case it isthe update algorithm that plays the dominant role asymptotically for largeL, while the multicanonical reweighting procedure sets the overall scale.

Therelative importance of the two approaches is just reversed at first-order phasetransitions when the infinite volume limit is considered.As a test case we have studied the model (2) with ǫ = 1. For this modelthe line of second-order phase transitions separating the broken and unbro-ken phase in the µ2 −g plane has recently been determined by Toral andChakrabarti [21].

Here we concentrate on the first-order phase transitionbetween the two ordered phases at g = 0.25 and µ2 = 1.30, which is suffi-ciently far away from the critical point at µ2c = 1.265(5) [21] to display thetypical behavior already on quite small lattices. A sensitive measure of thestrength of the transition is the interface tension σoo between the + and −phase, which turns out [16] to be σoo = 0.03459(49) (for comparison, aboutthe same value is found for the order-disorder interface tension in the two-dimensional q-state Potts model with q = 9, where σod = 0.03355 .

. .

[22]).We performed multicanonical simulations using the Metropolis update andthe W-cycle without post-sweeps for lattices of size V = L2 with L = 8, 16and 32. With the multigrid algorithm, due to the improved performance, wewere also able to study lattices of size L = 64.

Here each time series containsa total of 106 measurements taken every neth sweep, after discarding 104×nesweeps for thermalization. The number of sweeps between measurements,ne, was adjusted in such a way that in each simulation the measurements ofwm had an autocorrelation time of maximal 50, i.e., the length of each timeseries is at least 20, 000 τwm.A few of our results are collected in Table 1, where we give the integratedand exponential autocorrelation times of m and w(m)m as well as τ effac-cording to eq.

(3) for both update algorithms. We see that integrated andexponential autocorrelation times for m agree well with each other, show-ing that the corresponding autocorrelation function can be approximated bya single exponential.

For wm we obtain values for τ (0) that are consistentwith those for m within error bars. The integrated autocorrelation times,however, are significantly lower, implying that the autocorrelation functionis composed of many different modes.

We also observe that the differencebetween τwm and τ effcan be quite appreciable. From L = 8 to L = 64 the ra-tio τ eff/τwm = σ2muca/σ2can varies from about 1.9 to 4.6, reflecting the varying6

probability distribution shapes with increasing L.By fitting τ effto a power law, τ eff∝Lz, we obtain for both updatealgorithms an exponent of z ≈2.3, i.e., in this case it is thus the multigridupdate that reduces the overall scale. The autocorrelation times of the W-cycle are reduced by a roughly constant factor of about 20 as compared withthe Metropolis algorithm.

Of course, for a fair comparison we should alsotake into account that a W-cycle requires more elementary operations thana Metropolis sweep [8]. Such a work estimate, however, depends on manydetails of the implementation and it is hence difficult to give generally validfigures.

With our implementations on a CRAY Y-MP we obtained a realtime improvement factor of about 10.Table 1: Autocorrelation times for the multicanonical simulation using thestandard Metropolis (M) or multigrid W-cycle (W) update algorithm.L = 8L = 16L = 32L = 64τ (0)mM212(12)668(23)3120(200)−W11.30(32)37.2(2.0)148(11)746(62)τmM204.4(4.0)690(11)2984(63)−W10.88(12)34.69(76)150.0(4.0)758(37)τ (0)wmM209(12)655(31)2880(190)−W11.34(33)36.9(2.0)146(13)600(120)τwmM171.1(3.4)509.8(8.9)1840(40)−W9.82(11)27.58(59)96.6(2.4)374(23)τ effM322.7(6.1)1258(21)6050(120)−W18.51(20)67.4(1.3)321.9(7.6)1724(86)AcknowledgmentsW.J. thanks the Deutsche Forschungsgemeinschaft for a Heisenberg fellow-ship.7

References[1] For recent reviews, see, e.g., Dynamics of First Order Phase Transitions,edited by H.J. Herrmann, W. Janke, and F. Karsch (World Scientific,Singapore, 1992).

[2] G.M. Torrie and J.P. Valleau, Chem.

Phys. Lett.

28 (1974) 578; J. Comp.Phys. 23 (1977) 187; I.S.

Graham and J.P. Valleau, J. Phys. Chem.

94(1990) 7894; J.P. Valleau, J. Comp. Phys.

96 (1991) 193. [3] B.A.

Berg and T. Neuhaus, Phys. Lett.

B267 (1991) 249. For a review,see B.A.

Berg, in Ref. [1], p.311 [reprinted in Int.

J. Mod. Phys.

C3(1992) 1083]. [4] W. Janke in Ref.

[1], p.365 [reprinted in Int. J. Mod.

Phys. C3 (1992)1137]; and preprint HLRZ 57/92, J¨ulich (1992).

[5] B.A. Berg and T. Neuhaus, Phys.

Rev. Lett.

68 (1992) 9; W. Janke,B.A. Berg, and M. Katoot, Nucl.

Phys. B382 (1992) 649; B.A.

Berg,U. Hansmann, and T. Neuhaus, Phys.

Rev. B47 (1993) 497; Z. Phys.B90 (1993) 229; A. Billoire, B.A.

Berg, and T. Neuhaus, preprint SPhT-92/120, Saclay (1992); B. Grossmann and M.L. Laursen, in Ref.

[1], p.375[reprinted in Int. J. Mod.

Phys. C3 (1992) 1147]; and preprint HLRZ7/93, J¨ulich (1993); B. Grossmann, M.L.

Laursen, T. Trappenberg, andU.-J. Wiese, Phys.

Lett. B293 (1992) 175.

[6] K. Binder, in Monte Carlo Methods in Statistical Physics, edited by K.Binder (Springer, New York, 1979), p.1.See also the articles in Finite-Size Scaling and Numerical Simulations ofStatistical Systems, edited by V. Privman (World Scientific, Singapore,1990). [7] R.H. Swendsen, J.-S. Wang, and A.M. Ferrenberg, in The Monte CarloMethod in Condensed Matter Physics, edited by K. Binder (Springer,Berlin, 1991); C.F.

Baillie, Int. J. Mod.

Phys. C1 (1990) 91; A.D. Sokal,Monte Carlo Methods in Statistical Mechanics: Foundations and NewAlgorithms, Cours de Troisi`eme Cycle de la Physique en Suisse Ro-mande, Lausanne, 1989.8

[8] J. Goodman and A. D. Sokal, Phys. Rev.

Lett. 56, 1015 (1986); Phys.Rev.

D40, 2035 (1989). [9] D. Kandel, E. Domany, D. Ron, A. Brandt, and E. Loh, Jr., Phys.

Rev.Lett. 60, 1591 (1988); D. Kandel, E. Domany, and A. Brandt, Phys.Rev.

B40, 330 (1989). [10] W. Hackbusch, Multi-Grid Methods and Applications (Springer, Berlin,1985); S.F.

McCormick (ed. ), Multigrid Methods.

Theory, Applications,and Supercomputing (Dekker, New York, 1988). [11] K. Rummukainen, Nucl.

Phys. B390 (1993) 621.

[12] M. Creutz and B. Freedman, Ann. Phys.

132 (1981) 427. For reviewssee, e.g., B.J.

Berne and D. Thirumalai, Ann. Rev.

Phys. Chem.

37(1986) 401; N. Makri, Comp. Phys.

Comm. 63 (1991) 389.

[13] C. Alexandrou, Ph.D. Thesis, MIT, Cambridge (1985); C. Alexandrouand J.W. Negele, Phys.

Rev. C37 (1988) 1513.

[14] E.V. Shuryak and O.V.

Zhirov, Nucl. Phys.

B242 (1984) 393; E.V.Shuryak, Usp. Fiz.

Nauk. 143 (1984) 309 [Sov.

Phys. Usp.

27 (1984)448]. [15] R.G.

Miller, Biometrika 61 (1974) 1; B. Efron, The Jackknife, the Boot-strap and other Resampling Plans (SIAM, Philadelphia, PA, 1982). [16] W. Janke and T. Sauer, in preparation.

[17] R. Giachetti and V. Tognetti, Phys. Rev.

Lett. 55 (1985) 912; R.P.Feynman and H. Kleinert, Phys.

Rev. A34 (1986) 5080; W. Janke andH.

Kleinert, Chem. Phys.

Lett. 137 (1987) 162.

For a comprehensivereview, see H. Kleinert, Path Integrals in Quantum Mechanics, Statisticsand Polymer Physics (World Scientific, Singapore, 1990). [18] W. Janke, in Path Integrals from meV to MeV , proceedings, Bangkok,1989, ed.

V. Sa-yakanit et al. (World Scientific, Singapore, 1989).

[19] W. Janke and T. Sauer, Chem. Phys.

Lett. 201 (1993) 499; and preprintHLRZ 108/92, to be published in Path Integrals from meV to MeV ,proceedings, Tutzing, 1992.9

[20] Using the V-cycle the exponent z can only be reduced to z ≈1 forboth distributions. Also, the extra work required for the W-cycle withnpre = npost = 1 is not completely balanced by reduced autocorrelations.W.

Janke and T. Sauer, to be published. [21] R. Toral and A. Chakrabarti, Phys.

Rev. B42 (1990) 2445; see also A.Milchev, D.W. Heermann, and K. Binder, J. Stat.

Phys. 44 (1986) 749.

[22] C. Borgs and W. Janke, J. Phys. I (France) 2 (1992) 2011.10

Figure HeadingFig. 1: Effective autocorrelation times τ efffor the model (2) in d = 1 asa function of lattice size L with Lǫ = β = 10, µ2 = 1, g = 0.04 fordifferent Monte Carlo algorithms.

The canonical data are taken fromRef. [19].11

---

출처: arXiv:9305.016 • 원문 보기