From kea@sampson.ccsf.caltech.edu Mon Oct  3 15:23:10 1994
Date: Fri, 30 Sep 1994 18:25:23 -0700
From: Terri Canzian <kea@sampson.ccsf.caltech.edu>
To: paulc@npac.syr.edu
Subject: Is This What You're Looking For?


\documentstyle{article}
\def\gea{\mathrel{\raise .4ex\hbox{rlap{$>$}\lower 1.2ex\hbox{$\sim$}}}}
\begin{document}

\begin{center}
{\bf Monte Carlo Algorithms}
\end{center}

So, we have a valid Monte Carlo algorithm if:
\begin{itemize}
\item we have a means of generating a new configuration $B$ from a previous 
configuration $A$ such that the transition probability $W(A\rightarrow B)$ 
satisfies detailed balance.

\item the generation procedure is $\underline{\rm ergodic}$, i.e., every 
configuration can be reached from every other configuration in a finite number 
of iterations.
\end{itemize}

The Metropolis algorithm satisfies the first criterion for $\underline{\rm 
all}$ statistical systems.  The second criterion is model dependent, and not 
always true (e.g., at $T=0$).

\begin{center}
{\bf Metropolis Algorithm for the Ising Model}
\end{center}

For the Ising model, the change in the configuration is to try to change 
(``flip'') a single spin.  If we try to change many spins, $\Delta E$ will be 
large, so the probability to accept the change, $e^{-\Delta E/kT}$, will be 
small.

For a single spin flip, $\Delta E$ depends only on the spin values at the site 
and its nearest neighbors, i.e., update is $\underline{\rm local}$.

If $\Delta E \le 0$, we make the proposed change.

If $\Delta E > 0$, we make the change with probability $e^{-\Delta E/kT}$.

\begin{center}
{\bf Programming the Metropolis Algorithm on a Computer}
\end{center}

Update a single site at a time.  Choose sites at random (In most cases, this 
is not really necessary, and a simple contiguous loop over sites is 
sufficient).

A loop over all sites is referred to as one Metropolis $\underline{\rm 
sweep}$.

The probabilistic part of the algorithm is done using a random number 
generator (hence the name---Monte Carlo).

Generate a random number $r$ in [0,1).  A process with probability $p$ is 
performed if $r < p$.

The Metropolis algorithm is very simple to program.

\begin{center}
{\bf A Monte Carlo Simulation}
\end{center}

We generate configurations of the Ising model using a Monte Carlo algorithm, 
such as Metropolis.  Then what do we do with them.  We want to measure 
quantities of interest.  For instance,

\begin{tabular}{ll}
energy per spin & $E = - \frac{J}{V} \sum_{\langle i,j\rangle} S_i S_j$ \\
magnetization per spin & $M = \frac{1}{V} \sum_i S_i$ \hspace{1em} $V$=volume 
(\# sites) \\
correlation function & $\Gamma (n) = \frac{1}{V} \sum_i S_i S_{i+n} - \left( 
\frac{1}{V} \sum_i S_i \right)^2$ \\
correlation length $\xi$ & $\sim e^{-n/\xi}$ \\
specific heat & $\frac{C}{V} = \frac{1}{V} \frac{\delta E}{\delta T} = 
\left\langle \left( E - \langle E\rangle \right)^2 \right\rangle = 
\left\langle E^2\right\rangle - \langle E\rangle^2$ \\
susceptibility & $\frac{x}{V} = \frac{1}{V} \frac{\delta M}{\delta T} =  
\left\langle \left( M - \langle M\rangle \right)^2 \right\rangle = 
\left\langle M^2\right\rangle - \langle M\rangle^2$ \\
\end{tabular}

\begin{center}
{\bf Phase Transitions}
\end{center}

Phase transitions are described in part by their $\underline{\rm order}$:

\begin{itemize}
\item First order---discontinuity in the order parameter (e.g., magnetization) 
or energy, i.e., first derivative of the partition function.

\item Second order---divergence in the susceptibility or specific heat, i.e., 
second derivative of partition function.

\item Third order, etc.
\end{itemize}

Order parameter differentiates the two phases---usually zero in disordered 
phase, non-zero in order phase.

First order transitions are classified by quantities, such as the 
$\underline{\rm latent\ heat}$, the energy change at the critical point $T_c$.

\begin{center}
{\bf Critical Exponents}
\end{center}

Second order transitions are classified by their $\underline{\rm critical\ 
exponents}$, which measure how quantities diverge at the critical point.  
N.B.\ These are asymptotic $(T\rightarrow T_c)$ results.

\begin{eqnarray*}
M & \sim & \vert T-T_c \vert^\beta \\
X & \sim & \vert T-T_c \vert^\gamma \\
C & \sim & \vert T-T_c \vert^\alpha \\
xi & \sim & \vert T-T_c \vert^\nu \\
\end{eqnarray*}

In many cases, these exponents are $\underline{\rm universal}$, i.e., they do 
not depend on details of the model, but only on gross features such as the 
dimension of the space, the symmetries of the Hamiltonian.  This explains the 
quantitative success of very simple spin models, e.g., Ising.

\begin{center}
{\bf Problems, Subtleties, and Tricks of the Trade}
\end{center}

\begin{enumerate}
\item Thermalization

Must be sure that Markov process has reached ``equilibrium.''  Approach to 
equilibrium is exponential $\sim e^{t/\tau}$, but $\tau$ can be very large.  
We may do 100,000 sweeps at a certain $T$, but if $\tau$ is of the same order, 
then the results are bogus.

Check:  Start from hot (random) configuration and cool.  Start from cold 
(ordered) configuration and heat.  When the two have converged to about the 
same results, then the system is equilibrated.  (But what about first order 
P.T.s at $T_c$?)

Alternatively: Check to see whether measurements are converging to a constant 
average value.  Can do this by $\underline{\rm binning}$\@: break measurement 
values up into a number of large continuous (in time) segments (bins), and 
average over each bin.  Throw away results from early deviant bins. (Can also 
use binning to calculate errors---see later).

\item Ergodicity

Can be very difficult to prove.  Can get trapped in sub-sets of configuration 
space, periodic loops, etc.

$\underline{\rm Quasi-ergodic}$\@: may be a high energy barrier between 
subspaces, so it's possible to reach all states, but $\underline{\rm highly\ 
improbable}$ in a finite simulation.

Check:  Multiple runs from different initial configurations.

\item Fitting

Once we have measured, for example, specific heat $C$, need to do a multi-
parameter fit to data which usually has big errors.  This is very difficult!
\begin{displaymath}
C = A + B \left\vert T-T_c \right\vert^{-\alpha}
\end{displaymath}
This is a four parameter fit, so usually have so much freedom that a wide 
range of values are possible.  Usually need to assume $A\approx 0$, and get 
$T_c$ from other measurements.

Also, the above formula is only true asymptotically close to $T_c$.  But how 
close is that?  Not always clear.  May see ``apparent'' or ``pseudo-critical'' 
exponents, which change or disappear as get closer to $T_c$.

\item Finite Size Effects

When simulating any macroscopic system with a very large number of degrees of 
freedom, invariably make an approximation to deal with a smaller and/or 
discretized model system.  This introduces errors, or finite size effects.  
Have to understand these, and be able to extrapolate to an infinite system, 
usually by doing a number of simulations at different lattice sizes.

For spin models, have a finite lattice of $V=L^d$ sites.  But only get true 
phase transition (i.e., divergence) when $L\rightarrow\infty$.  For finite 
system, have rounded peaks rather than divergences.

\begin{center}
{\bf Finite Size Scaling}
\end{center}

In some instances, such as at second order phase transitions, the asymptotic 
(large $L$) behavior is known.

\begin{eqnarray*}
C_{\max} (L) & \sim & L^{\alpha/\nu} \\
\left\vert T_c (L) - T_c (\infty )\right\vert & \sim & L^{1/\nu} \\
M\left( T_c \right) & \sim & L^{-\beta/\nu} \\
X\left( T_c \right) & \sim & L^{\alpha/\nu} \\
\end{eqnarray*}

