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

The probability model

State space: Ω_{Λ,N} = Λ^{N} × S_{N}, where Λ = [0, L]^{3} with periodic boundary conditions. Point positions: X = (x_{1}, …, x_{N}) for x_{1}, …, x_{N} ∈ Λ.

The Hamiltonian arises from consideration of the Bose gas [Betz and Ueltschi]:

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

X’s (for now) are on the cubic unit lattice. Gibbs distribution: P (π) = e^{−H(π)}/Z. Not important for today’s talk: α_{ℓ}’s are cycle weights; r_{ℓ}(π) counts the number of ℓ-cycles in π.

High T: identity only. Low T: uniform distribution. Intermediate T: individual jump lengths remain short, but below a T_{c}, long cycles form. Quantification of ΔT_{c} as function of interaction strength is discussed in my dissertation (March 2010).

J. Kerl (Arizona)  A worm algorithm for random spatial permutations  February 22, 2010  2 / 11

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.

Metropolis step (analogue of single spin-flips for the Ising model): swap permutation arrows which end at nearest-neighbor lattice sites. This either splits a common cycle, or merges disjoint cycles:

As usual, the proposed change π → π′ is accepted with probability min{1, e^{−ΔH}}.

J. Kerl (Arizona)  A worm algorithm for random spatial permutations  February 22, 2010  3 / 11

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

Figure part 1: a long cycle on the torus almost meets itself in the x direction.

Part 2: after a swap-only step (above), one cycle winds by +1, and the other by −1. Metropolis steps create winding cycles only in opposite-direction pairs; total W_{x}(π) is still zero.

Part 3: if we reverse one cycle (zero-energy move), W_{x}(π) is now 2. This is the swap-and-reverse algorithm. This permits winding numbers of even parity in each of the three axes.

One idea to get all winding numbers: band updates. Compose π with a winding L-cycle τ. We obtain π′ = τπ, with winding number shifted by ±1 along a specified axis.

Problem: acceptance rate is approximately e^{−L}; too low.

Worm algorithm

Another idea: Adapt PIMC worm methods (Ceperley and Pollock 1986, and many others since), which break and re-join Brownian bridges. Here, open and re-close cycles. An open cycle will be able to wander around the torus, tunneling through the winding-number energy barrier which closed permutations have.

Use permutations on N + 1 points: the (N + 1)st is the wormhole point, w. Closed permutations have π(w) = w; open permutations have π(w) ≠ w.

Goal: invent 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: 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}.

J. Kerl (Arizona)  A worm algorithm for random spatial permutations  February 22, 2010  5 / 11

Extended random-cycle model

Extended lattice: Λ′ = Λ ∪ {w}. Extended energy:

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

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 holds.

• 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.

• This worm method will work for any Hamiltonian on the model of random spatial permutations — not only for the one given above.

Metropolis steps/sweeps for the worm algorithm

Irreducibility, aperiodicity, and detailed balance have been proved for Metropolis steps, using an H′ with γ term as above. Metropolis sweep: open, zero or more head/tail swaps, close. RVs 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)

Stopping time

Good news: 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.

Problem: 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. Interesting L (40-80 or so) are out of reach for the worm algorithm.

Other ideas

Other ideas for addressing the winding-number problem include the following:

• Be content with even winding numbers. Be sure to quantify the shift in critical temperature by means of multiple random variables, some of which rely on winding phenomena and some of which do not. This approach is taken in my dissertation.

• 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.

• Go to Athens … For more information, please visit http://math.arizona.edu/~kerl.

Thank you for attending!

Extra slide: random spatial permutations derived from Bose gas

Start with the Bose-gas Hamiltonian H(X) = − ∑_{i=1}^{N} ∇_{i}^{2} + ∑_{1≤i,j≤N} U(x_{i} − x_{j}), where U is a hard-core potential of radius a.

Write the partition function Z = Tr_{L_{sym}^{2}}(e^{−βH}) = Tr_{L^{2}} (P_{+}e^{−βH}) = Tr_{L^{2}} (e^{−βH}P_{+}) where P_{+} f(x_{1}, …, x_{N}) := \frac{1}{N!} ∑_{π∈S_{n}} M_{π}f(x_{1}, …, x_{N}) and M_{π}(fx_{1}, …, x_{N}) := f(x_{π(1)}, …, x_{π(N)}). That is, Tr_{L_{sym}^{2}}(e^{−βH}) = \frac{1}{N!} ∑_{π∈S_{N}} Tr_{L^{2}} (e^{−βH}M_{π}) .

Then:

• Interpret e^{−βH}M_{π} as an expectation over Brownian motions.

• Write e^{−βH}M_{π} as an integral operator, and find the kernel.

• Compute Tr (e^{−βH}M_{π}) in terms of Brownian bridges.

• Sum over π ∈ S_{N} to obtain Z = Tr_{L_{sym}^{2}}(e^{−βH}). Now Z is a sum over permutations π of e^{−H_{P} (X,π)}.

• Decouple the non-interacting from the interacting terms in the permutation Hamiltonian, so that we may write e^{−H_{P}^{(0)}(X,π)−H_{P}^{(1)}(X,π)}.

• Use a cluster expansion to drop all but 2-jump interactions, and find the logarithm of e^{−H_{P} (X,π)}. We recognize the random-cycle model with an explicit 2-jump interaction V.

J. Kerl (Arizona)  A worm algorithm for random spatial permutations  February 22, 2010  11 / 11