A worm algorithm for random spatial permutations (2010 Center for Simulational Physics Workshop, U. of Georgia)

John Kerl Department of Mathematics, University of Arizona

February 22, 2010

Abstract

Models of random spatial permutations arise in the study of Bose-Einstein condensation. Namely, permutations of sites occur with probabilities depending on lengths of permutation jumps, as well as on interactions between jumps. Below a critical temperature, one observes the onset of long permutation cycles in spite of short individual jump lengths. We have devised several Markov chain Monte Carlo algorithms for sampling from this probability distribution. In this note, we present one particularly promising technique: a worm algorithm. It admits an elegant correctness theory, asymptotically in the number of completed Metropolis sweeps. However, it suffers from a stopping-time problem: the CPU time needed to complete a sweep is strongly quadratic in the number of lattice points N.

1 The probability model

The model of random spatial permutations arises from consideration of the Bose gas [Sütő1, Sütő2, GRU, U07, BU07, BU08, SU09]. See also [BBHLV, SU09] for the contentious history of studies of the shift in critical temperature as a function of interaction strength. The state space is Ω_{Λ,N} = Λ^{N} × S_{N}, where Λ = [0, L]^{3} with periodic boundary conditions. Point positions are X = (x_{1}, …, x_{N}) for x_{1}, …, x_{N} ∈ Λ. The Hamiltonian is

H(π) = \frac{T}{4} ∑_{i=1}^{N} ∥x_{i} − x_{π(i)}∥^{2} + ∑_{ℓ=1}^{N} α_{ℓ}r_{ℓ}(π).

(The temperature scale factor T/4, not β/4, is surprising but correct for the Bose-gas derivation of the Hamiltonian.) We obtain a Gibbs distribution:

Y (Λ, X) = ∑_{σ∈S_{N}} e^{−H(X,σ)},  P (π) = P_{Λ,X}(π) = \frac{e^{−H(X,π)}}{Y (Λ, X)} .  (1.1)

Note that higher-energy states are lower-probability states. The α_{ℓ}’s are called cycle weights; r_{ℓ}(π) counts the number of ℓ-cycles in π. For my dissertation [Kerl], point positions X are held fixed on the cubic unit lattice, and the cycle weights are in the form of [Ewens]. That is, α_{ℓ} ≡ α are constant in ℓ. These weights are taken to be small: their effect is to perturb the critical temperature slightly; qualitative behavior of the model is not affected. Thus, the cycle-weight terms are not particularly important for this note; the reader may safely ignore them.

Figure 1: A spatial permutation on N = 26 points. There are 11 one-cycles, three two-cycles, one four-cycle, and one five-cycle. We say r_{1}(π) = 11, r_{2}(π) = 3, r_{4}(π) = 1, r_{5}(π) = 1, and r_{ℓ}(π) = 0 for all other ℓ.

Typical permutations are as follows. As T → ∞, any non-zero distance-dependent term carries a high energy penalty; only the identity permutation appears. As T → 0, distance dependence vanishes and we obtain the uniform distribution. At intermediate T, individual jump lengths remain short — but below a T_{c} ≈ 6.86, long cycles form. Quantification of ΔT_{c} as function of interaction strength is the main result of [Kerl].

2 Metropolis sampling

The expectation of a random variable S(π) is 𝔼[S] = ∑_{π∈S_{N}} P (π)S(π). The number of permutations, N!, grows intractably in N. The expectation is instead estimated by summing over some number M (10^{5} or 10^{6}) typical permutations. The sample mean ⟨X⟩_{M}is now a random variable with its own variance. The usual technical issues of Markov chain Monte Carlo (MCMC) methods are known and handled in my simulations and dissertation: thermalization time, proofs of detailed balance, autocorrelation, batched means, quantification of variance of sample means, and finite-size-scaling analysis of finite-lattice computational results. (See [Berg, LB, PV] for background on such methods as applied to models other than the model of random spatial permutations.)

The Metropolis step (the analogue of single spin-flips for the Ising model) is to swap permutation arrows which end at nearest-neighbor lattice sites. This either splits a common cycle, or merges disjoint cycles, as shown in figure 2. As usual [Berg, LB], the proposed change π → π′ is accepted with probability min{1, e^{−ΔH}}.

Figure 2: Swaps merge disjoint cycles and split single cycles. The left-hand permutation can be reached from the right-hand permutation via a swap, and vice versa.

