
% BB Nov 22 2008.

\documentclass[20pt,landscape,footrule]{foils}
\usepackage[dvips]{graphicx}
\usepackage[pdftex]{geometry}
\geometry{headsep=3ex,hscale=0.9}
\usepackage[pdftex]{color} % for the colored block
\date{\color{blue} Singapore MCMC Meeting, March 2004}
\usepackage[ignore]{pause}
%\usepackage{pause}
\usepackage{background}
%\usepackage{pp4slide}
%\MyLogo{AV HEP,CSIT}
\LogoOff
\begin{document}
\def\r{\color{red}}
\def\b{\color{blue}}
\def\s{\color{black}}

\foilhead{\r Parallel Computing} \b

After briefly discussing the often neglected, but in praxis frequently
encountered, issue of {\r trivially parallel} computing, we turn to 
parallel computing with information exchange. Our illustration is the 
{\r replica exchange method}, also called {\r parallel tempering.} 
Closely related (but less good) is the Jump-Walker (J-Walker) method. 

Parallel tempering has a low, but non-trivial, communication overhead 
between the processes. Information about the temperatures, energies 
and nearest neighbor locations of $\s n$ systems is communicated, and
logical decisions have to be made. This algorithm gives some kind of 
ideal start into the world of parallel processing and is also of major 
practical importance. In some applications it competes with the 
multicanonical algorithm, in others both algorithms supplement one 
another, and combinations of them have been explored.

\foilhead{\r Trivially parallel computing} \b

An example of a trivially parallel large scale endeavor is the {\it 
Folding@Home} distributed computing project, which achieved sustained 
speeds of about thirty teraflops in molecular dynamics simulations.

There are many tasks, which simply can be done in parallel, for instance 
{\r canonical MC simulations at different temperatures and/or distinct 
lattices.} For such applications the gain due to having many (identical) 
processors available scales linearly with the number of processors. In 
practice trivially parallel applications are frequently encountered. 
Another example are simulations of {\r spin glass replica.} It is 
far more cost efficient to run trivially parallel applications on PC 
clusters than on dedicated parallel machines (supercomputers), which 
feature expensive communication networks. % But, it appears always easier 
% to get CPU time funding on supercomputers than to find large-scale 
% trivially parallel clusters.

For MCMC simulations of single systems, which last (very) long in real 
time, parallelization is trivial, but there is some overhead due to 
multiple equilibration runs. Typically the {\r CPU time} consumption 
of such a simulation is divided into 
$$\s t_{tot}= t_{eq} + t_{meas}\ . $$
For parallelization of such a computation on $\s n$ independent nodes 
the single processor CPU time scales as
$$\s t_{tot}^1 = t_{eq} + t_{meas} / n\ .$$
{\r How many processors} should one choose (assuming a sufficiently large
supply)? A good rule of thumb is
$$\s n = {\tt nint}\,\left[ t_{meas}\over t_{eq} \right]\ .$$
The total processor time used becomes
$\s t_{tot}^n = n\,t_{eq}+t_{meas} \approx 2\,t_{meas}$, which limits
the CPU time loss due to multiple equilibrations to less than two
$$\s R_{loss} = {t^n_{tot}\over t_{tot} } \approx
  {2\,t_{meas} \over t_{eq}+t_{meas} } < 2\ , $$
and the gain is a reduction of the run time by the factor
$$\s R_{gain}^{-1} = {t^1_{tot}\over t_{tot} } \approx
  {2\,t_{eq} \over t_{eq}+ n\,t_{eq} } = {2\over 1+n}\ . $$
The worst case is $\s n=2$, for which the improvement factor in CPU 
time is $\s 2/3$. Obviously this parallelization makes most sense for 
$\s t_{meas}\gg t_{eq}$. 

\foilhead{\r Replica exchange method (parallel tempering)} \b

