a lot of writing

1 files changed, 354 insertions, 2 deletions

diff --git a/monte-carlo.tex b/monte-carlo.tex
index 5cf41fe..6509eab 100644
--- a/monte-carlo.tex
+++ b/monte-carlo.tex
@@ -7,6 +7,8 @@
 
 \usepackage[utf8]{inputenc}
 \usepackage{amsmath,amssymb,latexsym,mathtools,xifthen}
+\usepackage{algorithm}
+\usepackage{algorithmic}
 
 % uncomment to label only equations that are referenced in the text
 %\mathtoolsset{showonlyrefs=true}
@@ -25,6 +27,9 @@
 \def\dd{d} % differential
 \def\O{O}          % big O
 \def\o{o}          % little O
+\def\Z{\mathbb Z}
+\def\R{\mathbb R}
+\def\G{S}
 
 % subscript for ``critical'' values, e.g., T_\c
 \def\c{\mathrm c}
@@ -37,6 +42,8 @@
 \def\fS{\mathcal S}  % surface tension
 \def\fG{\mathcal G}  % exponential factor
 
+\def\H{\mathcal H}
+
 % lattice types
 \def\sq{\mathrm{sq}}
 \def\tri{\mathrm{tri}}
@@ -56,6 +63,17 @@
   \frac{\partial\tmp#2}{\partial#3\tmp}
 }
 
