Geometric allocation approach to accelerating directed worm algorithm

The worm algorithm is a versatile technique in the Markov chain Monte Carlo method for both classical and quantum systems. The algorithm substantially alleviates critical slowing down and reduces the dynamic critical exponents of various classical systems. It is crucial to improve the algorithm and push the boundary of the Monte Carlo method for physical systems. We here propose a directed worm algorithm that significantly improves computational efficiency. We use the geometric allocation approach to optimize the worm scattering process: worm backscattering is averted, and forward scattering is favored. Our approach successfully enhances the diffusivity of the worm head (kink), which is evident in the probability distribution of the relative position of the two kinks. Performance improvement is demonstrated for the Ising model at the critical temperature by measurement of exponential autocorrelation times and asymptotic variances. The present worm update is approximately 25 times as efficient as the conventional worm update for the simple cubic lattice model. Surprisingly, our algorithm is even more efficient than the Wolff cluster algorithm, which is one of the best update algorithms. We estimate the dynamic critical exponent of the simple cubic lattice Ising model to be $z\approx 0.27$ in the worm update. The worm and the Wolff algorithms produce different exponents of the integrated autocorrelation time of the magnetic susceptibility estimator but the same exponent of the asymptotic variance. We also discuss how to quantify the computational efficiency of the Markov chain Monte Carlo method. Our approach can be applied to a wide range of physical systems, such as the ${\left|\varphi \right|}^4$ model, the Potts model, the $\mathrm{O}\left(n\right)$ loop model, and lattice QCD.

Figure 1

Example of a configuration containing the present worm in the square lattice Ising model. The solid lines show the activated bonds, and the broken lines show the deactivated bonds. The solid circles indicate the worm head ($h$) and the worm tail ($t$), both of which break the loop constraint of the activated bonds.

Figure 2

Example of the worm scattering process for the square lattice case. The solid circle in each graph shows the worm head, and the arrow shows the moving direction of the head. When coming to a vertex (state $a$), the worm head scatters to another (or possibly the same) bond (states $b$, $c$, $d$, and $e$) with a certain probability.

Figure 3

Geometric allocation for the square lattice Ising model. There are two cases: (a) and (b). The solid (broken) lines show the activated (deactivated) bonds, and the solid circles show the worm head. The weight, or the measure, of each state is denoted by ${\pi }_i$ ($i=1,2,3,4$) apart from the normalization factor of the target distribution, and the allocated raw flow from $i$ to $j$ is denoted by $v_{ij}$. The detailed balance condition is satisfied in both cases: $v_{ij}=v_{ji}$. In the case of (a), we set $v_{12}={\pi }_4$, $v_{13}=v_{14}=\frac{1}{2}({\pi }_1-{\pi }_4)$, and $v_{34}=\frac{1}{2}(3{\pi }_4-{\pi }_1)$. This rejection-free allocation can be performed if $3{\pi }_4>{\pi }_1\iff T<2/ln2$. In the case of (b), we set $v_{12}=\frac{1}{2}({\pi }_1+{\pi }_4)$, $v_{13}=v_{23}=\frac{1}{2}({\pi }_1-{\pi }_4)$, and $v_{34}={\pi }_4$, which is possible at any temperature. The scale of the area ${\pi }_i$ is arbitrary; only the ratio ${\pi }_4/{\pi }_1$ matters.

Figure 4

Geometric allocation for the simple cubic lattice Ising model. There are four cases: (a), (b), (c), and (d). The six possible states are indexed such that (1, 2), (3, 4), and (5, 6) are the pairs of the states from and to which the worm head forward scatters like the square lattice case. The detailed balance condition is satisfied in all the cases: $v_{ij}=v_{ji}$. In the case of (a), we set $v_{12}={\pi }_6$, $v_{13}=v_{14}=v_{15}=v_{16}=\frac{1}{4}({\pi }_1-{\pi }_6)$, and $v_{34}=v_{56}=\frac{1}{4}(5{\pi }_6-{\pi }_1)$; in (b), $v_{12}=\frac{1}{2}({\pi }_1+{\pi }_6)$, $v_{13}=v_{23}=\frac{1}{2}({\pi }_1-{\pi }_6)$, and $v_{34}=v_{56}={\pi }_6$; in (c), $v_{12}=v_{34}=v_{56}={\pi }_6$ and $v_{13}=v_{15}=v_{35}=\frac{1}{2}({\pi }_1-{\pi }_6)$; in (d), $v_{12}=v_{34}=\frac{1}{4}(3{\pi }_1+{\pi }_6)$, $v_{15}=v_{25}=v_{35}=v_{45}=\frac{1}{4}({\pi }_1-{\pi }_6)$, and $v_{56}={\pi }_6$. This rejection-free allocation in (a) can be performed if $5{\pi }_6>{\pi }_1\iff T<2/ln(3/2)$. The allocations in (b), (c), and (d) are possible at any temperature.

