Restarting Search Algorithms with Applications to Simulated Annealing

25 Pages·2000·0.43 MB·English

Applied Probability Trust (February 11, 2000)
RESTARTING SEARCH ALGORITHMS WITH APPLICATIONS TO
SIMULATED ANNEALING
F.MENDIVIL,R.SHONKWILER,ANDM.C.SPRUILL,(cid:3) Georgia Institute of Tech-
nology
Abstract
Someconsequencesofrestartingstochasticsearchalgorithmsarestudied. Itis
shown under reasonable conditions that restarting when certain patterns occur
yieldsprobabilitiesthatthegoalstatehasnotbeenfoundbythenthepochwhich
converge to zero at least geometrically fast in n. These conditions are shown
to hold for (RSA), restarted simulated annealing, employing a local generation
matrix, a cooling scheduleTn (cid:24)c=n, and restarting after a (cid:12)xed number r+1
of duplications of energy levels of states when r is su(cid:14)ciently large. For (SA),
simulatedannealing,withlogarithmiccoolingtheseprobabilitiescannotdecrease
to zero this fast. Numerical comparisons between (RSA) and several modern
variationson(SA)arealsopresentedandinallcases(RSA)performsbetter.
RANDOMSEARCH;STOCHASTICPROCESSRENEWAL;CODINGSCHEDULE;
MINIMIZATION
AMS1991SUBJECTCLASSIFICATION:PRIMARY60K20;60G35;60J20;65K10
SECONDARY
1. Introduction
Evidence that stochastic algorithms can spend excessive time in states other than
the goal comes most frequently and easily from simulations and more rarely from
theoreticalarguments. Forexample,an\optimal"coolingschedule(seeHajek(1988))
for simulated annealing (SA) guaranteesthat the probability the searchprocess is in
the goal state tends to 1 as the number of epochs n tends to in(cid:12)nity; the expected
(cid:3)Postaladdress: SchoolofMathematics, GeorgiaInstitute ofTechnology, Atlanta,GA30332.
1
2 F. MENDIVIL,R.SHONKWILER,M. C. SPRUILL
time taken by (SA) to hit the goal can however be in(cid:12)nite (see Shonkwiler and Van
Vleck (1994) or Fox (1995)). Our main results are Theorems 3.1 and 4.1 and they
show, among other things, that under certain conditions, restarted search processes,
particularlysimulatedannealingwiththepropercoolingschedule,has(cid:12)niteexpected
hitting time. Restarting (SA) based on di(cid:11)erent criteria from ours has been studied
numerically as in Atkinson (1992)and Nakakuki and Sadeh (1994). The undesirably
slow convergence of (SA) has also motivated some research like that in Kolonko
(1995)andB(cid:19)elisle(1992)onthe randomadjustment of the coolingschedule, onnon-
randomadjustments as reported,for example,in van LaarhovenandAarts (1987)or
Andresen and Nourani (1998), and the thorough theoretical treatment of simulating
direct self-loop sequences in Fox (1995) and the truncated version in Fox and Heine
(1995).
Most of the literature to date has been devoted to the investigationof the limiting
probability the process is found in the goal states and sometimes (see for example
ChiangandChow(1988))to therateofconvergenceoftheseprobabilities. Following
Shonkwiler and Van Vleck (1994), we take a somewhat di(cid:11)erent approach. The best
an algorithm has done up to epoch n is easily recorded and accordingly our results
are stated in terms of the probabilities of not yet seeing the goal by the nth epoch
rather than in terms of the probability of being in the goal state at the nth epoch.
We call the probabilitiesof not yetseeing the goalbythe nth epochtail probabilities
(not to be confused with measure-theoretic notions of the same name as found for
example in Niermiro (1995)). Theorem 3.1 provides su(cid:14)cient conditions for general
search processes to guarantee geometric convergence of the tail probabilities to zero
under restarting after user speci(cid:12)ed patterns and Theorem 4.1 specializes these to
the case of restarting simulated annealing. Of course, in both cases expected hitting
times must be (cid:12)nite.
Fix an arbitrary stochastic algorithm for minimizing a function over the (cid:12)nite set
C, a function f to be minimized, and let G be the goal subset determined by f.
Let X be the candidate point in C dictated by the algorithm at epoch n for what
n
van Laarhovenand Aarts (1987)would term this problem instance. The sequence of
randomvariablesX is aC-valuedstochasticprocessandits lawofmotion,whichof
n
Restarting Search Algorithms 3
course should depend upon f, can be quite general. The process X is not assumed
n
to be Markov,just adapted to a (cid:12)ltration fF :n(cid:21)1g. One could study alternative
n
restartingrulesbuttheoneanalyzedhereinvolvesrepeatedinstancesofenergylevels
f(X )ofstatesandhencehittingtimesofthevectorprocessofr+1stateswithvalues
n
inCr+1 =C(cid:2)C(cid:2)(cid:1)(cid:1)(cid:1)(cid:2)C calledther-process. Uponthe eventofther-processlying
in a certain (cid:12)xed subset of Cr+1 the algorithm is restarted at a randomly selected
point, always using the same distribution as the one used in the initial choice.
Theanalysisofthisschemeisshowntobeamenabletotechniquesfromdiscretetime
renewaltheory. Hu etal.(1997)hadpreviouslyutilized renewalequationsto analyze
restartingfromspeci(cid:12)csubsetsofthedomainC. There,theyshowedthatundersome
conditions, the tail probabilities decrease geometrically to 0. Their results extended
and related earlier results, based on the Perron-Frobenius theory, which Shonkwiler
andVanVleck(1994)obtainedforMarkovalgorithms. Therestartingprocessstudied
by Hu et al. is implementable in practiceonlywhen there is specialknowledge of the
function to be minimized and the algorithm. Restarting from an arbitrary state can
be detrimental to the progress of the algorithm towards the goal and one may not
know, a priori, which states are the correctones to initiate restarting. Restarting on
repeated energy values alleviates this problem.
Fox’s work on self looping in (SA) and especially the truncated versionin Fox and
Heine (1995) is closely related to the method of restarting (SA) on the diagonal,but
there are some key di(cid:11)erences. Foremost perhaps is that Fox’s calculation of the
objective function at each neighbor of the points visited is not required.
Andresen and Nourani (1998) propose a constant thermodynamic speed cooling
schedulefor(SA)andcompareitsperformanceforsolvingapermanentproblemwith
a variety of other cooling schedules. We present a comparison of our method with
these several methods on a permanent problem also.
The paper is organized by section as follows. In section 2 the restarted process
andther-processarede(cid:12)nedandthetailprobabilitiesareshowntosatisfyarenewal
equation. In section 3 rate of convergence to zero of the tail probabilities is studied.
Section 4 deals with simulated annealing where su(cid:14)cient conditions for geometric
convergence of the tail probabilities are presented. The outcomes of simulations
4 F. MENDIVIL,R.SHONKWILER,M. C. SPRUILL
comparing(SA)with(RSA)arepresentedinsection5alongwithsomecommentson
the selection of the parameter r.
2. The restart and r-processes
Restarting when a sequence of states lies in a subset D of Cr+1 de(cid:12)nes a new
process onthe originalsearchprocess. Details ofthis de(cid:12)nition andthe fact that the
tail probabilities for the r-process satisfy a renewal equation are presented here.
Weassumethatthe goalsetGisanon-emptysubsetofthe (cid:12)nite setC. Fix r (cid:21)1,
let D be a subset of Cr+1 and de(cid:12)ne for subsets A of C the sets
D =f(x ;x ;:::;x )2D:x 2Ag:
A 1 2 r+1 1
(cid:8) (cid:2)
Denote by E the set E(t)= x2CnG:P (X ;X ;:::;X )2D jX =x and
n n+1 n+r n
(cid:3) (cid:9)
allhistoriesupto n (cid:21)tforalln 6=;forsome(cid:12)xedt>0,andletU =G[E. Here
and below, by informal expressions like P[(X ;X ;:::;X ) 2 D j X = x and
n n+1 n+r n
all histories up to n](cid:21)t, we mean the more pedantic
P[(X ;X ;:::;X )2DjF ](cid:21)t
n n+1 n+r n
on the event X = x. Introduce the following two conditions, where T = minfn :
n E
X 2Eg.
n
(A1) 1>P[T <T ]>0.
E G
(A2) There is a (cid:12)nite K (cid:21) 1 and a number (cid:30) 2 (0;1) such that uniformly for
x2CnG, and all n,
P[T >m+njX =x and all histories up to x](cid:20)K(cid:30)m;
U n
where the probability on the left hand side is that the (cid:12)rst epoch after n at
which the X process lies in U is greater than m+n.
Let fX(j) :n(cid:21)1g be iid copies of the process X and de(cid:12)ne
n
n o
N(j)=min n:(X(j);X(j) ;:::;X(j) )2D or X(j) 2G
n n+1 n+r n
or N(j)=+1 if the set above is empty.
Restarting Search Algorithms 5
Lemma 2.1 Under (A2), P[N(j)<1]=1.
Proof. We have
P[N(1)>mjX(1) =x]
1
X
= P[X(1) =x(m+r);:::;X(1) =x(3); X(1) =x(2)jX(1) =x];
m+r 3 2 1
wherethe sum extends over x(2);:::;x(m+r) with no segmentof length r+1 in D.
With U as above de(cid:12)ne
K(x(1);x(2);:::;x(m))=minfj (cid:21)1:x(j)2Ug
or +1. Then
Xm
P[N(1)>mjX(1) =x]= p (x)P[N(1)>mjK =j]+q (x)
1 j m
j=1
Certainly, since there are more paths from x at 1 to x(m) at m which do not touch
U than there are which don’t touch U and don’t have segments in D
q (x)(cid:20)P[T >mjX(1) =x](cid:20)K(cid:30)m−1
m U 1
and since given U (cid:12)rst hit at j (cid:20)m, N(1)>m entails,
P[N(1)>mjK =j](cid:20)1−t:
Therefore,
P[N(1)>mjX(1) =x](cid:20)K(cid:30)m−1+1−t:
1
ChooseM solargethatM >r+1andK(cid:30)M−1+1−t=(cid:26)<1. IfA istheeventofa
j
runoflengthatleastr+1ofstatesin C forepochsn2[(M+r)j+1;(M+r)j+M]
then conditioning shows that
2 3
j\+k
P 4 Ac5(cid:20)(cid:26)k;
i
i=j+1
so that P[A i. o.]=1.
j
6 F. MENDIVIL,R.SHONKWILER,M. C. SPRUILL
TherestartedprocessY andther-process(cid:5)(n)cannowbewellde(cid:12)nedasfollows.
n
Let u(0)=0 and for j (cid:21)1
Xj
u(j)= N(i):
i=1
For m(cid:21)1 let
J(m)=minfj >0:u(j)(cid:21)mg:
The r-process is
(cid:16) (cid:17)
(cid:5)(m)= X(J(m)) ;X(J(m)) ;:::;X(J(m))
m−u(J(m)−1) m−u(J(m)−1)+1 m−u(J(m)−1)+r
and
Y =X(J(m)) :
m m−u(J(m)−1)
Ther-process(cid:5)(n)andrestartprocessY areperhapsbestillustratedbythediagram
n
below, where the rightmost entries in a row are in parentheses indicating that a
segment of r+1 states lying in D has occurredthere and that it is the (cid:12)rstsuch run
in that row. In the interests of readability of the tabular form, write x(i;j) = X(i).
j
Here we take r =2.
x(1,1) x(1,2) (x(1,3) x(1,4) x(1,5))
x(2,1) x(2,2) x(2,3) (x(2,4) x(2,5) x(2,6))
x(3,1) (x(3,2) x(3,3) x(3,4))
x(4,1) x(4,2) x(4,3) x(4.4)
y(1) y(2) y(3) y(4) y(5) y(6) y(7) y(8) y(9) y(10) y(11) y(12) y(13)
One has from the realization depicted above the values of the r-process (cid:5)(m) for
m=1;:::;4:
(cid:5)(1)=(x(1;1);x(1;2);x(1;3)); (cid:5)(2)=(x(1;2);x(1;3);x(1;4))
(cid:5)(3)=(x(1;3);x(1;4);x(1;5)); (cid:5)(4)=(x(2;1);x(2;2);x(2;3))
Our convention throughout is that if the goal is encountered then the next r are
taken as identical (and our interest in the process is terminated) and otherwise the
(cid:12)rst hitting time (cid:28) is de(cid:12)ned by
U
(cid:28) =minfn(cid:21)1:(cid:5)(n)2D g:
U U
Writing
u =P[(cid:28) >n]; (2.1)
n G
Restarting Search Algorithms 7
one has upon decomposition of the event f(cid:28) >ng as
G
f(cid:28) >ng=f(cid:28) >ng[(f(cid:28) >ng\f(cid:28) =1g)
G U G U
[(cid:1)(cid:1)(cid:1)[(f(cid:28) >ng\f(cid:28) =ng)
G U
that
Xn Xn
un =bn+ P[f(cid:28)G >ng\f(cid:28)U =jg]=bn+ un−jfj;
j=1 j=1
where f =P[(cid:28) >n, (cid:28) =n] and b =P[(cid:28) >n]. Therefore, the tail probabilities
n G U n U
u for the r-process hitting times satisfy a renewal equation.
n
3. Finiteexpectedhittingtimeandgeometricdecreaseoftailprobabilities
Rapid convergence of the tail probabilities u = P[(cid:28) > n] is established in this
n G
section. It follows from (A1) that
X
P[(cid:28) >j;(cid:28) =j]<1;
G U
j(cid:21)1
P
for P[(cid:28) > j;(cid:28) = j] = P[(cid:28) < (cid:28) ] and the latter probability can be seen to
j(cid:21)1 G U E G
be in (0;1) as follows. The di(cid:11)erences in the times at which (cid:5) hits D are clearly iid
and each such time di(cid:11)erence is the result of one of precisely three types of events;
either the X process has a succession of r+1 states for which the (cid:5) process is in
D , or a succession of r+1 states for which (cid:5) is in D , or one of the states in G
Uc E
is hit (once) at that epoch. Let q, p , and p denote the respective probabilities of
E G
these three types of hits, with q =1−p −p . These are numbers de(cid:12)ned bythe X
E G
process. The conditional probability the stop was of type E given that the stop was
of type U is
p
P[(cid:28) <(cid:28) ]= E
E G
p +p
E G
which under the condition (A1) is clearly not zero or one.
Theorem 3.1 Under the conditions (A1) and (A2), E[(cid:28) ] < 1, E[(cid:28) ] =
U G
E[(cid:28)U] <1,andthereisaγ 2(0;1)anda(cid:12)niteconstantcsuchthatγ−nu !c
1−P[(cid:28)E<(cid:28)G] n
as n!1, with u given by (2.1).
n
8 F. MENDIVIL,R.SHONKWILER,M. C. SPRUILL
Proof. To begin, we prove that 1 < liminfk!1fk−1=k. Since the event ((cid:28)G >
n) \((cid:28)U = n) implies the event ((cid:28)U > n− 1), one has fn (cid:20) bn−1 for all n and
liminfk!1fk−1=k (cid:21) liminfk!1bk−1=k. We will show below that liminfk!1b−k1=k >
1. Assuming this for the moment, it would follow that both of the power series
P P
Ψ (z) = b zn and Ψ (z) = f zn have radius of convergence greater
b n(cid:21)0 n f n(cid:21)1 n
than 1. On a su(cid:14)ciently small neighborhood of 0 one has
Ψ (z)=Ψ (z)+Ψ (z)Ψ (z) (3.1)
u b f u
P
whereΨ (z)= u zn andu isgivenby(2.1). SincethepowerseriesΨ (z)has
u n(cid:21)0 n n f
P
radiusofconvergencegreaterthan1and f <1,thefunctionΨ (z)=(1−Ψ (z))
n(cid:21)1 n b f
is holomorphic in an open set containing the unit disc. Therefore, Ψ (z) has radius
u
of convergence greater than 1 and it follows that there is a γ 2 (0;1) and a (cid:12)nite
constant c such that γ−nu !c as n!1. Furthermore, (see also Feller (1968))
n
X
b
E[(cid:28) ]= u = :
G n (1−f)
n(cid:21)0
To complete the proof decompose the event ((cid:28) > n), using H = minfj : (cid:5)(j) 2
U
DCnUg, as
[n
((cid:28) >n)= (((cid:28) >n)\(H =j))[(H >n)\((cid:28) >n):
U U U
j=1
Then
Xn
P[(cid:28) >n]= P[((cid:28) >n)\(H =j)]+P[(H >n)\((cid:28) >n)]
U U U
j=1
Xn
= P[(cid:28) >nj(H =j)\((cid:28) >j)]P[(H =j)\((cid:28) >j)]+c
U U U n
j=1
Xn
= P[(cid:28) >n−j]k +c :
U j n
j=1
P
Since j(cid:21)1kj =P[(cid:5) hits DCnU before DU]<1, if we can show liminfk!1c−k1=k >
1 then we’ll have by application to P[(cid:28) > n] of the arguments above applied to
U
un, that liminfk!1bk−1=k > 1. But cn is the probability of no restarts in the (cid:12)rst
n epochs so consulting Lemma 2.1, for n = k(M +r) + d, 0 (cid:20) d < M + r −1,
c (cid:20)(cid:26)k (cid:20)K(cid:17)n, where K =(cid:26)1=(M+r)−1 >1 and (cid:17) =(cid:26)1=(M+r) <1.
n
Restarting Search Algorithms 9
The constant c in Theorem 3.1 may be zero and the rate of convergence could be
fasterthangeometric. SometimestheseriesΨ (z)andΨ (z)satisfyconditionsunder
b f
which the rate can be established precisely as follows.
Corollary 3.2 If (A1), and (A2) hold, if f is not periodic and there is a real
n
solution (cid:18) > 1 to Ψ ((cid:18)) = 1 satisfying Ψ ((cid:18)) < 1 then there is a (cid:26) 2 (0;1) and a
f b
(cid:12)nite positive constant c such that (cid:26)−nu !c as n!1.
n
Proof. Taking (cid:26)=1=(cid:18), setting v =(cid:18)nu , the v satisfy a renewal equation
n n n
Xn
vn =Bn+ vn−jFj;
j=1
with B = (cid:18)nb and F = (cid:18)nf . Under aperiodicity of the f , Theorem XIII.10.1
n n n n n
of Feller (1968) yields the conclusion and the precise rate of convergence if Ψ ((cid:18)) =
b
P
B <1.
n n
By restarting, the expected time to goal of a search process can be transformed
from in(cid:12)nite to (cid:12)nite. Multistart, (see Schoen (1991)) where under no restarting
the hitting time is in(cid:12)nite with positive probability, is an obvious example. Under
simple conditions like restarting according to a distribution which places positive
mass on each state, multistart trivially satis(cid:12)es the conditions (A1) and (A2) with
t=1. FurthermoretheconditionsofCorollary3.2holdanditprovidesaninteresting
formula for the Perron-Frobenius eigenvalue as the reciprocal of the root of a low
degreepolynomial(seeHuetal.(1997)). Moresubtleisthe(SA)Example4.3below.
Of course it is easy to construct, as in the next example, a search process which
makesarestartingschemeperform,unlessitisdonecorrectly,worsethantheoriginal
algorithm.
Example 3.3 Let the state space be f0;1;2g with 0 as the goal. Suppose for all
n(cid:21)2
8
>>>><1 if i=0 and j =k =1
P[Xn+1 =ijXn =j;Xn−1 =k]=>>>>:1 if i=1 and j =1 and k 6=1
1 if i=1 and j 6=1
10 F. MENDIVIL,R.SHONKWILER,M. C. SPRUILL
and P[X = i j X = j] = 1 if i = 1 and 0 otherwise. Restarting whenever a
2 1
repeatedstateoccursisclearlydetrimentalinthiscasesincetheoriginalprocessgoes
deterministicallydirectlytothe goalasquicklyaspossible,butmustrepeatthestate
1 before doing so. Restarting after 3 the same however does not hinder the original.
4. Application to simulated annealing
ThegeneralityoftheconditionsinTheorem3.1suggeststhepossibilityofdiscover-
ing particular algorithms or classes of algorithmsexhibiting these generalunderlying
features and hence also ones for which restarting results in the rapid decrease of the
tail probabilities. One instance of this is simulated annealing. Let the generation
matrix Q be (cid:12)xed (see any standard reference on (SA) such as van Laarhoven and
Aarts for terminology) and introduce the following assumption (A3) which we will
invoke at the appropriate time. Let fN(x) : x 2 Cg be a collection of subsets of C,
one for each distinct x, which we call neighborhoods of x, having the property that
foranyx;y inC thereisa(cid:12)nitesequencez(1);z(2);:::;z(j)ofpointsinC suchthat,
settingx=z(0)andy =z(j+1),z(i+1)2N(z(i))foralli=0;1;:::;j. Weassume
also that y 2N(x) entails x2N(y).
(A3)-ThegenerationmatrixQallowsonlytransitionstotheimmediateneighborsso
that q = 0 unless y 2 N(x), assigns positive mass to each of them so that q 6=0
xy xy
for y 2N(x)nfxg and q =0.
xx
Consider the minimization of a function f de(cid:12)ned on the (cid:12)nite set C. Let
a=min min ff(y)−f(x):f(y)−f(x)>0g;
x2Cy2N(x)
where N(x) is the collection of immediate neighbors in C of x. Since C is (cid:12)nite,
(cid:11)>0. Let the acceptance rule be the usual with
(A4) theprobabilityofacceptingatransitionfromxto y 6=xatepochnbeing given
by minf1;e−(f(y)−f(x))=c(n)g, where c(n) is strictly decreasing to 0 with increasing n
andtheprobabilityofatransitionfromxtoxisoneminusthesumoftheprobabilities
of transitions from x to all other states.
Thebinaryrelationx(cid:1)yiff(x)=f(y)de(cid:12)nedonC(cid:2)Cisanequivalencerelationso