Simple ``proof'':
\begin{eqnarray*}
\xi_{\max} (L) & \sim & L \hspace{1em} \mbox{correlation length must be less 
than length of lattice!} \\
\xi_{\max} (L) & \sim & \left\vert T_c (L) - T_c (\infty )\right\vert^{-\nu} \\
\mbox{Therefore} & & \left\vert T_c (L) - T_c (\infty )\right\vert \sim 
L^{1/\nu} \\
C_{\max} (L) & \sim & \left\vert T_c (L) - T_c (\infty )\right\vert^{-\alpha} \\
& \sim & L^{\alpha/\nu} \\
\end{eqnarray*}

These relations are extremely useful, and even allow us to extract exponents 
and $T_c (\infty )$.  But only if we are at ``large enough'' $L$.  If $\xi$ is 
very large (but not infinite), may get pseudo-exponents.  This is often the 
case at ``weakly'' first order transition (very small latent heat).  So it is 
often very difficult to distinguish between a first and second order phase 
transition, or even to say for sure if there is a phase transition at all!  
(Is $\xi = \infty$, or just $xi >> L$?).

In general, don't know the form of the finite size scaling, so just try to fit 
data to a $\underline{\rm scaling\ function}$ or go to large enough systems so 
that results $\rightarrow$ constant.

\item Boundary Conditions

Generally, use ``periodic boundary conditions,'' where the lattice wraps 
around on itself to form a torus.  This has been shown to give smallest finite 
size effects.

\item Acceptance

For the Ising model, there is only one possible change to the state\@: $S_i 
\rightarrow -S_i$.  More generally (e.g., continuous spin models) we are free 
to make the change as large or small as we like.

A small change means $\Delta E$ is small and so the $\underline{\rm 
acceptance}$ (proportion of proposed changes which are accepted) is near 1.  
But we move slowly through configuration space.  Large changes usually have 
large $\Delta E$ and thus low acceptance.  Aim for acceptance $\approx 0.5$.

\item Pseudo-random Numbers

Generating ``good'' random numbers can be a major problem in MC simulations.  
Many basic random number generators supplied with compters (especially 
parallel computers) are not sufficient.  Particularly true in the past, but 
unfortunately still true far too often today!

There are many tests for uniform distribution, correlations, etc.  Monte Carlo 
measurement of energ of Ising model is a good test---result is known exactly 
for $L^2$ lattice.  Some RNGs pass statistical tests, but fail the Ising model 
test!

Even some ``good'' RNGs are not good enough because their period of repetition 
is too small.

For example, multiplicative linear congruential generators
\begin{displaymath}
r_{i+1} = (a * r_i + c) \mbox{mod}\ m \\
\end{displaymath}
``good'' values of $a$, $c$, $m$ are known.  But for 32-bit word period is at 
$\underline{\rm most}$ $2^{31} \approx 10^{10}$.  Need period to be $>>$ 
number of random numbers used in simulation.

Big Ising simulation, $1024^2$ lattice, 10 million sweeps, uses $> 10^{13}$ 
random numbers.

Using 64-bit words, or combining two 32-bit MLCGs, can give periods $\gea 
10^{18} \approx$ \# nanoseconds per year.

Others (lagged Fibonacci, etc.) have even longer periods.

\item Critical Slowing Down

A problem with MC simulation is that the error decreases only as $1/\sqrt{N}$ 
(i.e., need to do 100 times more work to get one more significant figure 
accuracy).

However, things are even worse.  Error is only $1/\sqrt{N}$ for $N$ 
$\underline{\rm independent}$ measurements.  But MC configurations are usually 
correlated---$C_{i+1}$ depends on $C_i$.  Measure the autocorrelation function 
(similar to $\Gamma (n)$ except in computer time, not real space); for an 
operator (measurement) $A$:

\begin{eqnarray*}
\rho(t) & = & \frac{\langle A(t) A(0)\rangle - \langle A(0)\rangle^2}{\langle 
A(0)^2 \rangle - \langle A(0)\rangle^2} \\
& \sim & e^{-t/\tau_{\exp}} \\
\end{eqnarray*}

$\tau_{\exp}$ is the (exponential) auto-correlation time, and is equivalent to 
the $\tau$ which measures equilibration time of the Markov process.  Need to 
thermalize for time $t>> T_{\exp}$.

Integrated autocorrelation time
\begin{displaymath}
\tau_{\rm int} = \frac{1}{2} + \sum^\infty_{t=1} \rho (t)
\end{displaymath}
(normalized so $\tau_{\exp} \approx \tau_{\rm int}$ if $\rho\approx 
e^{-t/\tau}$)
is a measure of how correlations affect statistical errors.  Can show that
\begin{eqnarray*}
\mbox{error in mean} & = & \sqrt{\frac{2\tau_{\rm int}}{N}} \cdot \sigma \\
\sigma & = & \sqrt{\rm variance} \\
\end{eqnarray*}
Problem for MC is that $\tau$s $\underline{\rm diverge}$ at $T_c$!!  Critical 
slowing down.

To decorrelate a configuration, we need to make changes on a scale $>$ the 
spatial correlation length of the system.  For $\underline{\rm local}$ MC 
algorithms, such as Metropolis, the effect of a change is only seen by nearest 
neighbors.  Since the update is a stochastic process, the changes will 
propagate through the system like a random walk.  So the number of iterations 
required to propagate a distance $\xi$ goes like $\xi^2$.  At the critical 
point, $\xi \sim L$, so we have that $\tau \sim L^2$.  This is a major problem 
for MC simulations, since we need large $L$.
\end{enumerate}

\begin{center}
{\bf Dynamic Critical Exponents}
\end{center}

More generally, we have $\tau (L) \sim L^z$ where $z$ is called the dynamic 
critical exponent.

For local algorithm, such as Metropolis, $z$ is almost always $\approx 2$.

One of the most active areas of research in the field of MC simulation is 
trying to find new algorithms with $z<2$.  These are usually non-loca, multi-
scale algorithms, c.f. multigrid methods for solving PDEs, where the same kind 
of problem occurs.

\begin{center}
{\bf Acceleration Algorithms}
\end{center}

A number of different methods have been used to attempt to ``accelerate'' the 
dynamics of Monte Carlo updates.
\begin{itemize}
\item Over-relaxation
\item Multigrid
\item Fourier acceleration
\item Cluster algorithms
\end{itemize}

The first two are commonly used to improve numerical solutio of PDEs.  In some 
cases, or can reduce $z$ to $\approx 1$.  Multigrid Monte Carlo works well in 
certain cases. Cluster algorithms work very well, but currently are only 
applicable to a small class of spin models.

\begin{center}
{\bf Cluster Algorithms}
\end{center}

The idea behind cluster MC algorithms is very simple.  Local MC algorithms 
perform poorly because they update only one spin at a time.  Cluster 
algorithms use a clever way of finding large clusters of sites which can be 
flipped at once.  (Remember that trying to flip a group of sites chosen at 
$\underline{\rm random}$ gives a large $\Delta E$ and thus a small 
acceptance).

Cluster algorithms were invented by Swendsen and Wang in 1987, for Potts spin 
models.  They have since been generalized to other models, and other 
algorithms.

\begin{center}
{\bf The Potts Model}
\end{center}

The Potts model is an obvious generalization of the Ising model to an 
arbitrary number, $Q$, of discrete states.  Proposed by Domb to his grad 
student, Potts.

\begin{eqnarray*}
E & = & -J \sum_{\langle i,j\rangle} \delta \left( s_i, s_j \right) \\
\delta & = & 1 \hspace{1em} \mbox{if}\ s_i = s_j \rightarrow E = -J \\
& = & 0 \hspace{1em} \mbox{if}\ s_i \not= s_j \rightarrow E = 0 \\
\end{eqnarray*}
2-state Potts $\Leftrightarrow$ Ising

