diff options
| -rw-r--r-- | monte-carlo.tex | 356 |
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} |