3 Winding numbers: SO, SAR, and band-update algorithms

Part 1 of figure 3 shows a long cycle on the torus which almost meets itself in the x direction. In part 2, after a swap-only step (as described above), one cycle winds by +1, and the other by −1. Metropolis steps create winding cycles only in opposite-direction pairs; total W_{x}(π) (the number of wraps around the L-torus in the x axis) is still zero. In part 3, if we reverse one cycle (which is a zero-energy move), W_{x}(π) is now 2. Combining swap sweeps and cycle-reversal sweeps constitutes the swap-and-reverse algorithm. It permits Figure 3: Conservation of winding number in the swap-only algorithm, and a partial solution provided by the swap-and-reverse algorithm.

winding numbers of even parity in each of the three axes; it is used for all computational results presented in [Kerl].

One idea to get all winding numbers uses band updates. Namely, we compose π with a winding L-cycle τ. We obtain π′ = τπ, with winding number shifted by ±1 along a specified axis. However, the acceptance rate is found to be approximately e^{−L}, i.e. proposed band updates are effectively never accepted in an MCMC simulation.

4 Worm algorithm

Another approach to the winding-number problem is to adapt path-integral Monte Carlo (PIMC) worm methods. These were initially developed by Ceperley and Pollock, and have been used by many others since (e.g. [PC87, BPS06, PST98]). Such methods break and re-join Brownian bridges. Here, we want to open and re-close permutation cycles. An open cycle will be able to wander around the torus, tunneling through the winding-number energy barrier which closed permutations have.

The essential idea (figure 4) is to use permutations on

Figure 4: Open cycles as permutations on N + 1 points.

N + 1 points: the (N + 1)st is the wormhole point, w. Closed permutations have π(w) = w; open permutations have π(w) ≠ w. We invent an extended random cycle model with an energy function H′, Gibbs distribution P′, and Metropolis algorithm for S_{N+1} such that the marginal distribution on S_{N+1}, conditioned on closed permutations, matches the RCM Gibbs distribution. Then, random variables will be sampled only at closed permutations.

Theorem [Kerl]: Let S_{N} ↪ S_{N+1} by π(w) = w. If H(π) = H′(π) for all π ∈ S_{N}, then P′(π | π ∈ S_{N}) = P (π) for all π ∈ S_{N}.

The extended lattice is Λ′ = Λ ∪ {w}. The extended energy is

H′(π) = \frac{T}{4} ∑_{π(x_{i})≠w^{i=1}}^{N} ∥x_{i} − x_{π(i)}∥^{2} + ∑_{ℓ=2}^{N} α_{ℓ}r_{ℓ}(π) + γ1_{S_{N+1}\S_{N}} (π).

We define partition function Z′ and Gibbs distribution P′ as usual. As long as the energy function for the ERCM and the RCM agree on closed permutations, the marginality condition of the theorem holds. Thus, we are free to define energy terms in H′ for open permutations (the γ factor is just one possibility), as long as they vanish on closed permutations. Furthermore, this worm method will work for any Hamiltonian on random permutations — not only for the cycle-weight Hamiltonian given above.

5 Metropolis steps/sweeps for the worm algorithm; stopping time

Irreducibility, aperiodicity, and detailed balance have been proved [Kerl] for worm Metropolis steps, using an H′ with γ term. A Metropolis sweep consists of an open, zero or more head/tail swaps, and a close. Random variables are sampled only at closed permutations.

PSfrag replacements

x π(x)

w x

π(x)

w

Open at x w.p.

∼ 1 ∧ e^{−ΔH}

π^{−1}(w)

π(w)

w π^{−1}(w)

π(w)

w Close w.p. ∼ 1 ∧ e^{−ΔH}

x π(x)

Head swap at x w.p.

∼ 1 ∧ e^{−ΔH} Tail swap at x w.p.

∼ 1 ∧ e^{−ΔH}

x π(x)

π^{−1}(w) w x

π(x) π^{−1}(w)

w x

π(x)

w π(w)

w π(w)

Figure 5: Metropolis moves for the worm algorithm.

The good news is that examination of random-variable plots for L = 10, comparing SAR to worm, shows that similar results are produced — other than, of course, the winding-number histogram itself. The problem is that the the open worm tips wander around randomly within the L box, and fail to reconnect as L increases. Specifically, histograms show that the distribution of the wormspan ∥π(w) − π^{−1}(w)∥ peaks around L/2. SAR and worm CPU times are both ∼ aN + bN^{2}. SAR’s b is tiny; worm’s b is not (figure 6). Interesting L (40-80 or so) are out of reach for the worm algorithm.