The method of multiple Markov chains was introduced by Geyer and, 
independently, by Hukusima and Nemoto under the name {\bf replica 
exchange} method. The latter work was influenced by the {\bf simulated 
tempering}, which coined the name {\bf parallel tempering} for an 
exchange of temperatures in multiple Markov chains. But not that
simulated tempering is distinct method and actually a special case 
of the method of {\bf expanded ensembles}. 
A precursor of such ideas is a 1986 PRL by Swendsen and Wang, where 
a general replica exchange method is formulated (apparently its 
generality prevented one from identifying the important special 
case).

In practice one does not {\r exchange} configuration but {\r 
temperatures or similar parameters,} so that the amount of exchanged 
information is small. Therefore, the method is particularly well 
suited for parallel processing on clusters of workstations, which 
lack the expensive fast communication hardware of the dedicated 
parallel computers. The Jump Walker (J-Walker) approach does not
exchange configurations, but feeds replica from a higher temperature 
into a simulation at a lower temperature, and has consequently
problems with balance.

Parallel tempering performs $\s n$ canonical MC simulations
at different $\s\beta$-values with Boltzmann weight factors
$$\s w_{B,i} (E^{(k)}) = e^{-\beta_i E^{(k)}} = e^{-H}\,,\
     i=0, \dots , n-1\ , $$
where $\s \beta_0 < \beta_1 < ... < \beta_{n-2} < \beta_{n-1},$ 
and allows to exchange neighboring $\s\beta$-values:
$$\s \beta_{i-1} \longleftrightarrow \beta_i\ ~~{\rm for}~~
     i=1, \dots , n-1\ . $$
Their joint weight is
$$\s  e^{-\beta_{i-1}E_{i-1}-\beta_iE_i} =  e^{-H} $$
and the $\s \beta_{i-1}\leftrightarrow\beta_i$ exchange leads to 
\newpage

