
% Template for Elsevier CRC journal article
% version 1.0 dated 13 October 2009

% This file (c) 2009 Elsevier Ltd.  Modifications may be freely made,
% provided the edited file is saved under a different name

% This file contains modifications for Procedia Computer Science
% but may easily be adapted to other journals

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% This template uses the elsarticle.cls document class and the extension package ecrc.sty
%% For full documentation on usage of elsarticle.cls, consult the documentation "elsdoc.pdf"
%% Further resources available at http://www.elsevier.com/latex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The '1p' and 'times' class options of elsarticle are used for Elsevier CRC
\documentclass[1p,times]{elsarticle}

%% The `ecrc' package must be called to make the CRC functionality available
\usepackage{ecrc}

%% The ecrc package defines commands needed for running heads and logos.
%% For running heads, you can set the journal name, the volume, the starting page and the authors

%% set the volume if you know. Otherwise `00'
\volume{00}

%% set the starting page if not 1
\firstpage{1}

%% Give the name of the journal
\journalname{Physics Procedia}

%% Give the author list to appear in the running head
%% Example \runauth{C.V. Radhakrishnan et al.}
\runauth{J. Kerl}

%% The choice of journal logo is determined by the \jid and \jnltitlelogo commands.
%% A user-supplied logo with the name <\jid>logo.pdf will be inserted if present.
%% e.g. if \jid{yspmi} the system will look for a file yspmilogo.pdf
%% Otherwise the content of \jnltitlelogo will be set between horizontal lines as a default logo

%% Give the abbreviation of the Journal.
\jid{procs}

%% Give a short journal name for the dummy logo (if needed)
\jnltitlelogo{Physics Procedia}

%% Hereafter the template follows `elsarticle'.
%% For more details see the existing template files elsarticle-template-harv.tex and elsarticle-template-num.tex.

%% Elsevier CRC generally uses a numbered reference style
%% For this, the conventions of elsarticle-template-num.tex should be followed (included below)
%% If using BibTeX, use the style file elsarticle-num.bst

%% End of ecrc-specific commands
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The amssymb package provides various useful mathematical symbols
% JRK 2010-02-17 \usepackage{amssymb}
\usepackage{amssymb,amsmath,amsthm}
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%% \usepackage{lineno}

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%
%% \biboptions{comma,round}

% \biboptions{}

% if you have landscape tables
\usepackage[figuresright]{rotating}

% put your own definitions here:
%   \newcommand{\cZ}{\cal{Z}}
%   \newtheorem{def}{Definition}[section]
%   ...
\input{kerl_syms}

\usepackage[mathscr]{eucal}
\numberwithin{equation}{section}
\usepackage{graphicx}
\usepackage{psfrag}
\usepackage{afterpage}
\usepackage{subfigure}
\usepackage{floatflt}

% add words to TeX's hyphenation exception list
%\hyphenation{author another created financial paper re-commend-ed Post-Script}

% declarations for front matter

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for the associated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for the associated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for the associated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%%
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}
\dochead{}
%% Use \dochead if there is an article header, e.g. \dochead{Short communication}

\title{A worm algorithm for random spatial permutations}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{<author name>}
%% \address[label1]{<address>}
%% \address[label2]{<address>}

\author{John Kerl}

\address{PhD candidate, University of Arizona Department of Mathematics
(\texttt{kerl@math.arizona.edu})}

\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.  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}

\begin{keyword}
random spatial permutations \sep
Markov chain Monte Carlo \sep
Bose gas \sep
worm algorithm
%% keywords here, in the form: keyword \sep keyword

%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

\end{keyword}

\end{frontmatter}

%%
%% Start line numbering here if you want
%%
% \linenumbers

%% ================================================================
%% ================================================================
%% main text
%% ----------------------------------------------------------------
\section{The probability model}

The \emph{model of random spatial permutations}, or \emph{RSP model}, 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:su09,bib:baymetal} 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 the author's 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, and 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, bb = 0 0 200 110]{pmt1.pdf}
\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}.  This note focuses solely on a particular technique
for estimating expectations of random variables in the RSP model.

%% ----------------------------------------------------------------
\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 the author's simulations
and dissertation: thermalization time, proofs of detailed balance,
autocorrelation, quantification of variance of sample means, and
finite-size-scaling analysis of finite-lattice computational results.  (See
\cite{bib:berg,bib:landaubinder,bib:cggp} for background on such methods as
applied to models other than the RSP model.)

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}
\includegraphics[scale=0.5, bb = 0 0 400 100 ]{gkcycles.pdf}
\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}
\includegraphics[scale=0.6, bb = 0 0 500 120]{gkwinding.pdf}
\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
can 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}
\includegraphics[scale=0.3, bb = 0 0 400 120]{worm1slide.pdf}
\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 define an \emph{extended
random cycle model} with 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.  This is
justified by a theorem of \cite{bib:kerl}:  letting $\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 not only for the cycle-weight
Hamiltonian given above, but for any Hamiltonian on random permutations.

%% ----------------------------------------------------------------
\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 as
above.  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}
\includegraphics[scale=0.60, bb = 0 0 500 170]{wormmoves2x2.pdf}
\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.  SAR and worm CPU times are both approximately $aN + bN^2$.
For the SAR algorithm, $b$ is tiny; for the worm algorithm, it 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}
\includegraphics[scale=0.56, bb = 0 0 320 200]{sar_perf_small_L.pdf}
\includegraphics[scale=0.56, bb = 0 0 320 200]{worm_perf_small_L.pdf}
\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.
First, band updates, discussed above, have too-low acceptance rate.  Second,
one might temporarily pinch the torus geometry in the SAR algorithm, such that
the distance penalty for wrapping around the torus is decreased.  Third, one
may 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.
Fourth, the worm algorithm might be modified to somehow direct the worm.  For
example, 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.  Last, one may review
the PIMC literature again and seek other inspiration.

%% ================================================================
%% ================================================================




%% ----------------------------------------------------------------
%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% ----------------------------------------------------------------
%% References
%%
%% Following citation commands can be used in the body text:
%% Usage of \cite is as follows:
%%   \cite{key}         ==>>  [#]
%%   \cite[chap. 2]{key} ==>> [#, chap. 2]
%%

%% References with BibTeX database:

%\bibliographystyle{elsarticle-num}
%\bibliography{kerl_csp_2010_bibdb.bib}

%% Authors are advised to use a BibTeX database file for their reference list.
%% The provided style file elsarticle-num.bst formats references in the required Procedia style

%% For references without a BibTeX database:
%
\begin{thebibliography}{00}

\bibitem{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{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{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{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).

\bibitem{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{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{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{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{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{bib:ewens}
Ewens,~W.J.
\newblock {\em The sampling theory of selectively neutral alleles}.
\newblock Theor. Popul. Biol. 3, 87-112 (1972).

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

\bibitem{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{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{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{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{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).

\end{thebibliography}
\end{document}

%% ================================================================

%%
%% End of file `procs-template.tex'. 