$Q>2$ model has different critical exponents, interesting phase structure.  $Q 
\ge 5$ has first order transition in two dimensions.

\begin{center}
{\bf Percolation}
\end{center}

Swendsen and Wang derived their algorithm from a simple relation which 
transforms a Potts spin model into a bond percolation model.

What is percolation?

Suppose we have a two-dimensional lattice of sites, and we put a bond between 
every pair of neighboring sites with probability $p$.  Sites are bonded 
together into clusters.  At a critical value of $p$, will get $\underline{\rm 
percolation}$, i.e., one cluster will span the lattice.  For an infinite two-
dimensional lattice, $p_c = 1/2$.

\begin{center}
{\bf The Swendsen-Wang Algorithm}
\end{center}

SW found that by placing bonds between neighboring sites with the same spin 
value, the Potts model can be written in terms of interacting $\underline{\rm 
clusters}$ of spins.  For a certain choice of $p$,
\begin{displaymath}
p = 1 - e^{-\beta}
\end{displaymath}
the interaction energy is zero, so the clusters are independent.  This means 
we can update all the clusters (by changing the current spin values to a 
random new spin value) completely independently!

Since the clusters can be very large, the update is non-local, and $z$ can be 
very small (even zero!).

\begin{center}
{\bf The Wolff Algorithm}
\end{center}

Wolff introduced a variant of the SW algorithm, in which a site is chosen at 
random, and a single cluster is grown around the site.  This tends to favor 
larger clusters, and thus smaller autocorrelation times.  Wolff also 
generalized the cluster algorithms to continuous spin models.

Cluster algorithms work amazingly well for some models, e.g., two-dimensional 
Ising, where $z=0$.  Even for large $512^2$ lattice, $\tau\approx 10$, which 
is 0(1000) times smaller than $\tau$ for Metropolis.

\begin{center}
{\bf Connected Component Labeling}
\end{center}

The cluster algorithms are very different to Metropolis.  Main computational 
task is to take local information (bond connetions) and work out global 
information (clusters of sites).  This is an example of $\underline{\rm 
connected\ component\ labeeling}$.

Component labeling is commonly used in image processing, to join neighboring 
pixels into connected regions which are the shapes in the image.  Very 
efficient $\underline{\rm sequential}$ algorithms exist for this procedure.

\begin{center}
{\bf Parallel Algorithms}
\end{center}

$\underline{\rm Metropolis}$ algorithm for spin models is a clasic example of 
a problem which is easy to parallelize.  Regular data structures (simple 
domain decomposition, perfect load balancing), local interactions (local 
communication), lots of parallelism.

$\underline{\rm Swendsen-Wang}$ algorithm is a classic example of a problem 
which is very difficult to parallelize efficiently.  Clusters are irregular 
(fractal!) in size and shape (hard to load balance, lots of non-local 
communication).

$\underline{\rm Wolff}$ algorithm is a classic example of a problem which is 
virtually impossible to parallelize efficiently!

\begin{center}
{\bf Topics to be Covered in More Detail in Later Lectures}
\end{center}

\begin{itemize}
\item Parallel algorithms

Metropolis and Swendsen-Wang algorithms for SIMD and MIMD machines, coarse and 
fine grain machines.  Independent parallelism in MC simulations.  Parallel 
component labeling algorithms.

\item QCD and lattice field theory

How can we solve quantum field theories on a lattice, and why do we need 
TeraFLOPS to do it?

\item Random surfaces and quantum gravity and strings

Simulation of quantum gravity and superstrings using dynamically triangulated 
random surfaces---updating the lattice as well as the fields.  How do we do 
this in parallel?

\item Optimization

How can we use statistical mechanics and Monte Carlo to solve (approximately) 
general optimization problems.  Simulated annealing.

\item Other applications

Study of complex and disordered systems---polymers, protein folding, spin 
glasses and frustration, etc. Random numbers?
\end{itemize}

\begin{center}
{\bf Parallel Metropolis Algorithm}
\end{center}

The Metropolis algorithm for a spin model is well suited to parallelism, since 
it is regular and local.  The update of a site depends only on its nearest 
neighbors.

We cannot update a site and its neighbors at the same time, since they are 
dependent, and we would violate detailed balance.  However, we can partition 
the lattice into a checkerboard, and update all the black sites at once, since 
they are all independent.

\begin{center}
{\bf Red/Black or Checkerboard Update}
\end{center}

So, in order to satisfy detailed balance (and thus be guaranteed a valid monte 
Carlo procedure), we partition the lattice into a checkerboard of alternating 
red and black sites (note we can also do this for $d>2$ dimensions).  We can 
then update all the black sites in parallel, followed by all the red sites in 
parallel.

Note that this is a regular, SIMD, algorithm.  The checkerboard update is done 
with the help of a logical $\underline{\rm mask}$ (BLACK) which only activates 
the appropriate (virtual) processors which deal with black sites.  We then 
change the mask (BLACK=.NOT.BLACK) and update the other sites.

\begin{center}
{\bf Points to Note:}
\end{center}

\begin{itemize}
\item Use standard domain decomposition.
\item Generally want periodic boundary conditions, so use periodic shift 
operations (CSHIFT).
\item Very simple to code in Fortran 90.
\end{itemize}

Since this problem is simple, local, and regular, it is generally possible to 
parallelize very efficiently (even for the Ising model, where the calculation 
is very simple), $\underline{\rm provided}$ the number of sites per processor 
(VP ratio) is greater than 1.  Otherwise, the efficiency is $< 50$\%, since 
half the processors are idle.  Note that if the VP rati is $>2^d$, we no 
longer need bother with the checkboard scheme, since we no longer have the 
possibility of updating neighboring sites at once.

Notice that this algorithm, when written in Fortran with message passing 
(e.g., Express Fortran), looks $\underline{\rm exactly}$ the same as the 
sequential algorithm!

The only difference is in how the shift subroutine, which finds the 
neighboring spin values, is implemented.

For the sequential algorithm, shift just handles the periodic boundary 
conditions.

For the parallel algorithm, shift handles the boundary conditions and does 
message passig if the neighboring spin is on a different processor.

Also, note that here size is the size of the sub-lattice on each processor, 
not the total lattice size.

\begin{center}
{\bf Efficiency}
\end{center}

For the parallel Metropolis algorithm on a two-dimensional mesh of $N=P\times 
P$ processors, we use standard domain decomposition to distribute the 
$V=L\times L$ sites of a two-dimensional lattice over the processor grid, so 
that every processor has an $l\times l$ sub-lattice $(l=L/P)$.  The 
communication time for this local algorithm will be proportional to $l$ (the 
number of sites on the edge of the sub-lattice).

The calculation time will be proportional to $l^2$ (the number of sites 
covering the area or volume of the sub-lattice).  Thus, as long as $l$ is 
large enough, the communication/calculation ratio will be small, and the 
efficiency will be near 1.

\begin{center}
{\bf Coarse Grained Algorithm}
\end{center}

On a coarse grained message passing machine, there is a more efficient 
parallel Metropolis algorithm.  Again, we do a checkerboard partitioning of 
the lattice.  But now we do the data communicatin as large blocks of spins 
rather than single spins, by passing all the edge data to neighboring 
processors after each red or black update.  This reduces the latency time for 
communication enormously, and greatly improves efficiency.

\begin{center}
{\bf Measurements}
\end{center}

Measurements of standard quantities such as energy and magnetization are 
easily done in parallel.  Just calculate the quantities locally (e.g., energy 
per site), which can be done completely in parallel (no checkerboarding 
required), sum over sites in each processor, and then combine the partial sums 
on each processor, using for example SUM in CMFortran (Fortran 90) and 
EXCOMBINE in Express C and KXCOMB in Express Fortran (explicit message passing 
languages).  For example, energy$=$ SUM (energy-per-site)/volume.

\begin{center}
{\bf Non-local Measurements}
\end{center}

Some measurements are $\underline{\rm non-local}$, e.g., the correlation 
function
\begin{displaymath}
\Gamma (n) = \sum_i S_i S_{i+n}
\end{displaymath}