\begin{eqnarray} \nonumber \s
- \triangle H &\s = &\s \left( -\beta_{i-1} E^{(k)}_i -   
\beta_i E^{(k')}_{i-1} \right) - \left( -\beta_i E^{(k)}_i
  -   \beta_{i-1} E^{(k')}_{i-1} \right) \\ \nonumber
&\s = &\s\ \left( \beta_i-\beta_{i-1} \right)\,
      \left( E^{(k)}_i - E^{(k')}_{i-1} \right)
\end{eqnarray}
which is \r accepted or rejected according to the Metropolis algorithm,
\b {\it i.e.}, with probability one for $\s \triangle H\le 0$ and with 
probability $\s \exp (-\triangle H)$ for $\s \triangle H>0.$

The CPU time spent on the $\s \beta$ exchange should be less than 50\% of 
the total updating time. In practice it is normally much less (a few 
percent). The $\s \beta_i$ spacing should to be determined so that a 
reasonably large acceptance rate is obtained for {\bf each} pair. This 
can be done by a recursion, which is a slight modification of one by 
Kerler and Rehberg. 

Let us denote by $\s a^{pt}_i$ the {\bf acceptance rate of the $\s 
\beta_{i-1}\leftrightarrow\beta_i$ exchange}. Assume $\s a^{pt}_i>0$ 
for all $\s i=1,\dots,n-1$.  Iteration $\s m$ is performed with the 
$\s \beta$ values \s $\beta^m_i$, $\s i=0,\dots,n-1$. \b We define the 
$\s \beta^{m+1}_i$ values of the next iteration by
$$\s \beta_0^{m+1} = \beta_0^m ~~~{\rm and}~~~ \beta_i^{m+1} = 
     \beta_{i-1}^{m+1} + a_i^m\, ( \beta_i^m - \beta_{i-1}^m ) 
     ~~{\rm for}~~ i=1,\dots,n-1\ , $$
where
$$\s a_i^m\ =\ \lambda^m\, a^{pt,m}_i ~~~{\rm with}~~~
     \lambda^m = { \beta_{n-1}^m - \beta_0^m \over
     \sum_{i=1}^{n-1} a^{pt,m}_i\, (\beta_i^m - \beta_{i-1}^m) }\ .$$
For large acceptance rates $\s a^{pt}_i$ the distance $\s \beta_i^m
-\beta_{i-1}^m$ is increased, whereas it shrinks for small acceptance 
rates. The definition of the $\s a_i$ coefficients guarantees 
$\s \beta_{n-1}^{m+1}=\beta_{n-1}^m$, {\it i.e.} {\r the initial 
$\s \beta_{n-1}-\beta_0$ distance is kept.}

A problem of the recursion is still that the statistical noise stays 
similar in each iteration step. It may be preferable to combine 
all the obtained estimates with suitable weight factors $\s w^k.$ 
The final estimate after  $\s m$ iterations reads then
$$\s \beta_i = {\sum_{k=1}^mw^k\beta^k_i\over \sum_{k=1}^mw^k}\ .$$
Two reasonable choices for the weights are 
$$\s w^k = \min_i\, \{ a_i^{pt,k} \}\ , $$
which determines the weight from the worst exchange acceptance rate and 
suppresses the statistics of runs with a low $\s w^k$, and
$$\s w^k = {1\over \sqrt{\sigma^2}} ~~{\rm with}~~ \sigma^2 =
     \sum_{i=1}^n \left( {1\over a^{pt,k}_i} \right)^2\ , $$
which relies on what would be the correct statistical weighting of the 
acceptance rates, if they were statistically independent (what they 
are not).

\foilhead{\r Computer implementation} \b

We like to program the $\s \beta$ exchange using MPI. Initially the
indices of the $\s\beta$ values are chosen to agree with the rank
{\s\tt MY\_ID} of the process. The exchanges lead then to the situation 
for which the $\s \beta$ indices become an arbitrary permutation of the 
numbers $\s 0,\dots,n-1$. Two approaches are possible:

\begin{enumerate}

\item All the $\s n$ processes send the needed information (the 
energies of their configurations) to a master processor. The master 
handles the exchange and sends each process its new $\s \beta$ value
back (which can agree with the old value).

\item The involved processes handle the exchanges through point to 
point communication.

\end{enumerate}

We consider only the second approach, because for it the number of
sweeps per time unit has a chance to scale linearly with the number
of processors, whereas for the first approach the master processor 
will become the bottleneck of the simulation when many processors are
employed (when only few processors are used, it may be more efficient 
to let a master processor handle the exchanges).

Each process has to know the identities of the processes where the 
neighboring $\s\beta$ values reside. They are stored in a {\r neighbor 
array} {\s\tt NEIGH(2)}. The initial values are those defined in the 
{\s\tt phbwrite\_mpi.f} test program of {\s\tt 1MPI.} With each 
$\s\beta$ exchange the nearest neighbor and next nearest neighbor 
connections get modified and a {\r correct bookkeeping} is the
main difficulty of implementing the parallel tempering algorithm.
We cannot deal with all pairs simultaneously, but need one buffer
process between each pair. This leads to three distinct interaction
groups, starting either with \s $\beta_0$, $\beta_1$ or $\beta_2$, 
\b which we label by $\s {\tt INDEX} = 0,\,1,\,2$ as indicated for 
eight processes in the following figure ({\it e.g.}, for $\s 
{\tt INDEX}=0$ the buffer processes are {\s $i=2$ and~$i=5$)}: 

\newpage

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/parallel_t.pdf}
% \caption{\small .}
\end{center} \end{figure}

Our MPI Fortran subroutine {\s\tt p\_pt\_mpi.f} implements the $\s 
\beta$ exchange (as in the figure the variable {\s\tt INDEX} takes 
on the values {\s 0, 1 or 2}).

After the lowest $\s\beta$ value is determined, the higher $\s \beta$ 
values follow according to the scheme of the figure. The thus
identified processes send their {\s\tt iact}, {\s\tt IACPT\_PT} (for 
$\s a^{pt}_i$) and {\s\tt NEIGH(1)} values to their {\s\tt NEIGH(2)} 
neighbors, {\it i.e.} to the processes which accommodate their upper, 
neighboring $\s\beta$ values. These neighbor processes decide whether 
the $\s\beta$ exchange is accepted or rejected, relying on the 
logical-valued function 
$$\s {\tt LPT\_ EX\_ IA(BETA1,BETA2,IACT1,IACT2)}\ ,$$
which implements the Metropolis criterion. For either 
outcome the processes send two integers back to their lower neighbors 
(the originally initiating processes) and one integer to their upper
neighbors. If the exchange is {\r rejected,} the first element of each 
{\s\tt MPI\_SEND} is set to $\s -2,$ thus signaling to the receiver
that nothing changed. If {\r accepted,} the receivers store the 
addresses of their new neighbors in the {\s\tt NEIGH} array. In 
addition, two more actions are performed: First, the initiating 
processes inform their left neighbors about the change by sending 
a new neighbor address to {\s\tt NDEST3R}. This send is performed 
in any case, where the message $\s -2$ informs {\s \tt NDEST3R} that 
nothing changed. Second, the processes whose $\s \beta$ values 
changed re-calculate their heatbath weights and interchange also 
their parallel tempering acceptance rates {\s\tt IACPT\_PT}.

After a certain number of $\beta$ exchanges, the $\s\beta$ values are
updated by either of the subroutines
$$\s {\tt pt\_ rec0\_ mpi.f},\ {\tt pt\_ rec1\_ mpi.f}~~{\rm or}~~
     {\tt pt\_ rec2\_ mpi.f}\ . $$
In practice it may be unavoidable that one (or several) of the 
acceptance rates $\s a^{pt}_i$ are zero. Applying the Kerler-Rehberg 
recursion blindly gives the undesirable result $\s\beta_{i-1}'=\beta_i.$ 
To avoid this, we
$$\s {\rm replace\ all}~~ a^{pt}_i=0~~ {\rm values\ by} 
     ~~a^{pt}_{\min} = 0.5/N_{\rm update} $$
where $\s N_{\rm update}$ is the number of updates done. This works
when there are enough $\s \beta_i$ values to bridge the 
$\s [\beta_1,\beta_{n-1}]$ interval. Otherwise, the only solution 
to the problem of $\s a^{pt}_i=0$ values is to increase the number 
of processes.

We could gather the needed permutation of $\s \beta$ values, as well 
as their corresponding exchange acceptance rates, using point to 
point MPI communication. However, it is shorter to rely on the
\begin{eqnarray} \s\nonumber &\s {\tt MPI\_ ALLGATHER}&
  \s \!\!\! {\tt (SENDBUF,ISENDCOUNT,SENDTYPE,RECVBUF,} \\
  & &\s {\tt IRECVCOUNT,RECVTYPE,COMM,IERROR)} \nonumber
\end{eqnarray} 
instruction. The arguments of {\s\tt MPI\_ALLGATHER} are: {\s\tt 
SENDBUF,} the starting element of the send buffer; {\s\tt ISENDCOUNT} 
(integer), the number of elements in the send buffer; {\s\tt SENDTYPE,} 
the data type of the send buffer elements; {\s\tt RECVBUF,} the first 
element of the receive buffer; {\s\tt IRECVCOUNT} (integer), the number 
of elements in the receive buffer; {\s\tt RECVTYPE}, the data type of 
the receive buffer; {\s\tt COMM}, the communicator; {\s\tt IERROR,} 
the MPI error code.

When a very large number of processes is involved, one may want to
reprogram this in form of point to point interactions.
\newpage

Finally a usefull instruction, in particular for debugging purposes, is
$$\s {\tt MPI\_BCAST(BUF,ICOUNT,DATATYPE,ISOURCE,COMM,IERROR)} $$
which broadcasts the buffer {\s\tt BUF} from the process with rank 
{\s\tt ISOURCE} to all processes specified by the communicator {\s\tt 
COMM}. The number of elements in the buffer is given by {\s\tt ICOUNT,} 
the datatype by {\s\tt DATATYPE}, and {\s\tt IERROR} is the MPI error 
code.

\end{document}

