%\documentclass[12pt]{amsart}
\documentclass[10pt,oneside]{article}
\parindent=0pt \parskip=8pt

\usepackage{verbatim}
\usepackage{fullpage}

\usepackage{amssymb,amsmath,amsthm}
\usepackage[mathscr]{eucal}

\usepackage{makeidx}
%\numberwithin{equation}{subsection}
\numberwithin{equation}{section}
\input{isomath}
\input{mathenv}
\input{syms}
\usepackage{graphicx}
\usepackage{psfrag}
\usepackage{afterpage}
\usepackage{subfigure}
\usepackage{floatflt}

\makeindex

\usepackage{hyperref}
\usepackage{ifpdf}

\ifpdf
\pdfinfo{
/Author (John Kerl)
/Title  (A worm algorithm for random spatial permutations}
}
\fi

%--------------------------------------------------------------------
\begin{document}
\title{A worm algorithm for random spatial permutations \\
	(2010 Center for Simulational Physics Workshop, U. of Georgia)}
\author{John Kerl \\
	Department of Mathematics, University of Arizona}
\date{February 22, 2010}

\maketitle

\begin{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$.
\end{abstract}

%% ----------------------------------------------------------------
\section{The probability model}

The model of random spatial permutations arises from consideration of the Bose
gas \cite{bib:suto1,bib:suto2,bib:gru,bib:qmath,bib:bu07,bib:bu08,bib:su09}.
See also \cite{bib:baymetal,bib:su09} for the contentious history of studies of
the shift in critical temperature as a function of interaction strength.  The
state space is $\Omega_{\Lambda,N} = \Lambda^N \times \SN$, where
$\Lambda=[0,L]^3$ with periodic boundary conditions.  Point positions are
$\vecX = (\vecx_1, \ldots, \vecx_N)$ for $\vecx_1, \ldots, \vecx_N \in
\Lambda$.
The Hamiltonian is
\begin{align*}
	H(\pi) = \frac{T}{4} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
		+ \sum_{\ell=1}^N \alpha_\ell r_\ell(\pi).
\end{align*}
(The temperature scale factor $T/4$, not $\beta/4$, is surprising but correct
for the Bose-gas derivation of the Hamiltonian.)
We obtain a Gibbs distribution:
\begin{align}
\label{eqn:RCM_Gibbs_P}
	Y(\Lambda,\vecX) &= {\sum_{\sigma\in\SN} e^{-H(\vecX,\sigma)}}, &
	P(\pi) = P_{\Lambda,\vecX}(\pi)
	&= \frac{e^{-H(\vecX,\pi)}} {Y(\Lambda,\vecX)}.
\end{align}
Note that higher-energy states are lower-probability states.  The
$\alpha_\ell$'s are called cycle weights; $r_\ell(\pi)$ counts the number of
$\ell$-cycles in $\pi$.  For my dissertation \cite{bib:kerl}, point positions
$\vecX$ are held fixed on the cubic unit lattice, and the cycle weights are in
the form of \cite{bib:ewens}.  That is, $\alpha_\ell \equiv \alpha$ are
constant in $\ell$.  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.

\begin{figure}[!htp]
%\begin{floatingfigure}[r]{3in}
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{figures/pmt1.eps}
\caption[A spatial permutation on $N=26$ points.]
	{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(\pi)=11$,
	$r_2(\pi)=3$, $r_4(\pi)=1$, $r_5(\pi)=1$, and $r_\ell(\pi) = 0$ for all
	other $\ell$.
	\label{fig:pmt1}}
\end{center}
%\end{floatingfigure}
\end{figure}

Typical permutations are as follows.  As $T \to \infty$, any non-zero
distance-dependent term carries a high energy penalty; only the identity
permutation appears.  As $T \to 0$, distance dependence vanishes and we obtain
the uniform distribution.  At intermediate $T$, individual jump lengths remain
short --- but below a $T_c \approx 6.86$, long cycles form.  Quantification of
$\Delta T_c$ as function of interaction strength is the main result of
\cite{bib:kerl}.

%% ----------------------------------------------------------------
\section{Metropolis sampling}

The expectation of a random variable $S(\pi)$ is $\bbE[S] = \sum_{\pi\in\SN}
P(\pi) S(\pi)$.  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 $\pigm X$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
\cite{bib:berg,bib:landaubinder,bib: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
\ref{fig:socycles}.  As usual \cite{bib:berg,bib:landaubinder}, the proposed
change $\pi \to \pi'$ is accepted with probability $\min\{1, e^{-\Delta H}\}$.

\begin{figure}[!htb]
%\begin{floatingfigure}[r]{5in}
\begin{center}
\psfragscanon
\includegraphics[scale=0.5]{figures/gkcycles.eps}
\caption[Swaps merge disjoint cycles and split single cycles.]
	{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.
	\label{fig:socycles}}
\end{center}
%\end{floatingfigure}
\end{figure}

%% ----------------------------------------------------------------
\section{Winding numbers: SO, SAR, and band-update algorithms}

Part 1 of figure \ref{fig:sowinding} 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(\pi)$ (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(\pi)$ is now 2.  Combining swap sweeps and cycle-reversal sweeps
constitutes the swap-and-reverse algorithm.  It permits winding numbers of even
parity in each of the three axes; it is used for all computational results
presented in \cite{bib:kerl}.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{figures/gkwinding.eps}
\caption[Conservation of winding number in the SO algorithm.]
	{Conservation of winding number in the swap-only algorithm, and a
	partial solution provided by the swap-and-reverse algorithm.
	\label{fig:sowinding}}
\end{center}
\end{figure}

One idea to get all winding numbers uses band updates.  Namely, we compose
$\pi$ with a winding $L$-cycle $\tau$.  We obtain $\pi' = \tau\pi$, with
winding number shifted by $\pm 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.

%% ----------------------------------------------------------------
\section{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.
\cite{bib:pc87,bib:bps06,bib: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.
\begin{floatingfigure}[r]{3in}
\begin{center}
\psfragscanon
\includegraphics[scale=0.3]{figures/worm1slide.eps}
\caption{Open cycles as permutations on $N+1$ points.
	\label{fig:worm1}}
\end{center}
\end{floatingfigure}
The essential idea (figure \ref{fig:worm1}) is to use permutations on $N+1$
points:  the $(N+1)$st is the wormhole point, $w$.  Closed permutations have
$\pi(w)=w$; open permutations have $\pi(w) \ne w$.  We invent an
\emph{extended random cycle model} with an energy function $H'$, Gibbs
distribution $P'$, and Metropolis algorithm for $\SNp$ such that the marginal
distribution on $\SNp$, conditioned on closed permutations, matches the RCM
Gibbs distribution.  Then, random variables will be sampled only at closed
permutations.

Theorem \cite{bib:kerl}: Let $\SN \into \SNp$ by $\pi(w)=w$.  If $H(\pi) =
H'(\pi)$ for all $\pi\in\SN$, then $P'(\pi \mid \pi\in\SN) = P(\pi)$ for all
$\pi\in\SN$.

The extended lattice is $\Lambda' = \Lambda \cup \{w\}$.  The extended energy
is
\begin{align*}
	H'(\pi) &= \frac{T}{4} \sum_{\stackrel{i=1}{\pi(\vecxi) \ne w}}^N
		\|\vecxi - \vecxpii\|^2 + \sum_{\ell=2}^N \alpha_\ell r_\ell(\pi)
		+ \gamma 1_{\SNpc}(\pi).
\end{align*}
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 $\gamma$ 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.

%% ----------------------------------------------------------------
\section{Metropolis steps/sweeps for the worm algorithm; stopping time}

Irreducibility, aperiodicity, and detailed balance have been proved
\cite{bib:kerl} for worm Metropolis steps, using an $H'$ with $\gamma$ 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.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.50]{figures/wormmoves2x2.eps}
\caption{Metropolis moves for the worm algorithm.
	\label{fig:wormmoves}}
\end{center}
\end{figure}

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
\plainidx{wormspan} $\|\pi(w) - \pi^{-1}(w)\|$ peaks around $L/2$.  SAR and
worm CPU times are both $\sim$ $aN + bN^2$.  SAR's $b$ is tiny; worm's $b$ is
not (figure \ref{fig:worm_vs_sar_perf}).  Interesting $L$ (40-80 or so) are out
of reach for the worm algorithm.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.65]{figures/sar_perf_small_L.eps}
\includegraphics[scale=0.65]{figures/worm_perf_small_L.eps}
\caption[Scalability of SAR and worm algorithms.]
	{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.
	\label{fig:worm_vs_sar_perf}}
\end{center}
\end{figure}

%% ----------------------------------------------------------------
\section{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 $\pm 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.

%% ================================================================
\begin{thebibliography}{}

\bibitem[BBHLV]{bib:baymetal}
Baym,~G., Blaizot,~J.-P., Holzmann,~M., Lalo\"e,~F., and Vautherin,~D.
\newblock {\em Bose-Einstein transition in a dilute interacting gas}.
\newblock \texttt{arXiv:cond-mat/0107129v2}.
\newblock Eur. Phys. J. B 24, 107-124 (2001).

\bibitem[Berg]{bib:berg}
Berg,~B.
\newblock {\em Markov Chain Monte Carlo Simulations and Their Statistical
	Analysis}.
\newblock World Scientific Publishing (2004).

\bibitem[BPS06]{bib:bps06}
Boninsegni,~M., Prokof'ev,~N.V., and Svistunov,~B.V.
\newblock {\em Worm algorithm and diagrammatic Monte Carlo:  A new approach
	to continuous-space path integral Monte Carlo simulations.}
\newblock Physical Review E \textbf{74}, 036701 (2006).

\bibitem[BU07]{bib:bu07}
Betz,~V. and Ueltschi,~D.
\newblock{\em Spatial random permutations and infinite cycles}.
\newblock \texttt{arXiv:0711.1188}.
\newblock Commun. Math. Phys. 285, 469-501 (2009).

\bibitem[BU08]{bib:bu08}
Betz,~V. and Ueltschi,~D.
\newblock {\em Spatial random permutations with small cycle weights.}
\newblock \texttt{arXiv:0812.0569v1}.
\newblock Probabl. Th. Rel. Fields (2010).

%\bibitem[CGGP]{bib:cggp}
%Caracciolo,~S., Gambassi,~A., Gubinelli,~M., and Pelisetto,~A.
%\newblock {\em Finite-Size Scaling in the Driven Lattice Gas.}
%\newblock \texttt{arXiv:cond-mat/0312175}.
%\newblock Journal of Statistical Physics, vol. 115, Nos. 1/2, April 2004.

\bibitem[Ewens]{bib:ewens}
Ewens,~W.J.
\newblock {\em The sampling theory of selectively neutral alleles}.
\newblock Theor. Popul. Biol. 3, 87-112 (1972).

\bibitem[GRU]{bib:gru}
Gandolfo,~D., Ruiz,~J., and Ueltschi,~D.
\newblock {\em On a model of random cycles}.
\newblock \texttt{arXiv:cond-mat/0703315}.
\newblock Statist. Phys. 129, 663-676 (2007).

\bibitem[Kerl]{bib:kerl}
Kerl,~J.
\newblock {\em Critical behavior for the model of random spatial permutations}.
\newblock Doctoral dissertation, University of Arizona, May 2010 (in progress).

\bibitem[LB]{bib:landaubinder}
Landau,~D.P. and Binder,~K.
\newblock {\em A Guide to Monte Carlo Simulations in Statistical Physics}
	(2nd ed.).
\newblock Cambridge University Press (2005).

\bibitem[PC87]{bib:pc87}
Pollock,~E.L. and Ceperley,~D.M.
\newblock {\em Path-integral computation of superfluid densities.}
\newblock Physical Review B, vol. 36, no. 16 (1987).

\bibitem[PST98]{bib:pst98}
Prokof'ev,~N.V., Svistunov,~B.V., and Tupitsyn,~I.S.
\newblock {\em Exact, complete, and universal continuous-time worldline Monte
	Carlo approach to the statistics of discrete quantum systems.}
\newblock Journal of Experimental and Theoretical Physics, vol. 87, no. 2
	(1998).

\bibitem[PV]{bib:pv}
Pelissetto,~A. and Vicari,~E.
\newblock {\em Critical Phenomena and Renormalization-Group Theory.}
\newblock Physics Reports, 368 (6), p.549-727, Oct 2002.
\newblock \texttt{arxiv:cond-mat/0012164}.

\bibitem[S\"ut\H{o}1]{bib:suto1}
S\"ut\H{o},~A.
\newblock {\em Percolation transition in the Bose gas.}
\newblock J. Phys. A: Math. Gen. \textbf{26} (1993) 4689-4710.

\bibitem[S\"ut\H{o}2]{bib:suto2}
S\"ut\H{o},~A.
\newblock {\em Percolation transition in the Bose gas II.}
\newblock J. Phys. A: Math. Gen. \textbf{35} (2002) 6995-7002.

\bibitem[SU09]{bib:su09}
Seiringer,~R. and Ueltschi,~D.
\newblock {\em Rigorous upper bound on the critical temperature of dilute Bose
gases}.
\newblock \texttt{arXiv.org:0904.0050}.
\newblock Phys. Rev. B 80, 014502 (2009).

\bibitem[U07]{bib:qmath}
Ueltschi,~D.
\newblock {\em The model of interacting spatial permutations and its relation
to the Bose gas.}
\newblock \texttt{arXiv:0712.2443v3}.
\newblock Mathematical Results in Quantum Mechanics, pp. 225-272,
World Scientific (2008).

\end{thebibliography}

%--------------------------------------------------------------------
\end{document}