However, this can also be easily calculated in parallel.  We just take a copy 
of the spin variables and keep shifting them one site at a time.  After the 
first shift, everyone calculates their partial sum of $\Gamma (1)$, i.e., $S_i 
S_{i+1}$, then after the second shift everyone calculates $\Gamma (2)$, etc.  
After the spin has been shifted halfway around the lattice (only need to go 
halfway, since periodic boundary conditions $\Rightarrow \Gamma (n) = \Gamma 
(L-n)$), we do the sum over sites $\Gamma (n) = \mbox{SUM}\ \left( S_i S_{i+n} 
\right)$ for all $n$.

\begin{center}
{\bf Relation to DE Solvers}
\end{center}

Parallelizing the Metropolis algorithm for a two-dimensional spin model is 
very similar to parallelizing numerical methods for solving simple 
differential equations, such as Laplace's equation $\nabla^2 \phi = 0$.  This 
can be discretized onto a two-dimensional grid, with the update depending only 
on nearest neighbor points, e.g., an iterative scheme to solve this equation 
would replace $\phi_i$ by the average value of its four neighbors.  Here it is 
possible to update all sites at once (Jacobi algorithm).  However, get better 
convergence by using red/black update (like Gauss-Seidel).  Parllelize this 
just like coarse grain parallel Metroplis algorithm.

\begin{center}
{\bf Parallel Cluster Algorithms}
\end{center}

For the Swendsen-Wang cluster algorithm, the standard measurements can be done 
in parallel, just as for Metropolis.  Updating the spins can easily be done in 
parallel, as can setting up the bonds---all the sites just look at their 
nearest neighbors in the positive direction, and put a bond there with 
probability $\left(1-d^{-\beta}\right)\,\delta \left( s_i, s_j \right)$ for 
the Potts model.  Note that $\beta=\frac{J}{kT}$.

The only difficult problem in parallel is the $\underline{\rm connected\ 
comonent\ labeling}$, that is, assigning to all the sites a cluster label that 
is unique for each cluster.

\begin{center}
{\bf The Need for a Parallel Algorithm}
\end{center}

The reason we want to use cluster algorithms is that they offer such a large 
reduction in computational time over the traditional Metropolis algorithm.  
For example, for the Ising model on a $512^2$ lattice, the number of 
Metropolis iterations required to generate an independent configuratio is of 
order 1000, whereas for the cluster algorithms it is of order 10, that is, a 
100 fold improvement in speed.  Howver, we have seen that early attempts at 
running cluster algorithms on parallel supercomputers gave speed-ups only of 
order 10, whereas the Metropolis algorithm can easily have speed-ups of 
100--1000 on vector and parallel supercomputers, so th improvement of the 
cluster algorithm is lost unless we can find an efficient parallel algorithm.

\begin{center}
{\bf Regular and Local Vs. Irregular and Non-Local}
\end{center}

This type of problem occurs in many different guises in computational physics.  
In many applications, there are standard algorithms that are regular, simple, 
possibly local, and thus parallelize very efficiently.  However, it is this 
ver regularity and locality that makes them poor algorithms.  Complex physical 
problems and physical systems are usually irregular, non-local, and change 
rapidly with time, so simple regular and local algorithms tend to have 
problems, such as critical slowing down.  More complex algorithms, such as 
non-local, irregular cluster algorithms\@; multiscale and multigrid methods 
for DEs and spins models\@; adaptive, irregular grids for finite element 
calculations\@; hierarchical, adaptive N-body solvers\@; etc., all work much 
better, but are much more difficult to parallelize!

\begin{center}
{\bf Labeling or Growing a Single Cluster}
\end{center}

Let us first look at how sequential component labeling algorithms work.  We 
will start with the Wolff algorithm, which ``grows'' a single cluster.  
Amusingly enough, we can liken this to the growth of an ant colony, in the 
following manner:

Pick a site at random and place an ant on it.  The ant reproduces by placing a 
child on the four neighboring sites (in two dimensions) with probability $p$.  
This ``first generation'' of ants also looks at each of its four neighbors in 
turn, and if it is not already accupied, places a child of its own there with 
probability $p$.  The ant colony continues to grow unti the final generation 
produces no children.

\begin{center}
{\bf ``Ants-in-the Labyrinth'' (Depth-First Search)}
\end{center}

Clearly, we grow a Wolff cluster in exactly the same fashion.  Placing a baby 
ant on a new site corresponds to putting a bond between the two sites (with 
probability $p$), and the ant colony is the cluster.  This is straightforward 
to program---we merely need to keep track of which sites are in the cluster, 
and which are the newly labeled sites (baby ants) that need to be looped over 
in the next iteration.

This algorithm can also be used to do the connected component labeling for all 
sites in the Swendsen-Wang algorithm, given the configuration of bonds.  We 
just add an extra loop over all sites.  If a site has already been labeld 
(assigned to a cluster), we go to the next site.  If not, we increment the 
label counter and give it to our ant, which propagates it through the 
labyrinth of bonds connecting sites in the cluster.

\begin{center}
{\bf Fisher-Galler-Hoshen-Kopelmann (Breadth-First Search)}
\end{center}

First proposed by Fisher and Galler.  Implemented in a more efficient manner 
by $H and K$, for percolation studies.

Uses the idea of a $\underline{\rm rooted\ tree}$.  Each site has a pointer to 
another site in the cluster, and so on down a tree-like structure.  The site, 
which is the root of the tree, has the label for the cluster.  Sites find 
their labels by traversal of the tree. Prune the tree during the algorithm 
(make shortcuts!) to make the algorithm efficient.

\begin{center}
{\bf Parallel Ants}
\end{center}

If we want to parallelize the Wolff algorithm, we need to use a parallel 
``ants-in-the-labyrinth'' algorithm, since that is how the Wolff cluster is 
grown.  Unfortunately, this is a sequential process---each ant generation has 
to follow the previous one.  The only parallelism is $\underline{\rm within}$ 
a generation, that is, over the ``wavefront'' of the expanding cluster.  There 
is some parallelism, but it is very sparse, especially if the cluster being 
grown is very small (which is often the case for the Wolff algorithm).  So the 
Wolff algorithm is very inefficient on a vector or parallel machine.

Clearly, this is not what we want to do for a parallel Swendsen-Wang algorithm, 
either.

\begin{center}
{\bf Local Label Propagation (Self-Labeling Local Diffusion)}
\end{center}

The simplest and most obvious parallel component labeling algorithm proceeds 
as follows.  Suppose we have one site per processor.  Assign a unique label to 
each site (e.g., the processor number).  Every site looks at its neighbor in 
the positive $x$ direction.  If they have a bond between them, and the 
neighbor's cluster label is smaller than the label at the site in question, 
then it sets its label to be the same as its neighbor's.  If the neighbor's 
label is bigger, then it receives the smaller label.  This is repreated in the 
$y$ direction.

\begin{center}
{\bf SIMD Local Propagation}
\end{center}

Notice that this algorithm is purely SIMD, since everyone does the same thing 
at the same time.  Since it propagates the label only one lattice spacing at a 
step, the number of steps required to label the clusters is the maximum path 
length (along bonded sites) between any two sites in the same cluster.  For 
cluster algorithms this will be order $L$ (the length of the lattice), since 
at $T_c$ there is one large cluster.  This method will work better for 
labeling small connected components, and worse for labeling objects such as 
spirals.  This method is much too inefficient for cluster algorithms, 
basically because the parallel algorithm suffers from the same kind of 
``critical slowing down'' (trying to label a large cluster using local steps) 
as does the sequential Metropolis algorithm!  Need a 
$\underline{\rm non-local}$ parallel algorithm.

\begin{center}
{\bf MIMD Local Propagation}
\end{center}

On a coarse grain MIMD machine, this algorithm is much more effective.  Here, 
we can label the clusters on the sub-lattice on each processor using a fast, 
efficient sequential algorithm.  The local label propagation is now done only 
on the edges of the sub-lattices, to match up labels of clusters that cross 
processor boundaries.