Figure 5

Probability distributions of the difference between the coordinates of the two kinks (head and tail) for the $L=128$ square lattice model at the critical temperature in (a) the classical and (b) the present worm updates. The coordinates were measured 256 worm shifting or scattering steps after the insertion. The contours show the coordinates at which $P(x,y)=0.5,1.0,1.5,2.0$, and $2.5\times {10}^{-4}$. Because the two kinks are removed when meeting each other, the distribution is lowered near the center, which is more significant in the classical worm update. The removed worms are not shown here but counted in the normalization.

Figure 6

Tails of the probability distributions of the distance between the two kinks in the classical (open) and the present (solid) worm updates 64 (triangles), 128 (squares), and 256 (circles) local worm steps after the insertion, measured in the $L=128$ square lattice Ising model at the critical temperature. The distribution was measured at $r=\left|\mathbf{r}\right|$, where $\mathbf{r}=(x,y)$ and $\left|x\right|=\left|y\right|$. The tails are fitted to Gaussian distributions: $P\left(r\right)\propto e^{-r^2/2{\sigma }^2}$, where ${\sigma }^2$ is a parameter (variance). The inset shows the linear scaling of the estimated variance in the classical (circles) and the present (squares) worm updates as a function of the number of worm steps ($s$). The variance in the present worm update is approximately six times as large as in the classical worm update.

Figure 7

(a) The integrated autocorrelation time, (b) the variance, and (c) the asymptotic variance of the total energy estimator as a function of the system length of the simple cubic lattice Ising model at the critical temperature. The exponents of ${\tau }_{\mathrm{int},\widehat{E}}$ are estimated to be 0.28, 0.31, and 0.27, and those of $v_{\mathrm{asymp},\widehat{E}}$ are to be $-2.44$, $-2.35$, and $-2.46$ in the Wolff cluster (triangles), the classical worm (circles), and the present worm (squares) updates, respectively; the exponent of $v_{\widehat{E}}$ is estimated to be $-2.75$ in all the updates. The inset of panel (c) shows the ratios of the asymptotic variance in the classical worm (diamonds) and the Wolff cluster (pentagons) updates to the one in the present worm update. They are approximately 27 and 2.2 for large system sizes, respectively.

Figure 8

(a) The integrated autocorrelation time, (b) the variance, and (c) the asymptotic variance of the magnetic susceptibility estimator as a function of the system length of the simple cubic lattice Ising model at the critical temperature in the Wolff (triangles), the classical (circles), and the present worm (squares) algorithms. In the Wolff algorithm, we test two estimators using the spins (Wolff spin) and the cluster size (Wolff cluster) (see the main text for the detail of the estimators). The exponents of ${\tau }_{\mathrm{int},\widehat{\chi }}$ are estimated to be 0.150(9), $-0.50\left(1\right)$, $-0.731\left(7\right)$, and $-0.679\left(4\right)$; those of $v_{\widehat{\chi }}$ are to be 0.01, 0.58, 0.92, and 0.85; and those of $v_{\mathrm{asymp},\widehat{\chi }}$ are to be 0.18, 0.18, 0.22, and 0.18 in the Wolff spin, in the Wolff cluster, in the classical worm, and in the present worm updates, respectively. The inset of panel (c) shows the ratios of the asymptotic variance in the classical worm (diamonds) and the Wolff cluster (pentagons) updates to the one in the present worm update, which are approximately 23 and 1.6 for large system sizes, respectively.

Figure 9

Autocorrelation functions of (a) the total energy and (b) the magnetic susceptibility estimators in the classical (circles) and the present (squares) worm updates for the $L=4$ (open) and 8 (solid) simple cubic lattice Ising model. The horizontal axis is the rescaled time of the Monte Carlo dynamics in units of $L^3$ worm shifting or scattering steps.

Figure 10

Exponential autocorrelation times of the energy (solid) and the magnetic susceptibility (open) estimators as a function of the system length in the classical (circles) and the present (squares) worm updates for the simple cubic lattice Ising model. The inset shows the ratio of the autocorrelation time of the energy estimator in the classical algorithm to the one in the present algorithm.