+\newcommand\nn[2]{\langle#1#2\rangle}
+\newcommand\avg[1]{\langle#1\rangle}
+\newcommand\eavg[1]{\langle#1\rangle_{\mathrm e}}
+\newcommand\mavg[1]{\langle#1\rangle_{\mathrm m}}
+\def\e{\mathrm e}
+\def\m{\mathrm m}
+\newcommand\new[1]{\tilde{#1}}
+
+\renewcommand\vec[1]{\mathbf{#1}}
+\newcommand\unit[1]{\hat{\mathbf{#1}}}
+
 % used to reformat display math to fit in two-column ``reprint'' mode
 \makeatletter
 \newif\ifreprint
@@ -64,7 +82,9 @@
 
 \begin{document}
 
-\title{A cluster algorithm for an Ising model in an external field}
+\algsetup{indent=1em}
+
+\title{A cluster algorithm for the Ising model in an external field}
 \author{Jaron Kent-Dobias}
 \author{James P.~Sethna}
 \affiliation{Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY, USA}
@@ -77,12 +97,344 @@
 
 \maketitle
 
+\section{Introduction}
+
+The Ising model is a simple model of a magnet. It is comprised of $n$ sites
+arranged in a regular lattice. At every site $i$ is a magnetic spin that can
+be oriented up or down, corresponding to a local magnetization $s_i$ with
+value $1$ or $-1$, respectively. Configurations of the spins are uniquely
+described by the local magnetization of every site, 
+so the model's
+configuration space is $\G^n=\G\times\cdots\times\G$ for $\G=\{-1,1\}$
+and states of the model are given by $s=(s_1,\ldots,s_n)\in\G^n$. The
+Ising Hamiltonian $\H:\G^n\to\R$ gives the energy of a particular state
+$s\in\G^n$ by
+\[
+  \H(s)=-\frac12\sum_{i=1}^n\sum_{j=1}^nJ_{ij}s_is_j-HM(s)
+  \label{eq:ham}
+\]
+where $H$ is the external magnetic field, $M:\G^n\to\R$ is the magnetization
+and is defined for $s\in\G^n$ by
+\[
+  M(s)=\sum_{i=1}^ns_i,
+\]
+and the spin--spin coupling $J$ is defined by
+\[
+  J_{ij}=\begin{cases}
+    1 & \text{$i$, $j$ are neighboring sites}\\
+    0 & \text{otherwise.}
+  \end{cases}
+\]
+The quantities of interest in statistical physics are averaged out over
+timescales large compared to the dynamics of the system. If $A:\G^n\to\R$ is a
+physical quantity that depends on the model's configuration and $s(t)$ gives
+the configuration as a function of time $t$, then
+\[
+  \bar A=\lim_{T\to\infty}\frac1T\int_0^TA(s(t))\,\dd t
+\]
+Typically the average value of physical quantities is not actually computed by
+taking a time average. Instead, the ergodic hypothesis is applied, which
+stipulates that the system's dynamics access all states $s\in S^n$ with relative
+duration $e^{-\beta\H(s)}$, so that the time average can be written as the
+ensemble average, or the average over configuration space,
+\[
+  \avg A=\frac{\sum_se^{-\beta\H(s)}A(s)}{\sum_se^{-\beta\H(s)}}
+  =\frac1Z\sum_se^{-\beta\H(s)}A(s),
+  \label{eq:avg.obs}
+\]
+where $\beta=1/kT$ is the inverse temperature in energy units and
+\[
+  Z=\sum_se^{-\beta\H(s)}
+\]
+is the partition function.
+
+Monte Carlo methods attempt to estimate quantities of the form
+\eqref{eq:avg.obs} without doing the often intractable sum explicitly by
+sampling the set of possible states such that a state $s$ is sampled with
+probably $e^{-\beta\H(s)}$, appropriately weighting each state's contribution
+to the sum. A common way of doing this is the Markov-chain Monte Carlo
+approach, which steps through configuration space using a Markov process whose
+equilibrium distribution is Boltzmann. One of the most popular of these is the
+Metropolis--Hastings algorithm, shown in Algorithm \ref{alg:met-hast}. One
+spin in the current state $s$ is randomly chosen and perturbed. The new state
+is accepted if the resulting change in energy is negative, or with probability
+$e^{-\beta\Delta\H}$ if the change is positive. Thermodynamic averages of
+functions of the state can be made by periodically sampling from this loop for
+times much longer than the autocorrelation time of the Markov process.
+
+\begin{figure}
+  \begin{algorithm}[H]
+    \begin{algorithmic}
+      \REQUIRE $s\in\G^n$
+      \LOOP
+      \STATE $s'\gets s$
+      \STATE $i \gets\text{random integer between $1$ and $n$}$
+      \STATE $s'_i\gets-s_i$
+      \IF{$e^{-\beta(\H(s')-\H(s))}\geq\text{random number between $0$ and $1$}$}
+      \STATE $s\gets s'$
+      \ENDIF
+      \ENDLOOP
+    \end{algorithmic}
+    \caption{Metropolis--Hastings}
+    \label{alg:met-hast}
+  \end{algorithm}
+\end{figure}
+
+This algorithm is very generic and works well for many models. However, near
+the critical point of the Ising model it suffers from critical slowing down,
+where the autocorrelation time of the algorithm diverges with proximity to the
+critical point.
+The Wolff algorithm is such a method that is very efficient when
+$H=0$, i.e., in the absence of an external field. It is shown in Algorithm
+\ref{alg:wolff}. Briefly, whole clusters of spins are flipped in a way that
+correctly samples states from the Boltzmann distribution while never rejecting
+a flip attempt. A random spin is selected, then flipped. Each of the spin's
+neighbors that was initially oriented the same direction as the spin is then
+flipped with probability $1-e^{-2\beta}$, and the same is true recursively for
+each of those flipped spins' neighbors.
+
+The Wolff algorithm solves many of the problems of the Metropolis--Hastings
+algorithm, but only in the absence of a field, since it depends on the natural
+symmetry $\H(s)=\H(-s)$ which is true only when $H=0$. The Wolff algorithm
+can be extended to work in a field by combining it with aspects of the
+Metropolis--Hastings algorithm in what we will refer to as the hybrid
+Wolff/Metropolis--Hastings. In each step of this algorithm, a Wolff-style
+cluster is constructed exactly as above. Instead of moving on, however, the
+system is reverted to its previous state with probability
+$\max\{0,1-e^{\beta\Delta M H}\}$, where $\Delta M$ is the change in
+magnetization due to the cluster flip. This is precisely a
+Metropolis--Hastings step where the perturbation of $s$ to $s'$ is the Wolff
+cluster and, since the detailed balance with respect to the spin--spin energy
+is performed by the cluster creation, only the change in energy due to the
+external field is considered. For $H=0$ this reverts to the Wolff algorithm
+exactly. A strictly faster version of this algorithm is accomplished by
+evaluating the probability that the cluster-building is reverted at every step
+in its formation, so that large unlikely clusters abandoned before very much
+time is wasted in constructing them. This algorithm is given as follows.
+
+\begin{figure}
+  \begin{algorithm}[H]
+    \begin{algorithmic}
+      \REQUIRE $s\in\G^n$
+      \STATE \textbf{let} $q$ be an empty stack
+      \STATE \textbf{let} $i_0$ be the index of a random site
+      \STATE $q.\mathrm{push}(i_0)$
+      \STATE $\sigma\gets s_{i_0}$
+      \WHILE {$q$ is not empty}
+        \STATE $i\gets q.\mathrm{pop}$
+        \IF {$s_i=\sigma$}
+          \FOR {$j=1$ \TO $n$}
+            \IF {$1-e^{-2\beta J_{ij}s_is_j}\geq$ random number
+          between 0 and 1}
+          \STATE $q.\mathrm{push}(j)$
+            \ENDIF
+          \ENDFOR
+          \STATE $s_i\gets-s_i$
+        \ENDIF
+      \ENDWHILE
+    \end{algorithmic}
+    \caption{Wolff (Zero-Field)}
+    \label{alg:wolff}
+  \end{algorithm}
+\end{figure}
+
+
+\begin{figure}
+  \begin{algorithm}[H]
+    \begin{algorithmic}
+      \REQUIRE $s\in\G^n$
+      \STATE \textbf{let} $q_1$ and $q_2$ be empty stacks
+      \STATE \textbf{let} $i_0$ be the index of a random site
+      \STATE $q_1.\mathrm{push}(i_0)$
+      \STATE $\sigma\gets s_{i_0}$
+      \STATE stopped $\gets$ \FALSE
+      \WHILE {$q_1$ is not empty}
+        \STATE $i\gets q_1.\mathrm{pop}$
+        \IF {$s_i=\sigma$}
+          \FORALL {$j>0$ such that $J_{ij}\neq0$}
+            \IF {$s_is_j=1$ \AND $1-e^{-2\beta J_{ij}}\geq$ random number
+          between 0 and 1}
+          \STATE $q_1.\mathrm{push}(j)$
+            \ENDIF
+          \ENDFOR
+          \STATE $s_i\gets-s_i$
+          \STATE $q_2.\mathrm{push}(i)$
+          \IF {$1-e^{-2\beta\sigma H}>$ random number between 0 and 1}
+            \STATE stopped $\gets$ \TRUE
+            \STATE \textbf{break}
+          \ENDIF
+        \ENDIF
+      \ENDWHILE
+      \IF {stopped}
+        \WHILE {$q_2$ is not empty}
+          \STATE $i\gets q_2.\mathrm{pop}$
+          \STATE $s_i\gets-s_i$
+        \ENDWHILE
+      \ENDIF
+    \end{algorithmic}
+    \caption{Hybrid Wolff/Metropolis--Hastings}
+    \label{alg:h.wolff}
+  \end{algorithm}
+\end{figure}
+
+\section{Cluster-flip in a Field}
+
+Consider a system whose state $s=\{s_0,s_1,\ldots,s_n\}$ $s\in\Z_2^{n+1}$ is specified by the
+value of $n+1$ boolean variables with the Hamiltonian
+\[
+  \new\H(s)=-\sum_{i=0}^n\sum_{j=0}^n\new J_{ij}s_is_j,
+  \label{eq:new.ham}
+\]
+where
+\[
+  \new J_{ij}=\begin{cases}
+    1 & \text{$i,j>0$, $i$, $j$ are neighboring sites}\\
+    H & \text{$i=0$ or $j=0$}\\
+    0 & \text{otherwise}
+  \end{cases}
+\]
+gives the coupling. Since \eqref{eq:new.ham} is a discrete-spin Hamiltonian
+invariant under the transformation $s\to-s$, it directly admits analysis using
+the Wolff algorithm specified above. By rewriting the sums, however, we also
+see that
+\[
+  \new\H(s)=-\frac12\sum_{i=1}^n\sum_{j=1}^nJ_{ij}s_is_j-H\new M(s)
+  \label{eq:mod_hamiltonian}
+\]
+where
+\[
+  \new M(s)=s_0\sum_{i=1}^ns_i
+\]
+This is precisely the same Hamiltonian as \eqref{eq:ham} but with a redefined
+magnetization whose sign depends on the state of the zeroth spin. Our
+algorithm for sampling the Ising model in a field, then, is simply the Wolff
+algorithm described above with $\tilde J_{ij}$ replacing $J_{ij}$. These are
+formally equivalent, since every state $s\in S$ corresponds to two states
+$\tilde s\in\tilde S$, $\tilde s=\pm(1,s)$, such that $\H(s)=\tilde\H(\tilde
+s)$. We therefore have
+\begin{align}
+  \avg A_{\tilde S}
+  &=\frac{\sum_{\tilde s}e^{-\beta\tilde\H(\tilde
+  s)}A(s_0s)}{\sum_{\tilde s}e^{-\beta\tilde\H(\tilde s)}}\\
+  &=\frac{\sum_se^{-\beta\tilde\H(s_0,s)}A(s)+\sum_se^{-\beta\tilde\H(-s_0,-s)}A(s)}{\sum_se^{-\beta\tilde\H(1,s)}+\sum_se^{-\beta\tilde\H(-1,-s)}}\\
+  &=\avg A_S
+\end{align}
+$\tilde A(s_0,s)=A(s_0s)$. This is true for the mangetization $\tilde M$ as
+defined above, for instance.
+
+\begin{figure}
+  \begin{algorithm}[H]
+    \begin{algorithmic}
+      \REQUIRE $s\in\G^n$
+      \STATE \textbf{let} $q$ be an empty stack
+      \STATE \textbf{let} $i_0$ be the index of a random site
+      \STATE $q.\mathrm{push}(i_0)$
+      \STATE $\sigma\gets s_{i_0}$
+      \WHILE {$q$ is not empty}
+        \STATE $i\gets q.\mathrm{pop}$
+        \IF {$s_i=\sigma$}
+          \FOR {$j=0$ \TO $n$}
+            \IF {$1-\exp(-2\beta \tilde J_{ij}s_is_j)\geq$ random number
+          between 0 and 1}
+          \STATE $q.\mathrm{push}(j)$
+            \ENDIF
+          \ENDFOR
+          \STATE $s_i\gets-s_i$
+        \ENDIF
+      \ENDWHILE
+    \end{algorithmic}
+    \caption{Wolff (Non-zero Field)}
+    \label{alg:wolff-field}
+  \end{algorithm}
+\end{figure}
+
+
+At sufficiently small $H$ both our algorithm and the hybrid Wolff--Metropolis reduce
+to the ordinary Wolff algorithm. At very large $H$, the hybrid
+Wolff--Metropolis behaves exactly like Metropolis, where almost only one spin
+is ever flipped at a time with probability $\sim e^{-2\beta H}$, since the
+energy is dominated by the field.
+
+We measured the correlation time $\tau$ of these three algorithms at various fields
+and temperatures. This was done using a batch mean. Time was measured as ``the
+number of spins that have either been attempted to flip or attempted to add to
+a cluster.''
+
+$E_s=A_s+XB_s$
+\[
+  \pd{\avg B}X=\beta\big(\avg{B^2}-\avg B^2\big)
+\]
+\[
+  \chi=\pd{\avg M}H=\beta\big(\avg{M^2}-\avg M^2\big)
+\]
+
+Since \eqref{eq:mod_hamiltonian} has the form of a simple (ferromagnetic)
+spin-spin Hamiltonian, its ensemble average can be taken using the Wolff
+algorithm while never rejecting a step.
+
+
+\section{Magnetization Estimator}
+
+For any finite system, the average magnetization at zero field is identically
+zero at all temperatures. However, people often want to use finite simulations
+to estimate the average magnetization of the Ising model in the thermodynamic
+limit, where due to a superselection principle there is one of two average
+magnetizations at zero field below the critical temperature. The ergodic
+principle that defends the esemble average replacing a time average is no
+longer valid, since once the system has taken a certain magnetization it has
+zero probability of ever changing direction. This is typically
+accomplished by using $\avg{|M|}$. But what is the best way to estimate the
+magnetization of a finite system in the presence of a field?
+
+$\tilde M$ is the true magnetization of the finite-size system. We can
+estimate the magnetization of system in the thermodynamic limit by enforcing
+the superselection principle:
+\[
+  S_\e=\{s\mid M(s)>0\}
+\]
+\[
+  \eavg{A}=\frac{\sum_{s\in S_\e}e^{-\beta\H(s)}A(s)}{\sum_{s\in
+  S_\e}e^{-\beta\H(s)}}=\frac1{Z_\e}\sum_{s\in S_\e}e^{-\beta\H(s)}A(s)
+\]
+\[
+  S_\m=\{s\mid M(s)<0\}
+\]
+\[
+  \mavg{A}=\frac{\sum_{s\in S_\m}e^{-\beta\H(s)}A(s)}{\sum_{s\in
+  S_\m}e^{-\beta\H(s)}}=\frac1{Z_\m}\sum_{s\in S_\m}e^{-\beta\H(s)}A(s)
+\]
+For $H=0$, $\eavg M=\avg{|M|}$, defending the canonical measure of
+the ferromagnetic phase magnetization at zero field. This can be seen first by
+identifying the bijection $s\to-s$ that maps $S_\e$ to $S_\m$. Then, since
+$\H(s)=\H(-s)$, $M(s)=-M(-s)$ when $H=0$, $Z_\e=Z_\m$, $Z=2Z_\e+Z_0$, and
+\begin{align}
+  \avg{|M|}
+  &=\frac1{2Z_\e+Z_0}\sum_{s\in S}e^{-\beta\H(s)}M(s)\\
+  &=\frac1{Z_\e+\frac12Z_0}\sum_{s\in S_\e}e^{-\beta\H(s)}M(s)\\
+  &=\frac1{1+\frac{Z_0}{2Z_\e}}\eavg M
+\end{align}
+At infinite temperature, $Z_0/Z_\e\simeq N^{-1/2}\sim L^{-1}$ for large $L$,
+$N$. At any finite temperature, especially in the ferromagnetic phase,
+$Z_0\ll Z_\e$ in a much more extreme way.
+
+If the ensemble average over only positive magnetizations can be said to
+converge to the equilibrium state of the Ising model in the thermodynamic
+limit, what of the average over only \emph{negative} magnetizations,
+ordinarily unphysical? 
+
+Collapse of the magnetizations $\avg M$ and $\eavg M$ at low and
+high temperature
 
 \begin{acknowledgments}
-  Thatks!
+  Thanks!
 \end{acknowledgments}
 
 %\bibliography{monte-carlo}
 
+
+\begin{itemize}
+  \item metastable and stable state mixing at finite size?
+  \item tracking the wolff algorithm cluster size versus correlation length
+\end{itemize}
 \end{document}