Now we have good load balancing, and label of largest cluster is propagated 
across number of processors, not number of sites.  For an $L^2$ lattice on a 
$p^2$ processor array:
\begin{eqnarray*}
\mbox{communication} & \sim & p\left( \frac{L}{p} \right) = L \\
\mbox{calculation}   & \sim & \frac{L^2}{p^2} \\
\end{eqnarray*}

As long as $L>>p$, will get good efficiencies for this method.
\begin{displaymath}
\mbox{comm/calc} \sim \frac{p^2}{L} \hspace{1em} \mbox{\# proc} << L
\end{displaymath}

\begin{center}
{\bf Parallel Labeling Algorithms for Image Processing Applications}
\end{center}

There is a large literature of parallel algorithms for pixel connectivity and 
connected component labeling applications in image processing.  But many of 
these are useless for our application.
\begin{enumerate}
\item Pixel connectivity is a special case of component labeling, and some 
algorithms cannot be generalized.

\item Most computer science analsysis of algorithms concnetrates on 
$\underline{\rm worse\ case}$ complexity.  We are interested only in 
$\underline{\rm average}$ time complexity to label a large number of 
configurations.

\item Optimal asymptotic efficiencies are $\underline{\rm not}$ necessarily 
optimal on real machines.  May need $>> 10^6$ processors!

\item Algorithms may work well on small, regular objects, but not on large, 
irregular spin model clusters.
\end{enumerate}

\begin{center}
{\bf Global Equivalencing}
\end{center}

Cluster labeling can be looked on as a special case of an equivalence 
problem\@: given a list of equivalences (site $i \equiv$ site $j$), put all 
the elements into equivalence classes (cluster labels).  There are many 
algorithms known for solving this problem.

Parallel implementation:
\begin{enumerate}
\item Label clusters sequentially on each node.

\item Each processors looks at the edges of its neighbors and creates a list 
of equivalences.

\item Lists are sent to one node, which creates the global equivalence classes 
and returns the results to all nodes.
\end{enumerate}

\begin{center}
{\bf Hierarchical Method}
\end{center}

Problems with global equivalencing:
\begin{enumerate}
\item Equivalence step sequential bottleneck
\item One node needs enough memory for all cluster labels
\end{enumerate}

Solved (partially) by hierarchically merging the labels.

Same algorithm used by Embrechts, et al., for component labeling problem. Take 
$\log P$ steps rather than $P$.  Complicated to implement---could be improved.

\begin{center}
{\bf ``Multi-grid'' Algorithm (SIMD)}
\end{center}

Instead of just looking at neighboring sites, look at sites a distance $2^m$ 
away (this is fast on a hypercube).

If the label is the same at any point, put a connection between the sites.  If 
there are connections $A \stackrel{2^m}{\rightarrow} B 
\stackrel{2^m}{\rightarrow} C$, make a connection $A 
\stackrel{2^{m+1}}{\rightarrow} C$.

Now changes in the labels can be propagated long distances in one iteration.  
Expect labeling in $\sim \log_2 L$ iterations, each of which takes $\sim 
\log_2 L$ multi-grid steps, each taking time $o(1)$ on hypercube.  So time 
complexity of $O(\log L)^2$.

\begin{center}
{\bf Get/Send (Shiloach-Vishkin)}
\end{center}

Another non-local SIMD algorithm that we implemented was a method similar to 
an algorithm due to Shiloach and Vishkin.  It is based on the rooted tree 
approach of the FGHK algorithm (breadth-first search).  

Again, we start with every processor having a unique label (i.e., the root of 
its own tree).  Each iteration consists of the following steps:
\begin{enumerate}
\item Local label propagation

\item Send

If a site's label changes during step (1), it sends the new label to the 
processor labeled by this site's old label (i.e., the root of the current sub-
cluster), which picks the smallest of all the labels it receives as its new 
label.

\item Scan

Propagate labels fast along rows and columns.

\item Get

Every site gets the (possibly new) label from the root processor for its 
current label.  The get/send algorithm again has $\sim \log L$ iterations 
(rather than $L$), and $\nu (\log L)^2$ time complexity since non-local 
communications take time $O(\log L)$.

Get/send faster in practice than multigrid, both much faster than local label 
propagation.

Only implemented on SIMD machines, where performance is poor due to usual load 
balance problems.  Again, should be much better on MIMD machines (not 
implemented yet).
\end{enumerate}

\begin{center}
{\bf Independent Parallelism}
\end{center}

All is not lost for the parallelization of highly irregular Monte Carlo 
methods, such as the Wolff algorithm.  This is because domain decomposition is 
not the only level of parallelism available for MC algorithms.

Remember that the basic idea of MC simulation is that we are approximating an 
extremely large or infinite sum (or large or infinite dimensional integral) by 
a finite sum over important configurations. The error is predominantly 
statistical, that is, proportional to $1/\sqrt{N}$, where $N$ is the number of 
independent configurations.  So here is another way the simulation can be done 
in parallel---if different processors are generating independent 
configurations to add to the sum.

\begin{center}
{\bf Independent MC Simulations}
\end{center}

This is very easy to do.  We just set up different random number streams on 
each processor (which we would do in any case on a parallel machine, but here 
different ``processors'' could be different workstations in a network), set up 
different random initial configurations on each processor, and then run 
completely independent simulations on each processor.  So now we get a 
parallel efficiency of 100\%!!

This simple idea works very well for cluster algorithms, as long as the 
lattice is small enough that all the data will fit into the memory of a single 
processor, which may be a problem for parallel machines, but not too bad for 
workstation networks.  Note that this only works for MIMD machines, since 
cluster algorithm is irregular.

\begin{center}
{\bf ``Hybrid'' Parallelism}
\end{center}

We have seen that parallel S-W cluster algorithms for moderate sized lattices 
$(\sim 512^2)$ work fairly efficiently for moderate numbers of nodes $(\sim 
64)$ on a MIMD machine.  However, this does not prevent us from using 
massively parallel MIMD machines with $\sim 1000$ nodes for this problem.  We 
can use a hybrid of domain (data) parallelism with independent parallelism 
(over different MC runs) by, for example, running eight independent 
simulations of 64 nodes each on a 512-node machine.

We did this for the old nCUBE-1 at Caltech, and got performances equivalent to 
an eight-processor CRAY~Y-MP for this problem.

\begin{center}
{\bf More on Spin Models}
\end{center}

Recall that the Ising model has energy
\begin{displaymath}
E = -J \sum_{\lbrack i,j\rbrack} S_i S_j \hspace{1em} S = 
\stackrel{\pm}{\uparrow} \stackrel{1}{\downarrow}
\end{displaymath}

For $J>0$ we have a ferromagnet.  Energy is smallest when neighboring spins 
are aligned $\uparrow\uparrow$, $\downarrow\downarrow$, $+1 +1$, $-1 -1$.  
Ground state ($T=0$ lowest energy state) when all spins are aligned, or 
parallel.

For $J<0$ we have an antiferromagnet.  Energy is smallest when neighboring 
spins are opposite $\uparrow\downarrow$, $\downarrow\uparrow$, $+1 -1$, $-1 
+1$.  Ground state when all spins are antiparallel, that is, checkboard 
configuration for square two-dimensional lattice.  Only two possible ground 
states for Ising ferro- and antiferro-magnets on square lattice.

\begin{center}
{\bf Frustration}
\end{center}

Consider an Ising ferromagnet defined on a triangular lattice.  Similar to 
square lattice, except every site now has six nearest neighbors instead of 
four.  Triangular lattice ferromagnet has similar phase transition, same 
critical exponents, as square lattice ferromagnet.

This is $\underline{\rm not}$ the case for the triangular lattice 
antiferromagnet.  Behavior is very different, since not all neighboring sites 
can simultaneously have lowest energy.  Links are $\underline{\rm 
unsatisfied}$.  The mode is said to be $\underline{\rm frustrated}$.  Dealing 
with frustration is a difficult problem!

