% 
% 1993 npac reu program
% Final paper
%
\documentstyle[leqno,twoside,11pt,ltexproc]{article}
\setlength{\oddsidemargin}{+1in}
\setlength{\evensidemargin}{+.6in}
\setlength{\topmargin}{+5mm}
\raggedbottom

% This file allows the insertion of postscript files into 
% the document
%
\include{psfig}
\begin{document}
\cleardoublepage
\pagestyle{myheadings}

% Old title: Global Optimization With Simulated Annealing and
%	     Simulated Tempering 

\title{\bf A New Automatic Simulated Annealing--type Global
	   Optimization Scheme 
\thanks{
 Project conducted as part of the 1993 Research Experiences for 
 Undergraduates (REU) Program
 in High-Performance Computing conducted by the Northeast Parallel
 Architectures Center (NPAC) at Syracuse University.  The National
 Science Foundation provides primary funding for the 1993 NPAC REU
 Program, through NSF Grant CDA-9200577.  Additional funding is provided by
 the GE Foundation Faculty for the Future Program; the
 Vice President for Research and Computing, the Dean of the
 College of Engineering and Computer Science, the Dean of the
 College of Arts and Sciences, and NPAC, Syracuse University.}}

\author{M. Shore 	
	\thanks{Research Apprentice, 1993 NPAC REU Program;
	Majors: Computer Science, Music, Florida State University.}
\and P. D. Coddington
	\thanks{Project Leader, NPAC, Syracuse, NY.}
\and G. C. Fox 
	\thanks{Director, NPAC, Syracuse, NY; Professor of Physics,
	Syracuse University, Syracuse, NY.}
\and E. Marinari
	\thanks{Associate Director, NPAC, Syracuse, NY; Professor of
	Physics, Syracuse University, Syracuse, NY; Professor of
	Physics, Universit\`{a} di Roma {\em Tor Vergata}, Via E.
	Cernevale, 00173 Roma, Italy.}
\and C. H. Wu
	\thanks{Graduate Assistant, CIS, Syracuse University,
	Syracuse, NY.} }

\date{}
\maketitle

% The following places author's name next to page number on even pages and
% title on odd pages.
\markboth{Shore \& Coddington}
	 {A New Automatic Simulated Annealing--type Global
	  Optimization Scheme} 


\begin{abstract}
 
A study of global optimization schemes is presented.  Simulated
annealing, a general and proven technique for minimizing functions
with many coexisting states, is used as a foundation for the
development of a new, more automatic approach, called simulated
tempering.  This novel method upholds the eminent attribute of
simulated annealing -- the probabilistic guarantee of convergence upon
a global minimum.  It is unique, however, in that system equilibrium
is never disturbed.  Although simulated tempering is in its infancy,
it is promising, indeed, as the preliminary results suggest.  The
two--dimensional Ising Spin Glass Model, an NP--complete optimization
problem, is used as the test case.  Its theoretical formulation is
briefly addressed.

\end{abstract}


\section{Introduction}

Global optimization, also known as multivariate and combinatorial
optimization, is a technique for minimizing complex functions
dependent on many variables.  Such problems are prevalent in science
and, with the advent of high--performance computing, have attracted
considerable attention \cite{Ingber} \cite{Otten}.  Intuitively, the
formulation consists of a function to be minimized, otherwise known as
the {\em cost function}, and a number of constraints, which specify
the feasibility of options.  Consider, for example, the task of
finding the minimum value of a function which defines a surface is
3--space.  The value(s) of the function at the critical points suggest
possible optima of which one is convinced after applying the second
derivative test (luckily, there were no saddle points!)  Armed with
two or three discrete results, it is clear which one is smallest (or
largest).  This method may be characterized as {\em efficient}, i.e.
with it, one is able to solve all {\em instances}, or {\em
dimensions}, of the problem in a computation time bounded by a
polynomial the size of the instance.  (In the example above, the
instance would equal 3.)

In theory, analytical techniques could be generalized to {\em
N}--space.  In practice, however, there exits a class of optimization
problems for which analytical schemes, like the above, are not
realistic.  Take, for instance, a macroscopic crystal.  {\em N} in
this case would be on the order of $10^{23}$ -- even contemplation of
a numeric calculation is absurd!  Indeed, there are no known efficient
algorithms for optimizing many hard problems.  (The Travelling
Salesman Problem and variations of the Ising spin glass model are
classic examples.)  This class of problems is called {\em
NP--complete} (non--polynomial time complete) and is the motivation
for this work.

{\em Simulated annealing} is a popular and mature approach to
combinatorial optimization and forms the foundation for this project.
It is roughly analogous to the physical process of annealing, where a
melted substance is gradually cooled until a crystalline state is
achieved.  In a simulation, the ``frozen'' substance corresponds to
the minimal {\em cost}, or {\em ground state}, of the system.  (The
term {\em system} is used loosely, as this methodology is, in fact, a
general one.  For example, it can refer to a collection of cities in
the Travelling Salesman Problem or to a group of units in complex
circuit design.)  Simulated annealing is an iterative method which
works due to a simple premise: If the cost obtained by evaluating a
configuration of a system at temperature {\em t} is less than that of
temperature {\em t+1}, then the system is updated to that ``less
expensive'' state.  If, on the other hand, the cost of the proposed
configuration is higher, it is accepted probabilistically, i.e. it may
or may not be accepted.  This process has a spectacular implication.
By accepting ``upward'' moves some of the time, this algorithm is able
to ``jump'' out of local minima in its search!  In fact, simulated
annealing guarantees an optimal solution \cite{Geman}, if truly
gradual temperature decrements are observed.

As one may imagine, several ingredients are required to develop a
general--purpose global optimization technique.  At present, methods
that employ simulated annealing need to address a number difficult
questions.  What is an optimal rate for decreasing the temperature?
How long must the system remain at the current temperature before a
move to a lower one can be made?  What is a suitable criteria for
convergence?  Clearly, a working implementation would envelope the
simple simulated annealing algorithm in an array of complex,
computationally--intensive routines.  Indeed, this presents the main
obstacle for using this technique.

Recently, {\em simulated tempering}, a novel simulated annealing based
approach was introduced which possesses promising attributes
\cite{Marinari}, \cite{Wu}.  Similar to annealing, the tempering
scheme draws on an analogy with a physical process.  Instead of
gradual cooling however, the desired texture of some substance is
achieved by repeated cooling and heating.  Theoretically, this
technique performs several crucial operations automatically, that
require sophisticated implementation for simulated annealing , all the
while preserving the salient feature of this class of optimization
methods -- the probabilistic guarantee of convergence upon the global
minimum.

This work may be thought of as a preliminary step towards the
development of a general multivariate optimization scheme.  The study
is organized naturally in two phases: an exploration of simulated
annealing followed by the investigation of simulated tempering.  The
structure of this paper is dictated by the same format.  After a brief
presentation of the theoretical background for the experimental
setting, the two--dimensional Ising spin glass, simulated annealing
and simulated tempering are covered.  Results and observations are
discussed last.



\section{The Ising Spin Glass Model}

The Ising model, named after, Ernst Ising \footnote{This model was
actually originated by Ising's supervisor, Lenz.  However, it was
Ising who fist published results on the model and hence his name was
adopted.}, is certainly a monumental contribution to physics.  It aims
to describe {\em phase transitions}, phenomena which occur when a
small change in a parameter such as temperature causes a large--scale
change in the state of the system.  Although a thorough discussion of
the Ising model is well outside the scope of this work, it is
nevertheless important to describe the experimental setting for this
study.  In addition, the formulation of the Ising model conveniently
addresses several issues of statistical mechanics relevant to
simulated annealing--type optimization schemes.

In brief, the model consists of an arrangement of binary values, or
{\em spins}, on a lattice.  Every lattice site may correspond to an
electron, for instance, with a tiny magnetic field.  An independent
variable, $\varsigma_i$, is assigned to each lattice site, {\em i =
1,\dots, N}, which collectively represent the configuration space, {\em
X}, of the system.  An essential ingredient in the Ising model is the
sum over all possible configurations.  Since there are $2^{N}$
possible configurations, this sum is clearly enormous if {\em N} is at
all large.

The total energy of the system is described by the {\em Hamiltonian} 
(a natural choice for the cost function).  Only nearest--neighbor
interactions and interactions of lattice sites with an external field
contribute to this quantity.  A formal definition is in order.  For
each configuration, $\varsigma = (\varsigma_1,\dots, \varsigma_N)$

\begin{equation}
  H = H(\varsigma) = 
   - \sum_{<i,j>} E \varsigma_i \varsigma_j - \sum_i J \varsigma_i,
\end{equation}

\noindent
where {\em E} and {\em J} are the energies associated with
interactions with nearest--neighbors and with the external field,
respectively.  

The {\em partition function}, the central object in statistical
mechanics, also plays a fundamental role for the Ising model.  From
it, one may derive all of the important thermodynamic features of the
modeled physical system.  

\begin{equation}
  Z = Z(\beta, E, J, N) = \sum_{\varsigma = \pm 1} 
	\exp(-\beta H(\varsigma)). 
\end{equation}

\noindent
Here any dimensions the Hamiltonian may have are cancelled due to the
$\beta$ parameter.  Typically, $\beta = 1/kT$, where k is Boltzmann's
constant and {\em T} is temperature.  Since the temperature, {\em T},
is incorporated into $\beta$, they will be used interchangeably
throughout the remainder of this paper.  The reader may wish to
consult \cite{Cipra}, an excellent introduction to the Ising model.


\section{Simulated Annealing}

As was alluded to in the introduction, a general simulated annealing
implementation requires several finely tuned parameters.  Perhaps two
are most crucial: the {\em decorrelation} time, which refers to the
amount of time (i.e. the number of iterations) the system is kept at a
given temperature for sufficient ``thermalization,'' and the {\em
annealing schedule}, or the temperature scale, {\em T}, which
corresponds to the rate at which the system is ``cooled.''  Though the
latter has attracted considerable attention, the former is seldom
rigorously addressed.  The sections following the formulation of the
simulated annealing method describe an adaptive annealing schedule.
Due to the constraints of time, the decorrelation criterion was not
investigated.


\subsection{First formulation}

One of the early successful applications of the computer was the
simulation of systems with many degrees of freedom.  This was
motivated by experimental condensed matter and high energy physicists
who wished to study the interactions of many particles.
Unfortunately, to find the equations of state of models consisting of
perhaps millions of elements is not feasible on even today's
high--performance supercomputers.  The science of statistical
mechanics addresses this problem by providing heuristics for accurate
estimation of the behavior of {\em microsystems}, systems consisting
of a few elements, based on {\em macrosystems}, systems consisting of
many.

In a typical scenario, a model may be comprised of a number of spheres
in a box which interact by forces that decrease with distance.  The
energy content of the system could be calculated from the distances
between the particles.  A simulation may proceed as follows.  A new
configuration is constructed from the current one by imposing a random
displacement to a randomly selected particle at each step of the
simulation.  If the energy of this new state is lower than that of the
current state the displacement is accepted unconditionally and the
system is updated to the new configuration.  If the energy of the new
state is higher by $\Delta\varepsilon$ it is accepted with the
probability

\begin{equation}
  \exp \left( \frac{-\Delta\varepsilon}{kT} \right).
  \label{eq:prob}
\end{equation}

This is the fundamental procedure of simulated annealing and is known
as the Metropolis step.  It was shown \cite{Metropolis} \cite{Otten}
that many iterations of the Metropolis method of state generation lead
to a distribution of states in which the probability of finding the
system in microstate {\em i} with energy $\varepsilon_i$ is given by

\begin{equation}
  P(\varepsilon_i) = 
  \frac{\exp(-\varepsilon_i/kT)}{\Sigma_j\exp(-\varepsilon_j/kT)}. 
  \label{eq:boltz}
\end{equation}

\noindent
(Notice that the denominator is the partition function.)  

The probability function described by equation (~\ref{eq:boltz}) is
known as the Boltzmann density function and has several interesting
implications.  At very high temperatures, for instance, each state has
almost equal chances of being the current state.  Consider an
experiment which operates on a sample of size {\em N} from the space
of configurations, {\em X}.  It is clear that as $T \rightarrow
\infty$, the probability of finding the system in any of {\em N}
states approaches $\frac{1}{N}$.  For low temperatures, however, only
those configurations which result in low energies have a high
probability of being selected, as should be clear by similar
reasoning.


\subsection{The annealing schedule}

Since the formulation of simulated annealing as a general--purpose
global optimization scheme by Kirkpatrick {\em et al.}
\cite{Kirkpatrick} and Cern\'{y} \cite{Cerny}, the annealing
schedule has been the subject of much contemplation \cite{Andresen}
\cite {Hajek} \cite{Hoffmann} \cite{Nulton} \cite{Rees}
\cite{Salamon}.   This is not surprising since the rate at which a
system is cooled is not merely a performance consideration, but also
one of correctness of the algorithm.  It is illustrative to draw
further on the physical process of annealing.  Consider, for instance,
a melted substance which is gradually cooled and thus driven toward a
crystalline state.  If the temperature is decreased too quickly, the
resulting crystal may have many defects or even lack all crystalline
order.

Most approaches to dynamically adjusting the annealing schedule
utilize a concept known as {\em constant thermodynamic speed}.
Intuitively, a condition is imposed that forces the level of
statistical uncertainty between two contiguous configurations to be
constant.  In our experimental implementation statistical uncertainty
was obtained by calculating standard deviations of an ensemble of
energies associated with configurations of a ``thermalized'' system
(i.e. one which is ready to move to a lower temperature).  Nulton and
Salamon \cite{Nulton} offer an excellent discussion of this idea.

Following is a summary of the way the said adaptive annealing schedule
works.

\begin{figure}[htb]
\begin{enumerate}
  \item Thermalize the system at temperature, {\em T = t}.
  \item Calculate the standard deviation of an ensemble of energies
	resulting at {\em T = t}, {\em sd(E(t))}.
  \item Calculate the subsequent temperature, {\em T = t + 1}, by
\begin{equation}
  T(t+1) = T(t) \exp \left( \frac{-\lambda T(t)}{sd(E(t))} \right).
  \label{eq:adapt}
\end{equation}
\end{enumerate}
\caption{The algorithm for the adaptive annealing schedule.}
\end{figure}

\noindent
The $\lambda$ and the number of energies in the ensemble are constants
currently determined {\em a priori}.  (A feasible method for obtaining
these values dynamically is proposed in section 5.3.)


\section{Simulated Tempering}

It should be clear that the fundamental problem with the simulated
annealing method is that the system is driven out of equilibrium with
every temperature decrement.  Though sophisticated techniques for
generating temperature schedules and decorrelation times have been
investigated, the algorithm suffers in two ways.  First, its
implementation becomes significantly more complex, and second,
additional computational time is incurred.  Simulated tempering, a new
method developed by Marinari and Parisi \cite{Marinari}, overcomes
these difficulties more elegantly and proficiently.


\subsection{Formulation}

The development of the simulated tempering algorithm was motivated by
models with phase transitions of order one or greater.  Such systems
are characterized by coexisting states at $T_c$ (i.e.  ``freezing'').
The original experimental setting, for instance, involved the Random
Field Ising Model which has a first order phase transition, i.e. two
simultaneous states at $T_c$: a high and low energy state \cite{Wu}.
Such problems present some difficulty for the annealing approach due
to a strictly decreasing temperature scale (the algorithm will tend to
get stuck in {\em one} of these states).  The sampling scheme of
simulated tempering, however, is not constrained by a decreasing
temperature schedule.  A simulation of tempering, in fact, treats the
temperature parameter, {\em T}, stochastically.  In actuality, the
temperature selection procedure is quite similar to the Metropolis
method of phase space evolution.

The tempering procedure results in important advantages over
annealing.  Namely, the system {\em always} remains in equilibrium and
a broader exploration of local minima is possible.  It is realized in
the following way.  A discrete variable, {\em m, (m = 1,\dots, M)}, is
added to the configuration space, {\em X}.  Intuitively, one may
regard this step as an introduction of another stochastic parameter
which represents a thermal quantity.  This allows for a Monte Carlo
procedure to be performed for temperature selection, and thus
represents a deviation from one fundamental standard simulated
annealing premise -- a continual decrease in temperature.  The
probability distribution is given by

\begin{equation}
   P(X,m) \propto \exp (-H_t(X,m)), 
\end{equation}

\noindent
where $\beta$ is incorporated into the variable, $H_t$ ($H_t$ is not
actually the Hamiltonian because of the {\em m} parameter.  In the
context of this work, $H_t$ will stand for a new term, the tempering
Hamiltoninan), resulting from the new configuration space,

\begin{equation}
  H_t(X,m) \equiv \beta_m H(X) - g_m.
  \label{eq:gm}
\end{equation}

\noindent
The $g_m$ parameter will be considered in more detail in section 5.2.
For the present, it is sufficient to state that it is a normalizing
energy value that is required to make equation (~\ref{eq:gm}) true due
to the introduction of {\em m} into the new configuration space.

It is clear that this formulation of $H_t$ results in the Boltzmann
distribution for $\beta_m = \beta$ (i.e. in the case of simple phase
space, {\em X}).  In the relevant instance of discrete $\beta$'s,
however, the probability of selecting a particular one is described by

\begin{equation}
  P(m) \equiv \exp(-(\beta_m f_m - g_m)),
 \label{eq:gmprob}
\end{equation}

\noindent
where $f_m$ represents the corresponding free energy.

Indeed, equation (~\ref{eq:gmprob}) is proportional to 

\begin{equation}
  Z(\beta_m) \exp(g_m).
  \label{eq:prop}
\end{equation}

\noindent
Since the aim of the tempering method is to freely sample
configurations resulting from various temperatures, a careful
estimation of $g_m$ will prove to be fruitful.  In fact, all of the
{\em p(m)'s} will equal $\frac{1}{M}$ if $g_m = -\ln(Z_m)$, an
equality obtained by trivial manipulation of equation
(~\ref{eq:prop}).  This process incorporates temperature
into the stochastic procedure of phase space evolution.

\begin{figure}[htb]
\begin{enumerate}
  \item Thermalize system at $\beta_1$.
  \item Let $\beta_m$ be the current temperature.  Change to
	neighboring $\beta$'s in the following way.
  \begin{itemize}
	\item Generate random number, $\rho$, ~~~~ $0 < \rho < 1$.
	\item If $\rho < \frac{1}{2}$ attempt to move to $\beta_{m-1}$;
	otherwise attempt to move to $\beta_{m+1}$.
	\item {\em Accept} the new temperature with probability
	$\exp(-\Delta H_t)$.
  \end{itemize}
  \item Repeat step 2.
\end{enumerate}
\label{fig:acc}
\caption{Temperature sampling by simulated tempering.}
\end{figure}

In the algorithm described in Figure 2, $\beta$ acceptance simply
means that the system is ready to move to a new temperature.  This is
a subtle point that merits comment, as the tempering scheme requires
that the last action of the second step of the algorithm be taken some
of the time.  At present, the aim is to develop an implementation such
that acceptance be approximately equal for increasing and decreasing
temperatures.  In addition, the {\em acceptance rate}, or a percentage
of accepted temperatures over time, for increasing and decreasing
temperatures is believed to be optimal when it is around 50 \%.  (The
reader may wish to consult
\cite{Marinari} for the initial formulation of the simulated tempering
method.)


\section{Results and Discussion}

As is so often that case, time has presented a hindrance.  Indeed, a
major ingredient to the annealing--type approaches to multivariate
optimization, the decorrelation criterion, has yet to be investigated.
For instance, it would be rather gratuitous to consider performance
(e.g.  timings and iterations) when the current annealing code uses a
conservative estimate for the decorrelation (thermalization) time and
the tempering implementation ignores the issue all together and has
not yet been tested with the dynamic temperature schedule.  In
addition, since the primary effort revolved around the method and not
the problem, details of the results concerning the discovered ground
states of the spin glass models will not be considered in this text.
Suffice it to say that both the simulated annealing and simulated
tempering experiments consistently converged on minimum states with
energies of around -362.0 and -1430.0 for the 16 $\times$ 16 and 32
$\times$ 32 size models, respectively.  At any rate, this study
presents several truly enlightening discoveries.

\vspace{0.5in}
\begin{center}
\begin{tabular}{| c || c | c | c || c | c |} \hline
Problem & Thermalization  &  Total   &  Metropolis & Initial & Final  \\ 
size    & sweeps          &  sweeps  &  sweeps     & $\beta$ & $\beta$\\
\hline \hline
16 x 16 & 500             & 700      & 10          & 0.30    & 1.0180 \\ 
32 x 32 & 500             & 1000     & 50          & 0.20    & 1.5498 \\ 
\hline
\end{tabular} \\
\vspace{0.2in}
Table 1. {\em Runtime parameters for the annealing and tempering
simulations \footnote{The final $\beta$ values fluctuated slightly due
to variations in the temperature schedules.}.} 

\vspace{0.3in}
\begin{tabular}{| c | c |} 
\hline
$\lambda$ & Ensemble size \\ 
\hline \hline
10.0      & 50.0 \\ 
\hline
\end{tabular} \\
\vspace{0.2in}
Table 2. {\em Runtime parameters for the adaptive annealing schedule
described in section 3.2}
\end{center}


\subsection{A fast optimal annealing schedule?}

In order to converge upon a global minimum with probabilistic
certainty a logarithmic cooling schedule must be observed
\cite{Ingber}.  In practice, however, researchers seldom adhere to
this fundamental requirement because of the incurred computational
cost and commonly use exponential schedules, e.g.,

\begin{equation}
	T_{k+1} = \alpha T_k, ~~~ 0 < \alpha < 1, 
	\label{eq:exp_sched}
\end{equation}

\noindent
justified only by expediency.  For many applications the true global
optimum is not actually necessary -- oftentimes one {\em close} to the
global optimum is sufficient.  Despite this observation, since
annealing is a relatively delicate process, one can imagine how this
implementation aspect may result in inappropriate results.

This brings about the first major assertion of this work: the dynamic
annealing schedule outlined in section 3.2 is, in the context of the
two--dimensional Ising spin glass, exponential in nature, as Figure
~\ref{fig:sched_comp} depicts.  This finding clearly merits rigorous
investigation, as its possible ramifications are dramatic.  Indeed, it
is likely that one will be able to apply the {\em fast} simulated
annealing algorithm, i.e. one utilizing an exponential temperature
schedule, {\em and} be guaranteed an optimal result.

\begin{figure}[htb]
   \centerline{
        \psfig{figure=sched_comp_32.ps,width=3.0in} }
        \caption{Comparison of dynamic and exponential annealing schedules}
        \label{fig:sched_comp} 
\end{figure}

It is important to note, however, that by no means is the dynamic
cooling schedule superfluous.  There are two primary reasons for this
fact.  First, annealing schedules vary with problem size, and second,
the temperatures selected by the constant thermodynamic speed
annealing schedule may not resemble an exponential curve for other
hard optimization problems.  For the present intents and purposes and
in light of the fact that generating temperatures dynamically is a CPU
intensive process, an efficient implementation may utilize the dynamic
scheme only for the initial pass.  The $\alpha$ parameter in equation
(~\ref{eq:exp_sched}) for generating the exponential schedule can then
be extrapolated and used for generating an exponential cooling scale
in subsequent runs.


\subsection{Implementation considerations for simulated tempering}

The primary question for the tempering scheme is how to achieve the
kind of temperature acceptance described in section 4.1.  Since it is
the change in the Hamiltonian, $\Delta H$, that governs temperature
acceptance, a reasonably accurate method is needed to ``guess'' where
a suitable next temperature will lie.  This is truly a sensitive
matter, as if the guess is erroneous enough, temperature acceptance
will be greatly affected.  Consider an example where a move is
attempted from $\beta_m$ to $\beta_{m+1}$.  The change in the $H_t$ is
calculated in the following way.

\begin{equation}
  \Delta H_t = E \Delta \beta_m - \Delta g_m
  \label{eq:dh}
\end{equation}

\noindent
where {\em E} is the instantaneous energy value, {\em H(X)}, $\Delta
\beta_m = \beta_{m+1} - \beta_m$ and $\Delta g_m = g_{m+1} - g_m$.
The reader may recall from section 4.1 that new $\beta$'s are accepted
with probability $e^{-\Delta H_t}$.  Due to the fact the the quantity
$\Delta H_t$ in (~\ref{eq:dh}) is directly dependent on the value of
$\Delta g_m$, careful approximations are a must.

Theoretically, the criterion for new $\beta$ acceptance is quite
simple.  In fact, all that is required is that $\Delta g_m$ be an
intermediate energy value \cite{Marinari}, say $E(g_m)$, such that

\begin{equation}
  E(\beta_m) < E(g_m) < E(\beta_{m+1}).
  \label{eq:dg_req}
\end{equation}

\noindent


The condition described in equation (~\ref{eq:dg_req}) is realized in
a straightforward fashion: $E(g_m)$ is determined from the average of
the energies at $\beta_m$ and $\beta_{m+1}$ (or $\beta_m$ and
$\beta_{m-1}$, as the case may be).  From this formulation one is able
to increase the accuracy of the estimation by subdividing the interval
between $\beta$'s \cite{Wu}.  A general equation is given by

\begin{equation}
   \Delta g_m = E(g_m) = \frac{\Delta \beta_m}{n} 
   \left\{ \sum_{i=1}^{n} (E_{i,m}) - \frac{1}{2} (E_{1,m} + E_{n,m})
   \right\},
\end{equation}

\noindent
where {\em n} is the number of subdivisions.

The transition probabilities for visiting five $\beta$ values are
illustrated in Figure ~\ref{fig:relax}.  The reader may recall from
the formulation of the simulated tempering algorithm that two aspects
are desirable: similar acceptance values for increasing and
decreasing temperatures, and an overall acceptance rate of 50 \%.  In
that light, considerable improvement is seen in the {\em n = 5} case.

\begin{figure}[htb]
  \centerline{
     \psfig{figure=n1.ps,width=2.5in} 
     \psfig{figure=n5.ps,width=2.5in} }
     \caption{Transition probabilities for a 32 $\times$ 32 Ising spin
	glass model.}
     \label{fig:relax}
\end{figure}

The current tempering implementation further improves the $\Delta g_m$
estimate.  Since temperatures are revisited, previous $\Delta g$'s are
averaged in with the current one.  This procedure presents an
additional subtle advantage.  Since the Metropolis method for energy
calculation is probabilistic, the possibility of bogus results is real
and quite inevitable.  The new method of weight assignment to the
$\Delta g$ parameter reduces the ill--effects of the procedure over the
course of the simulation.  It is important to add that these results
were obtained from a new implementation of the tempering algorithm.
In addition to systematically updating the $\Delta g$ parameter, this
experimental version generated the temperature schedule (and thus the
$\Delta g$'s) internally.  Prior implementations required a
temperature list and a corresponding $\Delta g$ list produced by
simulated annealing runs.

\begin{figure}[htb]
\begin{enumerate}
  \item Thermalize the system at $\beta_0$.
  \item Generate a new $\beta$.
  \item Perform the tempering step as in 2 of Figure 2.
  \begin{itemize}
    \item If attempt to move to $\beta_{m+1}$ not previously visited,
	generate a new $\beta$.
  \end{itemize}
  \item Repeat step 3.
\end{enumerate}
\label{fig:impl}
\caption{The new simulated tempering implementation.}
\end{figure}

As Figure 5 points out, the actual simulation of tempering only occurs
on pre--existing temperatures.  In fact, new temperatures are only
added when the algorithm attempts to move to a {\em lower} temperature
not previously visited.  Furthermore, the current method assumes that
$\beta_0 = 0$ (i.e. $t_0 = \infty$) and thus requires a downward move
in the beginning of the simulation.  This also implies, of course,
that a move {\em below} $\beta_0$ would never take place.


\subsection{Simulated tempering and the temperature schedule}

To this point, the discussion considered theoretical and
implementation details.  It is consequential to recall the underlying
motivation, however: the development of an {\em efficient} and {\em
effective} global optimization method.  How can simulated tempering
contribute to this end?

Consider the graphs in Figure ~\ref{fig:acc_comp}.  The trend depicted
by the transition probabilities at different $\alpha$ values is
self--evident.  Indeed, this heuristic provides feasible automation.
Currently, the constants in the formulation of both the dynamic and
the exponential temperature schedules are not only named {\em a
priori}, but are initially determined quite crudely (i.e. guessed).
The $\beta$ acceptance rates obtained through the transition
probabilities from tempering simulations can be easily utilized to
adjust the scheduling parameters at runtime.  In the case of an
exponential temperature schedule, for example, this would entail the
manipulation of the $\alpha$ variable until the desired performance is
achieved.

\begin{figure}[htb]
  \centerline{ 
    \psfig{figure=tp_0.85.ps,width=2.3in}
    \psfig{figure=tp_0.89.ps,width=2.3in}
    \psfig{figure=tp_0.92.ps,width=2.3in} } 
  \caption{Acceptance comparison for different $\alpha$ values for
	a 32 $\times$ 32 model.  n, the subdivision value, is 5.  The 
	$\alpha = 0.89$ case is an almost ideal realization of the 
	requirements described in section 4.1.}
  \label{fig:acc_comp}
\end{figure}


\subsection{Overlap}

In sampling ground states of systems with perhaps millions of local
minima, how can one be certain that the algorithm is not visiting the
{\em same} (or only a few) local minima throughout the simulation?
Indeed, this is a legitimate and interesting question for all
applications of combinatorial optimization.  Fortunately, there is an
elegant and straightforward solution.

The aim here is to derive a quantitative way to observe how a system
evolves.  Intuitively, this is achieved by superimposing two
configurations and noting their similarities and differences.
Consider two local state configurations, $\omega^\alpha$ and
$\omega^\beta$, occurring at times {\em t} and {\em t+1},
respectively.  Let {\em V} be the total volume of the system (with the
Ising spin glass model, the volume is simply the number of lattice
sites).  The {\em overlap} \cite{Bachas}, $\omega$, is then given by

\begin{equation}
  Q_{\alpha \beta} = \frac{1}{V} \sum_{i=1}^{V} 
	\omega_i^\alpha \omega_i^\beta, 
\end{equation}

\noindent
It is clear that if $\omega^\alpha$ and $\omega^\beta$ are exactly the
same, $Q_{\alpha \beta}$ would equal 1.  Conversely, if the
configurations were exactly opposite, --1 would be the result.  

\begin{figure}[htb]
  \centerline{ 
    \psfig{figure=overlap.ps,width=3.5in} }
  \caption{A plot of overlap as a function of $\beta$.  As a point of
	interest, this data was obtained from a simulation of
	annealing.  In a simulation of tempering one would expect a 
	more undulating graph.}
  \label{fig:overlap}
\end{figure}

A module for calculating this quantity has been implemented and tested
in the context of simulated annealing and has truly provided an
abundance of information.  At high temperatures, for example, one
might expect {\em Q} to fluctuate about 0.  As the system freezes,
however, the overlap calculation would approach 1.  Figure
~\ref{fig:overlap} illustrates this phenomenon.  (To see why this is
so, recall from section 3.1 that the probability of incurring changes
to the configuration space at low temperatures is small.)  In
addition, the notion of overlap can serve as a robust dynamic
criterion for the termination of the algorithm.  Given a proper
cooling schedule and sufficient decorrelation times, one can conclude
with confidence that the system is in a ground state if the overlap
calculation from contiguous configurations approaches 1 to within some
small value, $\epsilon$.  ($\epsilon$ would, of course, have to be
determined experimentally.)



\section{Conclusions}

This work has brought to light several important contributions to the
study of multivariate optimization.  First, it has shown
experimentally that the exponential annealing schedule is optimal when
applied to the two--dimensional Ising spin glass model.  Next, results
from an improved implementation of the simulated tempering algorithm
were presented.  Satisfactory temperature transition probabilities
clearly indicate that the algorithm is functioning properly.  A
captious point that has not yet been overtly stated is the following.
Previous attempts of simulating tempering on larger Ising spin glass
models have been problematic (mainly due to poor estimation of the
$\Delta g$ parameter) \cite {Wu}.  The current implementation,
however, has performed equally well on the 32 $\times$ 32 size lattice
as on the 16 $\times$ 16 and has shown trends to that suggest
comparable performance for even larger problems.  (Larger system sizes
were not investigated because the current tempering implementation is
lacking the autocorrelation function.)  Finally, the usefulness of
simulated tempering was illustrated with regard to the determination
of an optimal cooling schedule.

Although these findings further confirm the conjecture that simulated
tempering is a truly viable method for combinatorial optimization,
there is, as was previously mentioned, crucial aspects that have still
to be investigated.  The realm of the Ising spin glass model, for
instance, calls for additional research: the simulation of even larger
models and the determination of an {\em actual} ground state are two
avenues which are already under consideration.  For the latter,
parallelization and the decorrelation gauge are obvious requisites.
Lastly, before a statement about the generality of the said scheme can
be made, the method needs to be applied to other hard optimization
problems, such as the Travelling Salesman Problem.


% The '*' suppresses section numbering.
\section*{Acknowledgements}

The primary author wishes to express his gratitude to two especially
influential people: Professor Edward Bogucz and T. J. Willis.  Dr.
Bogucz is the embodiment of the most noble ideals of academia (a realm
plagued with bureaucracy and politics).  T.J., in addition to many
intriguing discussions with regard to this work, has been my
late--night compadre and, perhaps more importantly, has shown me the
that there is only one way to drink coffee -- {\em black}.


\begin{thebibliography}{99}

\bibitem{Andresen}
Andresen, B., Hoffmann, K. H., Mosegaard, K., Hulton, J., Pedersen,
J., M. and Salamon, P., {\em On lumped models for thermodynamic
properties of simulated annealing problems}, J. Phys., vol. 49,
1485--1492, 1988.

\bibitem{Bachas}
Bachas, C. P., {\em Spin--glass model of crystal surfaces}, Phys. Rev.
Lett., vol. 54, no. 1, 53-55, 1985.

\bibitem{Cerny}
Cern\'{y}, V., {\em A thermodynamical approach to the travelling
salesman problem: An efficient simulation algorithm}, Report, Comenius
University, Bratislava, Czechoslovakia.

\bibitem{Cipra}
Cipra, B. A., {\em An introduction ot the Ising Model}, Math. Oper.
Res., 937-959, 1987.

\bibitem{Geman}
Geman, S. and Geman, D., {\em Stochastic Relaxation, Gibbs Distribution
and the Bayesian Restoration in images}, IEEE Trans. Patt. Anal. Mac.
Int., vol. 6, 721--741, 1984.

\bibitem{Hajek}
Hajek, B. {\em Cooling schedules for optimal annealing}, Mathematics
of Operations Research, vol. 13, no. 2, 311--329, 1988.

\bibitem{Hoffmann}
Hoffmann, K. H. and Salamon, P., {\em The optimal simulated annealing
schedule for a simple model}, J. Phys. A: Math. Gen., vol. 23,
3511--3523, 1990.

\bibitem{Ingber}
Ingber, L. {\em Simulated annealing: Practice versus theory},
Statistics and Computing, 1993, to appear.

\bibitem{Kirkpatrick}
Kirkpatrick, S., Gellat, C. D., Jr. and Vecchi, M. P., {\em Optimization
by simulated annealing}, Science, vol. 220, 671--680, 1983.

\bibitem{Marinari}
Marinari, E. and Parisi, G. {\em Simulated tempering: A new Monte
Carlo scheme}, Europhys. Lett., vol. 19, 451, 1992.

\bibitem{Metropolis}
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, E. {\em
Equation of state   calculation be fast computing machines}, J. Chem.
Phys., vol. 21, 1087--1092, 1953.

\bibitem{Mirkin}
Mirkin, G., Vasudevan, K. and Cook, F. A., {\em A comparison of
several cooling schedules for simulated annealing implemented on a
residual statistics problem}, Geophys. Res. Lett., vol. 20, no. 1,
77--80, 1993.

\bibitem{Nulton}
Nulton, J. D. and Salamon, P., {\em Statistical mechanics of
combinatorial optimization}, Phys. Rev. A, vol. 37, no. 4, 1351--1356,
1988. 

\bibitem{Otten}
Otten, R. H. J. M., van Ginneken, L. P. P. P., {\em The Annealing
Algorithm}, Kluwer Academic Publishers, Dordrecht, 1989.

\bibitem{Rees}
Rees, S. and Ball, R. C., {\em Criteria for an optimum simulated
annealing schedule for problems of the travelling salesman type}, J.
Phys. A: Math. Gen., vol. 20, 1239--1249, 1987.

\bibitem{Salamon}
Salamon, P., Nulton, J. D., Harland, J. R., Pedersen, J., Ruppeiner,
G. and Liao, L., {\em Simulated annealing with constant thermodynamic
speed}, Comp. Phys. Comm., vol. 49, 423--428, 1988.

\bibitem{Wu}
Wu, C. H., {\em On a study of simulated tempering as an optimization
method}, Technical Report, NPAC, Syracuse University, 1993.

\end{thebibliography}

\end{document}



