path: root/monte-carlo.tex

%
%  Created by Jaron Kent-Dobias on Thu Apr 20 12:50:56 EDT 2017.
%  Copyright (c) 2017 Jaron Kent-Dobias. All rights reserved.
%
\documentclass[aps,prl,reprint]{revtex4-1}

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

% I want labels but don't want to type out ``equation''
\def\[{\begin{equation}}
\def\]{\end{equation}}

% math not built-in
\def\arcsinh{\mathop{\mathrm{arcsinh}}\nolimits}
\def\arccosh{\mathop{\mathrm{arccosh}}\nolimits}
\def\ei{\mathop{\mathrm{Ei}}\nolimits} % exponential integral Ei
\def\re{\mathop{\mathrm{Re}}\nolimits}
\def\im{\mathop{\mathrm{Im}}\nolimits}
\def\sgn{\mathop{\mathrm{sgn}}\nolimits}
\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}

% scaling functions
\def\fM{\mathcal M}  % magnetization
\def\fX{\mathcal Y}  % susceptibility
\def\fF{\mathcal F}  % free energy
\def\fiF{\mathcal H} % imaginary free energy
\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}}
\def\hex{\mathrm{hex}}

% dimensions
\def\dim{d}
\def\twodee{\textsc{2d} }
\def\threedee{\textsc{3d} }
\def\fourdee{\textsc{4d} }

% fancy partial derivative
\newcommand\pd[3][]{
  \ifthenelse{\isempty{#1}}
    {\def\tmp{}}
    {\def\tmp{^#1}}
  \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
\@ifclasswith{revtex4-1}{reprint}{\reprinttrue}{\reprintfalse}
\makeatother

\begin{document}

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

\date\today

\begin{abstract}
  An abstract.
\end{abstract}

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