\begin{center}
{\bf Spin Glasses}
\end{center}

Generalize the interactions in an Ising spin model, so interaction strength 
becomes a variable, dependent on the different sites:
\begin{displaymath}
E = - \sum_{\langle i,j \rangle} J_{ij} S_i S_j
\end{displaymath}

Consider the model where the (magnitude of) strength is constant $\left( 
\left\vert J_{ij} \right\vert = J \right)$, but varying in sign, i.e.,
\begin{displaymath}
J_{ij} = +J \mbox{or}\ -J, \hspace{1em} J>0
\end{displaymath}

Suppose we choose the sign at random for every $i,j$ link.  This will 
introduce frustration, since not all links around certain plaquettes can be 
satisfied.

Randomness and frustration are hallmarks of real systems called spin glasses.

\begin{center}
{\bf Properties of Real Spin Glasses}
\end{center}

Certain metallic alloys have magnetic interactions which oscillate in value as 
a function of the separation of the different atoms.

These allows can be created in an $\underline{\rm amorphous}$ state---the 
positions of the atoms are random and disordered like a glass, rather than 
regular and ordered like a crystal.

The magnetic interactions between spins are thus sometimes ferromagnetic, 
sometimes antiferromagnetic.  This is what is known as a spin glass.

Spin glasses have a number of interesting properties:
\begin{itemize}
\item transition to a ``glassy'' phase at low temperatures,

\item susceptibility has a cusp rather than a divergence at $T_c$,

\item very large number of metastable states (local minimum of energy) in spin 
glass phase,

\item long equilibration times (logarithmic rather than exponential relaxation 
to equilibrium)---``remnant magnetication'' $M\sim 1/\log t$,

\item degeneracy (or near-degeneracy) of ground state, with many 
configurations having (nearly) the lowest energy value,

\item finding ground states is a difficult (NP hard) problem,

\item main properties of real spin glasses (e.g., amorphous magnets with 
oscillating interaction strength) can be described fairly well by $\pm J$ spin 
glass Ising model on a regular crystalling lattice,

\item to study spin glasses using computer simulation, need a method of 
producing low $T$ configurations without getting stuck in a higher energy 
metastable configuration.
\end{itemize}

\begin{center}
{\bf Monte Carlo Simulation of Spin Glasses}
\end{center}

The interesting behavior in a spin glass occurs at very low temperature, so we 
want to do a Monte Carlo simulation of an Ising spin glass at low $T$.

Suppose we adopt the usual MC approach---generate an initial random 
configuration of spins, and then start running the simulation at the 
temperature $T$ where we want to study the system.  After some thermalization 
time $\tau$, the system will reach thermal equilibriu, and we can start taking 
measurements as usual.  Right?  Well, not exactly.

\begin{center}
{\bf Quasi-Ergodicity}
\end{center}

Recall that a Monte Carlo algorithm will work only if it is ergodic, that is, 
it is possible to reach any state from any other state in a finite number of 
iterations.

The Metropolis algorithm is ergodic for $T>0$ since there is always a finite 
probability $e^{-\Delta E/kT}$ of flipping a spin.

However, the algorithm may not be ergodic at $T=0$, since moving from one 
ground state to another in a frustrated spin model may require higher energy 
intermediate states.  For $T$ very small, the probability of climbing out of a 
local minimum is non-zero, but may be so small that in any real simulation it 
will never happen.  The system is then only quasi-ergodic.

\begin{center}
{\bf Thermalization of Spin Glasses}
\end{center}

Thus, for a system like an Ising spin glass with a rough ``energy landscape'' 
with many local minima, starting from a random initial configuration and 
running at a very low temperature $T$, the system will invariably get trapped 
in a local minimum, and never reach a thermalized state (at least not in the 
finite time of a real simulation).

This is analogous to ``quenching'' a metal, that is, starting from a very hot 
(disordered) state, and rapidly cooling to a low temperature state.  For real 
metals, this will also produce a local minimum of the energy, leading to a 
brittle, non-crystalling state.

\begin{center}
{\bf Simulated Annealing}
\end{center}

In order to avoid the metastable states produced by quenching, metals are 
often cooled very slowly, which allows them time to order themselves into 
stable, structurally strong, low energy configurations.  This is called 
``annealing.''

Annealing gives the system the opportunity to jump out of local minima with a 
reasonable probability while the temperature is still relatively high.

We can adopt the same approach in Monte Carlo simulations.  We start with a 
random configuration at a very high temperature, and then reduce the 
temperature ``very slowly'' until we reach the desired low temperature.  This 
should result in a thermalized configuration.

\begin{center}
{\bf Finding Ground States}
\end{center}

In many cases, we are interested in finding the ground state of a spin model, 
that is, the zero temperature, lowest energy state.  In many regular models, 
it is possible to deduce the form of the ground state configurations.  
However, for disordered frustrated systems, such as spin glasses, finding a 
ground state is very difficult (in fact it is NP hard---i.e., there is no 
polynomial time algorithm for finding a solution).

In practice, ground states are found using simulated annealing.  This involves 
starting from a number of different random initial configurations, and slowly 
cooling them down to zero temperature.  This produces local energy minima 
that are close to, and possibly equal to, the global minimum.

\begin{center}
{\bf Optimizatino}
\end{center}

Finding the ground state of a spin glass is an example of an optimization 
problem:

Find the minimum (or maximum) value of a $\underline{\rm cost\ function}$ 
$C(s)$ for all possible states $s$.

For the spin glass, the cost function is the energy, and the states are all 
$2^N$ possible configurations of an $N$ site lattice.

Kirkpatrick, et al. realized that the method of simulated annealing that they 
were using for finding ground states in spin models could be generalized to 
solving $\underline{\rm arbitrary}$ optimizatin problems by associating the 
cost function with the energy, and introducing a fictitious temperature.

\begin{center}
{\bf Optimization Using Simulated Annealing}
\end{center}

For a general optimization problem, the temperature is just a parameter that 
governs the probability of increasing the cost function at any step, via the 
usual Metropolis algorithm form $e^{-\delta C/T}$, where $\Delta C$ is the 
change in the cost function due to a change in the configuration.  Just as for 
the spin glass, having a non-zero temperature allows the procedure to jump out 
of local minima.  Zero temperature corresponds to a ``steepest descent'' type 
algorithm, where only changes that do not increase the energy are accepted.

Simulated annealing works well for a number of combinatorial optimization 
problems.

\begin{center}
{\bf Finding an Optimal Solution}
\end{center}

We expect that as long as the temperature is reduced ``slowly enough,'' that 
is, in small enough increments, and allowing time for thermalization at each 
$T$, then simulated annealing should allow us to reach the ``ground state,'' 
that is, the optimal solution.  But what is ``slow enough''?

Geman and Geman showed that if the temperature is reduced as
\begin{displaymath}
T_k = \frac{T_0}{\log k}
\end{displaymath}
for $T_0$ large enough (so that the configurations at $T_0$ are close to 
random), then we are ``statistically'' guaranteed to find the optimal value.

\begin{center}
{\bf Statistically Optimal Solutions}
\end{center}

This result means that if we anneal slowly enough, we have a non-zero 
probability of obtaining the optimal solution.  Thus, if we perform enough 
different annealings, starting at different random initial configurations, we 
will eventually obtain an optimal solution.  Will generally give good (if not 
optimal) solution.  The trouble is, although the probability of finding a 
solution is non-zero, it may be incredibly small (c.f. quasi-ergodicity), so 
in practice we may never find the optimal solution.  And even if we do, we 
cannot know it is optimal without using some other analytic techniques (e.g., 
branch and bound).

\begin{center}
{\bf Advantages of Simulated Annealing}
\end{center}

Although this may seem like a rather weak statement, guaranteeing a 
statistically optimal solution for arbitrary cost functions is something no 
other optimization technique can claim.

Simulated annealing:
\begin{enumerate}
\item can deal with arbitrary systems and cost functions

\item statistically guarantees finding an optimal solution

