Control and Cybernetics Adaptive simulated annealing (ASA)

22 Pages·2013·4.74 MB·English

Control and Cybernetics
vol. 25 (1996) No. 1
Adaptive simulated annealing (ASA): lessons learned
by
Lester Ingber
Lester Ingber Research,
P.O.B. 857, McLean, VA 22101, U.S.A.,
e-mail: [email protected]
Abstract: Adaptive sirnulated annealing (ASA) is a global opti
mi~ation algorithm based ou an associated proof that the parameter
space can be ::mmplecl much more dfic:iently than by using other
previous simulated annealing algorithms. The author's ASA code
has been publicly available for over two years. During this time the
author has vohmteercd to help people via e-mail, and the feedback
obtained has been used to further develop the code. Some lessons
karued, in particular some which arc relevant to other simulated
annealing algorithms, arc described.
1. Introduction
Adaptive silllulatcd auncaling (ASA) is a global optimi~ation algorithm that
relics on randomly irnportauce-sarnpling the parameter space, i.e., in coutrast
to utili~ing deterministic approaches often used by OR and mathematical pro
gramming people.
Siuce the public release of the very fast simulated rcaunealing (VFSR) code
(Ingber, 1!)8!)), now called ASA (Ingber, 1!)!)3a), I have volunteered to help users
via e-mail ou code-specific problems they might enco1mter. The popularity of
the code, roughly measnrcd by the si~e of the ASAJist of people specifically
rerptestiug Lo he placed on a monthly list of updates-- now at about 500 names·
incrcascd about a factor of two after an article in The Wall StTeel; Journal
(Wofsey, 1!)!)3). I have had contact with thousands of users via e-mail, a few
via other clwmrels of cornrnuuication, <t!l(l it is dear that these arc just a small
fraction of people that have at least tried to use the code. Therefore, I conclude
tlw.t the code is quite stable and free of rnajor bugs.
This paper deals with some of the lessons learned from this interaction. I
thiuk that at least some of these lessons rnay be useful to other developers of
sirnnlated mmealiug (SA) code, as well as to many users. A previous paper
has described SA in the context of its usc across many disciplines, and the
34
L. INGBER
variation iu use in actual practice versus the algorithm offered in theoretical
papers (Ingber, 1993b).
Section 2 gives a brief mathematical description of the ASA algorithm. Sec
tiou 3 describes some of the options in the code, and explains why they are
useful. Section 4 addresses some comparison tests and controversial claims,
emphasizing that in fact there are many shades of SA. Section 5 gives a brief
conclusion.
ASA source code in C-language is publicly available via anonymous ftp (Ing
ber, 1993a), or via the world wide web (WWW) at
http:/ /www.alumni.caltech.edu/rvingberf Problems with the code can be ad
dressed to [email protected]:altech.edu. Requests to be placed on the ASA mail
ing list should be addressed to [email protected]:altech.edu.
2. ASA algorithm
2.1. "Standard" simulated annealing (SA)
The Metropolis Monte Carlo integration algorithm (Metropolis et al., 1953)
was generali:.-;ed by the Kirkpatrick algorithm to include a temperature schedule
for efficient searching (Kirkpatrick et al., 1983). A sufficiency proof was then
showu to put an lower bound on that schedule as 1/ log( I;), where tis an artificial
time rneasme of the annealing schedule (Geman and Gernan, 1984). However,
independent credit usually goes to several other authors for indepell(lently de
veloping the algorithm that is now recogni:.-;ed as simulated annealing (Cerny,
HJ82; Pincus, 1970).
2.1.1. Boltzmann annealing (BA)
Credit for the first simulated <:umealiug is generally giveu to a Monte Carlo
irnportaw:e-sampling technique [or doing large-dirnensioual path iutegrals aris
iug in statistical physics problems (Nictropolis ct a!., 1953). This rnetlwd was
generalized to fitting uon-couvex cost-functions arising in a variety of problems,
e.g., finding the optimal wiring for a densely wired computer chip (Kirkpatrick
ct a!., Hl83). The choices of probability distributions described in this sectiou
arc generally specified as Bol tzmaun annealing (J3A) (Szu and Hartley, 1!)87).
The method of sirrmlat.ed aunealing consists of three functional relation~hips.
1
l. qr (:1:): Probability deusity of state-space of D parameters .7: = {.T ; i =
1, D}, where the ~ubscript T ~ignific~ a pararnctcrization popularly re
ferred to as the "temperature".
2. h.(6..E): Probability for acccptauc:e of new cost-functiou given the .iust
previous value.
3. T(k): schedule of "annealing" the "tcrnperatuw" T in aunealing-time
steps k, i.e., of changing the volatility or fluctuations of mw or both of the
35
Adaptive Hilllnlated annealing {ASA): lcs:;ons lcan1cd
The acceptance probability is based on the chances of obtaining a new state
with "energy" Ek+l relative to a previous state with "energy" Ek,
h(~E) +
exp( -Ek+dT) exp( - Ek/T)
1
1 + exp( ~E /T)
~ exp( -~E/T), (1)
where ~E represents the "energy" difference between the present and previous
values of the energies (considered here as cost functions) appropriate to the
physical problem, i.e., ~E = Ek+1 - Ek· This essentially is the Bolt~mann
distribution contributing to the statistical mechanical partition function of the
system (Binder and Stauffer, 1985).
This sampling algorithm also can be described by considering: a set of states
labeled by x, each with energy e(1:); a set of probability distributions p(x); and
the energy distribution per state d(e(x)), giving an aggregate energy E,
·~.::>(x)d(e(x)) =E. (2)
.''J;
The principle of maxirni~ing the entropy, S,
S =-LP(x) ln[p(x)/p(x)], (3)
X
where x represents a reference state, using Lagrange multipliers (Mathews and
Walker, 1!)70) to constrain the energy to average value T, leads to the most
likely Gibbs distribution G( :~; ),
1
G(:~:) = -z cxp( -H(:~:)/T), (4)
in terms of the normali~ing partition function Z, and the Hamiltonian H oper
ator as the "energy" function,
L
Z = exp( - H(:t:)/T). (5)
For such distribntions of states and acceptance probabilities defined by fum;
Lions such as h(~E), the equilibrium principle of detailed balance holds. I.e.,
the distributions of states before, G(xk), and after, G(:~;k+I), applying the ac
ceptance criteria, h(~E) = h(Ek+l - Ek) arc the same:
(G)
This is sufficicut to establish that all states of the r:;ystem cau be ::;arnpled,
36
L. INGBER
the temperature is changed, and so, at best, this must be done carefully awl
gradually.
An important aspect of the SA algorithm is to pick the ranges of the panun
ctcrs to be searched. In practice, computation of continuous systems requires
some discretization, so without loss of much generality for applications described
here, the space can be assumed to be cliscrctized. There arc additional con
straints that arc required when dealing with generating awl cost functions with
integral values. Many practitioners usc novel tedmiqucs to narrow the range
as the search progresses. For example, based on functional forms derived for
rnany physical systems belonging to the class of Gaussian-Markovian systems,
one could choose an algorithm for q,
(7)
where 6..r = .r - .T0 is the deviation of :~: from To (usually taken to be the jm:t
previously chosen point), proportional to a "moment urn" variable, and whc~rc
T is a measure of the fluctuations of the Boltzmann distribution q in the D
dimensional T-space. Given g(6..r-), it has been proven (Gernan and Geman,
1084) that it sufflccs to obtaiu a global minimum of E(:~:) if Tis selected to be
uot faster than
T(!-.:) = To (8)
Ink'
with To "large enough".
For the purposes of this paper, a heuristic demonstration follows, to show
that (8) will sulflce to give a global rninirnum of E(.1:) (Szu and Hartley, 1987).
In order Lo statistically a;;sme, i.e., requiring many trials, that auy point in
:~:-space c:u1 be sampled infinitely often in annealing-time (IOT), it suffices to
prove that the products of probabilities of not generating a state :~: IOT for all
amH:aliug-timcs succc:ssivc Lo k yield zero,
0
IT( )
1-.% = A: ..... oo 0. (0)
A:
This is equivalent to
L =".. ...
q" oo, (10)
00
k
as sec11 by taki11g the log of (0) and Taylor expa11ding i11 qk:· The problem tlwu
reduces to fi11ding T(k) to satisfy (10). Note that, giveu a very large space to
sample, often at best only a "weak" ergodicity cau be assmncd for this proof,
awl auy such crgodic:ity even for well-umlerstood physical systems is an opc11
area of research (Ma, 1085).
For 13A, if T(k:) is selected to be (8), then (7) gives
00 00 00
) nk > ) cxp(- lnk) = ) 1/k = oo. (11)
Adaptive :-=inndatcd annealing (ASJ\): lc:-=:-=ou:-= lc~\I'I!Cd 37
Al thougl1 t.lJCre arc tiound physical principles underlying the choices of (7)
and (1) (Metropolis ct al., 1!JG3), it was noted that this method of finding the
global mininuun in :~;-space was not limited to physics examples requiring bona
fide "tcmper:=ctnres" and "energies". Rather, this methodology can be re::cdily
extended to any problem for which a re::csonablc probability density h.(!:::..:~:) can
be fonr11d::tted (Kirkpatrick et al., 1983).
2.1.2. Simulated quenching (SQ)
l\llany resean:hers have found it very attractive to take advantage of the case
of coding and implementing SA, utili,.;ing its ability to handle quite complex
cost functions and constraints. However, the long time of execution of standard
Dolt,.;ntalllJ-type SA has many times driveu these projects to utili,.;e a tempera
tnrc sdwdulc too fast to satisfy the sufficiency conditions required to establish
a tmc (weak) ergodic search. A logarithmic ternperatnre schedule is cousistent
with the Dolt~mauu algorithm, e.g., the temperature schedule is taken to he
T, _ T lu k:o
(12)
k:- 0 lnk:'
where T is the "temperature", k: is the "time" index of annealing, and k: is
0
some startinp; iudcx. This can be written for large k: as
In k:of:::..h:
- To /;:(In k:)2, k » 1
, lnko
T" - 7o- .,.,--:-::: (13)
/;:(Ink: )2 ·
However, some researchers using the Dolt.,.;mmm algorithm usc cxponeutial sched
ules, <~.p;.,
f:::..T
- , = (r:- 1)/:::..k,/;: » ]
1
/,:
T,, = To exp((r: - 1)k), (14)
with expediency t.he only rcasmt given. vVhilc perhaps someday less stringent.
necessary conditions may he developed for the Dolt,.;rnauu algorithm, this is not
now th<~ state of affairs. The question arises, what is the value of this clear
misus<~ of the claim to use SA t.o help solve these problems/systems (Ingl><~r,
1!J!J3h)? Delow, a variant of SA, adaptive simulated annealing (ASA) (Inp;hcr,
l!JS!J; Ingber, l!J!J3a), in fact docs justify au exponential annealing sdtedulc, but
38
L. INGBER
2.1.3. Fast annealing (FA)
Although there are many variants and improvements made on the "standard"
Boltzmann algorithm described above, many textbooks finish just about at this
point without going into more detail about other algorithms that depart from
this explicit algorithm (van Laarhoven and Aarts, 1987). Specifically, it was
noted that the Cauchy distribution has some definite advantages over the Boltz
mann form (Szu and Hartley, 1987). The Cauchy distribution they define is
T
g(.6.x) = (l:,x2 + T2)(D+l)/2' (15)
which has a "fatter" tail than the Gaussian form of the Boltzmann distribution,
permitting easier access to test local minima in the search for the desired global
minimum.
It is instructive to note the similar corresponding heuristic demonstration,
that the Cauchy g(.6.x) statistically finds a global minimum. If (8) is replaced
by
T(k) = ~0, (16)
then here
Loo9 k To Loo 1
~ l:,xD+l k = 00· (17)
ko ko
Note that the "normalization" of g has introduced the annealing-time index k,
giving some insights into how to construct other annealing distributions. The
method of FA is thus seen to have an annealing schedule exponentially faster
than the method of BA. This m~?thod has been tested in a variety of problems
(Szu and Hartl<:!y, 1987).
2.2. Adaptive Simulated Annealing (ASA)
In a variety of physical problems we have a D-dimensional parameter-space.
Different parameters have different finite ranges, fixed by physical consider
ations, and different annealing-time-dependent sensitivities, measured by the
curvature of the cost-function at local minima. BA and FA have g distributions
which sample infinite ranges, and there is no provision for considering differences
in each parameter-dimension, e.g., different sensitivities might require different
annealing schedules. These are among several considerations that gave rise to
Adaptive Simulated Annealing (ASA). Full details are available by obtaining
the publicly available source code (Ingber, 1993a).
ASA considers a parameter o{ in dimension i generated at annealing-time
k with the range
(18)
39
Adaptive simulated annealing (ASA): lessons learned
calculated with the random variable yi,
a.k+1 a.~+ yi(Bi- Ai),
yi E [-1, 1]. (19)
Define the generating function
gD 1 =g D . .
(20)
.9r(Y) = 2(JyiJ + T;) ln(1 + 1/Ti) gT(y'),
where the subscript i on Ti specifies the parameter index, and the k-dependence
in Ti(k) for the annealing schedule has been dropped for brevity. Its cumulative
probability distribution is
g
·Yl ·YD D
Gr(Y) = ./_ ... ./_ dyt1 ... dytngr(yt) = C~(yi),
1 1
Gi ( i) = ~ sgn(yi) ln(1 + JyiJ/T;)
(21)
T y 2 + 2 ln(1 + 1/Ti) .
yi is generated from a u.i from the uniform distribution
u! E U[O, 1],
(v.; -
yi = sgn ~) T;[(1 + 1/T;)I2 11- 1]. (22)
1J.i-
It is straightforward to calculate that for an annealing schedule for T;
T;(k) = To; exp( - c;k1fn ), (23)
a global minima statistically can be obtained. I.e.,
1 ]1
00 oo [D
L.9k:::::; L IT 2JyiJci k = oo. (24)
ko ko t=1
It seems sensible to choose control over ci, such that
Toiexp(-mi) when kf = expn;,
c·,, mi exp( -ni/D ), (25)
where mi and ni can be considered "free" parameters to help tune ASA for
specific problems.
It has proven fruitful to use the same type of annealing schedule for the
acceptance function has used for the generating function g, i.e., (23) and (25),
but with the number of acceptance points, instead of the number of generated
points, used to determine the k for the acceptance temperature.
40
L. INGBER
In one implementation of this algorithm, new parameters ak+l arc generated
from old parameters o:i, by generating the 1/'s until a set of D an') obtained
satisfyiug the rauge constraints. In another alternative supported in ASA, useful
for some constraint problems, the 1/'s arc generated sequentially for each test
of the cost function.
2.2.1. Reannealing
\i\Theuever doing a multi-dimensional search in the course of a real-world non
linear physical problem, inevitably one must deal with different changing sensi
tivities of the ryi in the search. At any given annealing-time, it seems sensible to
attempt to "stretch out" the range over which the relatively insensitive parame
ters arc bc~ing seardwd, relative to the range::; of the more sensitive parameters.
It has proven fruitful to accomplish this by periodically rescaling the anneal
ing-time k, csscutially reaunealing, every hundred or so acceptance-events (or
at some user-defined modulus of the number of accepted or generated states),
iu terms of the sensitivities s; calculated at the most current minimum value of
the cost fuw:tion, 1_,
(26)
In terms of the largest s; = .Smax, a default rescaling is performed for each k; of
each parameter dimension, whereby a new index k'; is calculated from each k;,
k; --+ k' ;,
(27)
T;o is set to uuity to begiu the search, which is ample to span each parameter
dimension.
The acceptance tern perature is similarly rescaled. Since the initial acccp
taw:c temperature is set equal to an initial trial value of 1_, this is typically
very large relative to the current b<~st minirnmn, which may tend to distort the
scale of the region c:mrently being sampled. Therefore, when this rescaling is
performed, the initial acceptance temperature is rc:oct to the maximum of the
most current minimum and the best current miuimum of 1_, and the anncaliug
time index associated with this tcmpcratme is reset to give a uew temperature
equal to the minimum of tlw current cost-function and the absolnte values of
the current best awl last minima.
Also generated arc the "standard deviatious" of the theoretical forms, cal
culated as [D21_/(Doi)2]-112, for each parameter o;. This gives an estimate of
the "noise" that accompanies fits to stochastic data or functions. At the end of
the nm, the off-diagonal elcrneuts of the "variance matrix" arc calculated for all
parameters. This inverse curvature of the theoretical cost function can provide
tl]
2.2.2. Quenching
Auotlwr adaptive fcatme of ASA i~ its ability to perform qneuchiu12: in a HlC
thodical fasl1iou. Tl1i~ is applied by uoti1112: that the tcmperatme scbednle above
cau IH' reddiued as
1i (k;) 1'O i exp ( - r:; /,:;Q ;/D) '
r:;, 'f/1.; exp( - n;_Q;j D), (28)
in tenus of tlu~ "qnencbing factor" Q;. The above proof fails at (2tJ) if Q; > 1
as
LITD L
Jj/,:Q;/D = 1//;:Q; < 00. (20)
/.· k
Tl1is siJliple cnJcnlatiou shows how the "curse of dimensiouality" arises, awl
also gives a po~~ible way of living witl1 Ll1is di~case. Iu ASA, the infllleuce of
];u·l-!:e diuJeusious heconws clearly focnssed on the exponeutial of tl1<~ power of
/,: being 1/ D, as tl1e am1ealing req11ircd to properly smnpk the space becomes
prohibitively slow. So, if we cauuot conllnit reso11rces to properly Sct!llpk the
space er!-!:odically, then for son1e systems perhaps the ucxt best proced me wonld
be Lo tmn ou qnencbiug, w!Jerehy Q; <:<UJ hccoJnc 011 the order of the size of
nllJJJber of din1ensions.
1/
Tlw scale of the power of D ten1pcnttme sdwdnlc nscd for the acceptance
fnl!(:tiou cau he altered in a sii11ilm fn~l1ion. How<'V<~r, tl1i~ doc~ uot affect the
aune;tlinl-!: proof of ASA, <l!Jd so t!Jis Ittay nsed witl1011L dantal-!:illl-!: tlw (weak)
er12:odici ty property.
2.2.3. ASA applications
Tlw nhov<~ ddiiJC~ this 111d hod of ::tdapl".ive ~inndated a.nuealiug (ASA), previ
OIIsly called very fnst sii111Iiated rcauucaliu!!; (VFSll) ( 1I JI!;ber, 1080) ouly ImiJJCd
sl!di to <·m!l.rast it LlH' previons Iiidlwd of fast aunealiu12: (FA) (Sw and Hartley,
lD87). Tl1<~ aiJlte::tliul-!: sdwdnles for tlw tentperatmes T; decrease cxpoueutinlly
ill F\.11IJCaJii11!;-til!le /,:, i.e., 'fj = 'Jio cxp( -r:;/,:lfD). Qf <:Olll'~e, tl!C frtttcr the taiJ of
Ll1c /2:Cl1cratiu,e; fnnctiou, the SIIIFt.ilcr the ratio of acceptance to 1-!:Cilerated poiuts
in Llw Iii:. I-Towcvcr, in practice, it is fo111!d that for a 12:iveu gencrati11J.; fnw:tion,
t.l1is ratio is approxiiiJatel~r con~taut as tile lit fiud~ a glolml IlliuiiilllltJ. Tlwre
fore, for a la.rgc p;uaJJtcl.cr ~pace, tlw dficiell(:y of tl1e fit is dcten11iw:d by i.lw
aJl!I<:aiillJ.; ~c!H:dllk of tlw geueratiug f1I11ctiou.
A IIIFt.jor differeuc<: between ASA aud JJA algoritluiJ~ is that the ergodic
saiJipliug takes place iu au n+ 1 dimcn~ioual spa<:<:, i.e., in tcn11~ o[ n paramet.cr~
and tl1<: cost fJIII(:tion. Iu ASA the exponential aJHJealiug sclwdtdcs perwit
n:solln:es t.o he spent adapt.ivcly ou reaiJJH:aliug and ou paciug t.hc couvcrge11<:<'
iu all diJtt<~IJsioll~, e!lSIIJiii/2: ai!lplc global seardJiug iu Ll1e lir~t. phase~ of scarclJ
42
L. JNGBER
I have used ASA in several systems, ranging from combat analysis (Ingber,
Fujio, and Wehner, 1991; Ingber and Sworder, 1991), to finance (Ingber, 1990;
Ingber, Wehner et al., 1991), to neuroscience (Ingber, 1991), to a set of test
problems (Ingber and Rosen, 1992), to a new technique combining the power of
SA with the physics of large-scale systems (Ingber, 1992).
An optimization algorithm typically is just one tool used in a major project,
and many users of ASA use this code as just one such tool. Furthermore,
the author has made it a policy not to give out any information he receives
from users unless they specifically permit him to do so. Thus, many people
ask questions and give feedback on the code they might not give otherwise. In
order to get maximum feedback without unduly bothering researchers, the e
mail ASA_list is not open forum, but rather an efficient moderated medium to
gather information. Some published acknowledgments to use of the code are in
the asa_papers file of the ASA archive; the disciplines range from physics (Brown
et al., 1994; Tang et al., 1995) to neural networks (Cohen, 1994; Indiveri et al.,
1993) to difficult imaging problems (Wu and Levine, 1993) to finance (Wofsey,
1993), where most of the latter uses are kept proprietary.
3. ASA options
3.0.4. Options
ASA likely is the most powerful and flexible SA code presently available, be
cause the code has benefited from the feedback of the many users, and their
feedback has been used to add much to the code beyond the basic ASA algo
rithm described above.
The code has two basic modules in the ASA C-code, a user and an asa
module. All OPTIONS in the code have been tested to work with templates
provided in the user module. Feedback has developed a code which seems to run
well across many platforms, e.g., PC's, Macs, Crays, many UNIX workstations,
etc.
The emphasis in development of ASA has been to add power and flexibility
wherever possible. To make these extra features and code accessible to non
expert programmers, a "meta-language" of OPTIONS is used. Many of these
OPTIONS can be set in the provided Makefile, an asa_opt data file from which
to read in information, arguments passed to the compilation procedures, or in
the user module files (user.c, user.h , and asa_user.h, the latter being a bridge
between the asa and user modules).
A price is paid for this continual development of a more powerful and more
flexible code. The new user is presented With many OPTIONS, on the order of
a hundred. In many cases, when the ASA default OPTIONS work fine, only the
user's own call to his/her cost function is required. However, if these defaults
are not suitable for a particular system, then the user can become bewildered
by the many OPTIONS. If not much is known a priori about the system to be