Figure 6: Scalability of SAR and worm algorithms. CPU times for 10^{4} SAR sweeps and 10^{3} worm sweeps are shown as a function of N = L^{3} for L = 5 to 12. SAR time is nearly linear in N; worm time is strongly quadratic in N. Interesting L (40-80) are unattainable.

6 Other ideas

Other ideas for addressing the winding-number problem include the following. (*) Band updates (above) have too-low acceptance rate. (*) Temporarily pinch the torus geometry somehow in the SAR algorithm, such that the distance penalty for wrapping around the torus is decreased. (*) Temporarily reduce and restore the temperature T in the SAR algorithm — this is an annealing method. This approach brings with it a performance problem: re-thermalization would need to be performed after each annealing step. (*) Modify the worm algorithm to direct the worm somehow. At the time the worm is opened, add a distance weight of ±L in the x, y, or z direction which will be removed by a wrap around the torus, increasing or decreasing that winding-number component by 1. Our attempts to do this have not satisfied detailed balance. (*) Review the PIMC literature again and seek other inspiration.

References

[BBHLV] Baym, G., Blaizot, J.-P., Holzmann, M., Laloë, F., and Vautherin, D. Bose-Einstein transition in a dilute interacting gas. arXiv:cond-mat/0107129v2. Eur. Phys. J. B 24, 107-124 (2001).

[Berg] Berg, B. Markov Chain Monte Carlo Simulations and Their Statistical Analysis. World Scientific Publishing (2004).

[BPS06] Boninsegni, M., Prokof’ev, N.V., and Svistunov, B.V. Worm algorithm and diagrammatic Monte Carlo: A new approach to continuous-space path integral Monte Carlo simulations. Physical Review E 74, 036701 (2006).

[BU07] Betz, V. and Ueltschi, D. Spatial random permutations and infinite cycles. arXiv:0711.1188. Commun. Math. Phys. 285, 469-501 (2009).

[BU08] Betz, V. and Ueltschi, D. Spatial random permutations with small cycle weights. arXiv:0812.0569v1. Probabl. Th. Rel. Fields (2010).

[Ewens] Ewens, W.J. The sampling theory of selectively neutral alleles. Theor. Popul. Biol. 3, 87-112 (1972).

[GRU] Gandolfo, D., Ruiz, J., and Ueltschi, D. On a model of random cycles. arXiv:cond-mat/0703315. Statist. Phys. 129, 663-676 (2007).

[Kerl] Kerl, J. Critical behavior for the model of random spatial permutations. Doctoral dissertation, University of Arizona, May 2010 (in progress).

[LB] Landau, D.P. and Binder, K. A Guide to Monte Carlo Simulations in Statistical Physics (2nd ed.). Cambridge University Press (2005).

[PC87] Pollock, E.L. and Ceperley, D.M. Path-integral computation of superfluid densities. Physical Review B, vol. 36, no. 16 (1987).

[PST98] Prokof’ev, N.V., Svistunov, B.V., and Tupitsyn, I.S. Exact, complete, and universal continuous-time worldline Monte Carlo approach to the statistics of discrete quantum systems. Journal of Experimental and Theoretical Physics, vol. 87, no. 2 (1998).

[PV] Pelissetto, A. and Vicari, E. Critical Phenomena and Renormalization-Group Theory. Physics Reports, 368 (6), p.549-727, Oct 2002. arxiv:cond-mat/0012164.

[Sütő1] Sütő, A. Percolation transition in the Bose gas. J. Phys. A: Math. Gen. 26 (1993) 4689-4710.

[Sütő2] Sütő, A. Percolation transition in the Bose gas II. J. Phys. A: Math. Gen. 35 (2002) 6995-7002.

[SU09] Seiringer, R. and Ueltschi, D. Rigorous upper bound on the critical temperature of dilute Bose gases. arXiv.org:0904.0050. Phys. Rev. B 80, 014502 (2009).

[U07] Ueltschi, D. The model of interacting spatial permutations and its relation to the Bose gas. arXiv:0712.2443v3. Mathematical Results in Quantum Mechanics, pp. 225-272, World Scientific (2008).