\item is relatively easy to code

\item generally gives a ``good'' solution
\end{enumerate}

This makes S.A. an attractive option for optimization problems where 
``heuristic'' (specialized or problem specific) methods are not available.

\begin{center}
{\bf Disadvantages of Simulated Annealing}
\end{center}

\begin{itemize}
\item Too slow

Repeatedly annealing with a $1/\log k$ schedule is very slow, especially if 
the cost function is expensive to compute

\item May be ``overkill''

For problems where the energy landscape is smooth, or there are few local 
minima, SA is overkill---simpler, faster methods (e.g., gradient descent) will 
work better.  But usually don't know what the energy landscape is.

\item Heuristic methods can be superior

Heuristic methods, which are problem-specific or take advantage of extra 
information about the system, will often be better than general methods.  But 
SA is often comparable to heuristics.
\end{itemize}

\begin{center}
{\bf The Annealing Schedule}
\end{center}

To quantify how fast we perform the annealing, we need to specify:
\begin{description}
\item[$T_k$:] the range of temperaturesss
\item[$N_k$:] the number of Metropolis iterations at each $T_k$
\end{description}

This is known as the annealing schedule.  $T_k$ and $N_k$ are parameters that 
can be tuned to improve the performance of SA for different problems.  Finding 
a good annealing schedule is key to producing good solutions in a reasonable 
time.

Often don't need true optimal solution---just a ``good'' solution, say within 
1\% of the lowest value.  Tradeoff between quality of solution and time to 
find it.

\begin{center}
{\bf Simulated Quenching}
\end{center}

The Geman and Geman result is a bound on the annealing schedule in order to 
guarantee an optimal solution.  For some problems, we can do this with a 
faster annealing schedule, and in man cases when we don't need an optimal 
solution, a faster annealing may be more efficient.

A popular annealing schedule is one that is exponential rather than 
logarithmic:
\begin{eqnarray*}
T_{k+1} & = & CT_k \hspace{1em} 0<C<1 \\
T_k & = & \frac{T_0}{\exp\lbrace (1-C)k\rbrace} \\
\end{eqnarray*}
This is sometimes known as ``simulated quenching.''  Usually, gives quite good 
results as long as $C$ and $N_k$ are chosen appropriately.

\begin{center}
{\bf Determining the Schedule from Statistical Mechanics Ideas}
\end{center}

Ideally, one would like a method for determining the annealing schedule based 
on the behavior of the specific system. We can obtain some guiding principles 
by considering the statistical mechanics of the problem.  One good rule of 
thumb is that we would like to minimize the variation from equilibrium.  
During the annealing, we can measure the energy (cost) and specific heat 
(variance of the energy) at each temperature value.  These values can be used 
to help determine the annealing schedule.  Note that since we are constantly 
changing $T$, it is not easy to obtain equilibrium values of $E$ and $C_H$ for 
each temperature.

\begin{center}
{\bf A ``Physical'' Annealing Schedule}
\end{center}

\begin{description}
\item[$N_k$:] We would like to have $N_k$ be such that we reach thermal 
equilibrium at $T_k$ before moving to $T_{k+1}$.

Can do this using standard methods for testing for equilibration (binning).  
Or alternatively, can use clever techniques for estimating $\tau_{\exp}$ from 
measurements of $E$, and then set $N_k \approx \tau_{\exp}$.

\item[$T_k$:] To avoid moving too far from equilibrium, could choose $T_{k+1}$ 
so that fluctuations in $E\left( T_k \right)$ overlap those of $E\left( 
T_{k+1} \right)$.
\end{description}

\begin{eqnarray*}
C & = & \frac{dE}{dT} \approx \frac{\Delta E}{\Delta T} \\
T_{k+1} & \approx & T_k - \frac{\Delta E}{C_k} \\
\end{eqnarray*}

\begin{center}
{\bf Examples of Optimization Problems}
\end{center}

\begin{itemize}
\item Graph partitioning

In designing computer chips, we have a number of connected elements.  If there 
are too many elements to fit on one chip, they need to be partitioned onto 
(for example) two chips.  Since wires between chips are expensive, we want to 
minimize the number of connections between chips.  Since silicon area is also 
expensive, we want to balance the number of elements on each chip.

Introduce a variable $S_i$ which is $+1$ if the element is on the first chip, 
and $-1$ if it is on the second.  Introduce a connectivity matrix $C_{ij} = 1$ 
if $i$ and $j$ are connected, 0 otherwise.

The penalty for connections between chips is
\begin{displaymath}
H = -\sum_{ij} C_{ij} S_i S_j
\end{displaymath}

We want roughly the same number of elements in chip 1 $(S_i = +1)$ and chip 2 
$(S_i = -1)$, so we want to make the magnetization $\sum_i S_i$ as close as 
possible to zero.  Can do this by adding a penalty term to the cost function
\begin{eqnarray*}
H & = & - \sum_{ij} C_{ij} S_i S_j + \mu \left( \sum_i S_i \right)^2 \\
& = & \mbox{constant}\ - \sum_{ii} J_{ij} S_i S_j \\
\end{eqnarray*}
where $J_{ij} = C_{ij} - 2\mu$

This is the form of an Ising spin glass!  As we might expect, SA works well 
for this problem.
\end{itemize}

\begin{center}
{\bf The Traveling Salesman Problem}
\end{center}

The classic optimization problem can be very simply stated---a salesman has to 
visit $N$ cities, and wants to take the shortest possible route.  Here, the 
cost funtion is the length of the tour.

Here, the change to a ``configuration'' (a particular tour) is not as simple 
as a single spin flip in the spin glass problem, or moving an element from one 
side to the other in the graph partitioning problem.

Each tour can be presented as a permutation of the numbers 1 to $N$, which 
represent the cities.  The simplest change to the tour is to swap paris of 
cities, and measure the change in the tour path.

\begin{center}
{\bf The Update Step}
\end{center}

Here, the update is performed by swapping the order of two adjacent cities in 
the tour (B and C).  There is no reason why we have to swap adjacent cities in 
the tour---larger, non-local swaps may be more effective at moving through 
configuration space, especially at high $T$.  At low $T$, large moves are not 
so useful since they will mainly be rejected.  The type of move can depend on 
the temperature.  Choosing efficient, ergodic **THE FAX DIDN'T PICK UP THIS 
NEXT LINE** .... for SA.

\begin{center}
{\bf Parallel Simulated Annealing}
\end{center}

SA can be easily parallelized on a coarse grain machine by using independent 
parallelism, that is, using a different random initial condition and different 
random number streams on each processor.  Each processor will find an 
independent local minma, and we choose the smallest of these as our best 
solution.  Alternatively, we may want to pass information between processors 
during the annealing.  For example, replicate ``good'' configurations.

For large problems, we may be happy with a good approximate solution, so we 
just parallelize the problem. For example, domain decomposition for spin 
glass, map tour onto a ring of processors **FAX DIDN'T PICK UP**

\begin{center}
{\bf Simulated Tempering}
\end{center}

Simulated tempering was invented by Parisi and Marinari to study disordered 
Ising spin models (spins in a random external magnetic field).  These models 
have a first order phase transition, so there are two coexisting states at 
$T_c$---a high energy and low energy state.  Standard Metropolis algorithm 
tends to get stuck in one of the these states.

Tempering does a Monte Carlo update of the $\underline{\rm temperature}$---try 
to change $T$, do a Metropolis accept/reject depending on $\Delta (E/kT)$.  
Sweep slightly above and below $T_c$, moves you in an out of high and low $E$ 
states, so do correct sampling of configurations at $T_c$.  **FAX DIDN'T PICK 
UP**

\begin{center}
{\bf Tempering for Optimization}
\end{center}

One of the main problems with simulated annealing is that it is difficult to 
find an optimal annealing schedule, and that changing $T$ drives the system 
out of equilibrium.

If we use tempering instead, the system is always in equilibrium.  Also, we 
can determine the temperature changes (the equivalent of the annealing 
schedule) automatically, by varying the $\Delta T$s so that the acceptance of 
a change is $\sim 50$\%.  $\Delta T$s may be too small, random walk in $T$ may 
not be optimal---need to check this out.

\begin{center}
{\bf Lattice Gauge Theory}
\end{center}

Lattice gauge theory, and in particular lattice QCD, is probably the 
application area with the largest use of (vector and) parallel supercomputers, 
with the possible exception of CFD and weather forecasting.

It was the driving force behind the construction of a number of major centers 
for parallel computing applications (Caltech, Los Alamos, Edinburgh, 
Southampton, ...), and the construction of a number of powerful parallel 
computers (Caltech, Columbia, APE, Fermilab, ...).  One of the DOE ``Grand 
Challenge'' areas.  Currently driving the ``TeraFLOPS Initiative,''' to build 
a 1~TLOPS machine by 1995. **NOT SURE ABOUT LAST THREE WORDS - FAX CUTOFF**

\begin{center}
{\bf Background to QCD}
\end{center}

Probably the major achievements in physics this century have been the 
development of relativity and quantum mechanics, and the synthesis of these 
two theories in the 50's into quantum field theories of fundamental particles 
and forces.

\begin{description}
\item[1950's:] Quantum electrodynamics (QED), quantum field theory describing 
electromagnetism (photons, electrons, positrons).

\item[1960's:] Quantum chromodynamics (QCD), quantum field theory the strong 
nuclear (``color'') force (quarks, gluons, hadrons).

\item[1970's:] Unified theory of the weak and electromagnetic forces (VWS) 
(neutrinos, W. Z. Higgs).
\end{description}

\begin{center}
{\bf Lattice QCD}
\end{center}

QED is a fantastically successful theory.  Some properties of electrons can be 
measured, and calculated using QED, to precisions of one part in $10^8$, and 
agree within errors.  Calculations can be done analytically using perturbation 
theory (series expansion in a small parameter---usually the strength of the 
force or interaction).

QCD is similarly a very impressive theory that has helped make sense of the 
bewildering array of subatomic particles discovered in colliders (``atom 
smashers'').  Unfortunately, since the interaction strength is so much larger 
than for QED, perturbation theory does not work well---not at all in some 
cases, maybe to one part in 10 in others.

Since QCD cannot be solved analytically, people turned to approximate methods, 
simpler semi-empirical models, etc., to do non-perturbative calculations.  
Wilson developed a technique for approximating QCD and other quantum field 
theories by discretizing them onto a lattice in a way that preserved the 
fundamental properties and symmetries of the theory, and also provided a 
simple way of dealing with divergences due to self-interactions, which are a 
problem with QFTs at very short distances.

Discretizing the theory onto a finite lattice also made it amenable to 
computer simulation.

\begin{center}
{\bf Lattice QCD and Parallel Computing}
\end{center}

The first numerical studies of lattice QCD were done around 1979.  Early work 
used standard mainframe computers of the time, such as the VAX.  It soon 
became apparent that these machines could only be used for very preliminary 
studies, and that many orders of magnitude more computing power was required 
to get useful results to compare with experiment.

In the late 1970s, vector supercomputers were becoming available, and it was 
realized that they were ideal for lattice QCD, which is very regular, and very 
numerically intensive.

In the early 1980s, parallel computers were being introduced, and it was 
realized that these would do just as ***FAX CUT OFF LAST LINE***

\begin{center}
{\bf QCD and Parallel Computing}
\end{center}

Lattice field theory, and in particular lattice QCD, are ideal candidates for 
using parallel computing.
\begin{itemize}
\item Extremely numerically intensive.  Estimated that $\sim 1$~TLOPS year of 
calculation required for definitive results from full theory.

\item Regular lattice structure (SIMD).

\item Calculation/communication ratio is very large.

\item Operations are (mostly) local.

\item Scalable---use more processors to do larger problems.

\item Simulations take so long that researchers are willing to spend many 
months programming and optimizing their codes in difficult software nad 
hardware environments. ***FAX CUT OFF LAST LINE***
\end{itemize}

\begin{center}
{\bf Experimental QCD and Parallel Computing}
\end{center}

As an aside, it is worth noting that not only people working on the 
theoretical QCD are using parallel computers.

Experimentalists, who study particle interactions in large colliders, need to 
model these interactions (which are probabilistic quantum processes) using 
Monte Carlo simultion, for which (we have seen) parallel computers are ideal.

Measurements of interactions are written to tape for later analysis.  since 
there are far too many ``events'' to store them all, need a real-time system 
to rapidly check each event, and only save ``interesting'' ones.  this is done 
using ``processor farms'' of hundreds of processors.  The SSC will need 
thousands of processors and Tbytes ***FAX CUT LAST LINE***

\begin{center}
{\bf Quantum Field Theory}
\end{center}

In classical (Newtonian) mechanics, the motion of a body is strictly 
determined by Newton's laws of motion.  These laws can be derived from a very 
simple idea called the ``action principle''---particles move so as to minimize 
the $\underline{\rm action}$ $S$, which is related to the energy of the body.  
The pat taken by a particle between two points A and B is the path of least 
action. 

In quantum mechanics, everything is probabilistic, and all paths 
$\underline{x}(t)$ are possible.  Feynman showed that the probability of a 
particle going from A to B can be calculated from a ``path integral'' (or 
functional integral)
\begin{displaymath}
K(A,B) = \sum_{{\rm path}s} e^{i\,s({\rm path})}
\end{displaymath}

***THE REST OF THE EQUATION CUT OFF BY FAX***

\begin{center}
{\bf Field Theory and Statistical Mechanics}
\end{center}

In QFT, the fundamental objects are not paths of particles, but quantum fields 
$\Phi (x)$ describing the particles.  $x=\left(\underline{x},t\right)$ is a 
position in space-time.

Probabilities are calculated using a ``generating functional''
\begin{displaymath}
Z = \int \lbrack D\Phi \rbrack e^{i s\lbrack\Phi\rbrack}
\end{displaymath}
It is possible to make a transformation Minkowski space-time to regular 
Euclidean four-dimensional space.  This gives:
\begin{displaymath}
Z_E = \int_{\rm fields} \lbrack D\Phi \rbrack e^{-s_E\lbrack\Phi\rbrack}
\end{displaymath}
Note that this is formally identical to the partition function in statistical 
mechanics:
\begin{displaymath}
Z = \int_{\rm configs} \lbrack DC \rbrack e^{-H\lbrack C\rbrack/kT}
\end{displaymath}
***THE REST OF THE EQUATION WAS CUT BY FAX***

The fundamental quantity in statistical mechanics is the correlation function
\begin{displaymath}
\langle S(x) S(0)\rangle \sim e^{-x/\xi}
\end{displaymath}
which expresses correlations in the system (between spins in a spin model) in 
terms of the correlation length $\xi$.

The fundamental quantity in field theory is the ``propagator'' $\langle \Phi 
\left( x_1 \right) \Phi \left( x_2 \right) \rangle$, which describes the 
motion of the particle between $x_1$ and $x_2$ (c.f. $K(A,B)$ in quantum 
mechanics).  If we look at the propagation of the particle in time, we have 
that
\begin{displaymath}
\langle \Phi (t) \Phi (0)\rangle \sim e^{-mt} \hspace{1em} \Phi (t) = \int d^3 
\underline{x} \Phi \left( \underline{x}, t \right)
\end{displaymath}
where $m$ is the ``mass gap,'' the **LAST LINE CUT OFF BY FAX**

There is a deep relationship between quantum field theory and statistical 
mechanics, which shows itself in the mathematical similarities between 
partition function and generating functional, correlation function and 
propagator, correlation length and inverse mass.

What this means is that when we discretize field theories onto a lattice, we 
can apply all the well-known numerical techniques (Monte Carlo simulation, 
importance sampling, finite size scaling) and theoretical ideas (theory of 
phase transitions and critical phenomena) from statistical mechanics to the 
seemingly completely different area of quantum field theory.


\end{document}

