Table Of ContentRacingdownhill: optimizationand the random-field Isingmodel
D. Clay Hambrick
Department of Physics, Harvey Mudd College, Claremont, CA 91711 and
Department of Physics, Syracuse University, Syracuse, NY 13244
Jan H. Meinke and A. Alan Middleton
5
Department of Physics, Syracuse University, Syracuse, NY 13244
0
(Dated:January10,2005)
0
2 Thepush-relabelalgorithmcanbeusedtocalculaterapidlytheexactgroundstatesforagivensamplewitha
n random-fieldIsingmodel(RFIM)Hamiltonian. Althoughthealgorithmisguaranteedtoterminateafteratime
a polynomialinthenumberofspins,implementationdetailsareimportantforpracticalperformance. Empirical
J resultsforthetimingindimensionsd=1,2, and3areusedtodeterminethefastestamongseveralimplemen-
2 tations. Directvisualizationoftheauxiliaryfieldsusedbythealgorithmprovidesinsightintoitsoperationand
1 suggestshowtooptimizethealgorithm.RecommendationsaregivenforfurtherstudyoftheRFIM.
]
n I. INTRODUCTION in a capacitated graph [9, 10, 11]. Ogielski’s computations
n
of the RFIM ground state properties [11] utilized the push-
-
s relabel(PR) algorithmformax-flowintroducedbyGoldberg
i In systems with quenched disorder, changes in the ran- andTarjan[12]. ThePRalgorithmisinpracticethemostef-
d
dom background take place over time scales that are much
. ficientalgorithmavailableformanyclassesofproblems[13].
t longer than the time scale for evolution of the primary de-
a Although its generic implementation has a polynomial time
m grees of freedom. In magnetic systems, the spin degrees bound,itsactualperformancedependsontheorderinwhich
of freedom interact in an effectively frozen random en-
- operations are performed and which heuristics are used to
d vironment determined by substitutional disorder or vacan- maintain auxiliary fields for the algorithm. Even within this
n cies. The random-field Ising model (RFIM), defined as
polynomialtimebound,thereisapower-lawcriticalslowing
o ferromagnetically-coupled spins subject to a spatially vary-
down of the PR algorithm at the zero-temperature (T = 0)
c ing magnetic field, is a prototypicalmodelfor magnets with
[ transition[11,14,15].
quencheddisorderthathasbeenstudiedsince the1970s(for This paper presents results that are useful for minimizing
1 reviews, see, e.g., [1]). In dimensionsd > 2, the RFIM has
CPUtimeinRFIMground-statesimulations.Webeginwitha
v a transition between ferromagnetic (FM) and paramagnetic
briefreviewoftheRFIM,includingitsdefinitionandphases,
9 (PM) states; this transition can be found by varying either
in Sec. II. We then discuss the implementations for the PR
6
the temperature or the disorder, at sufficiently small values
2 algorithmsthatwestudy. ThePRalgorithmredistributesthe
of the remaining control parameter. Fishman and Aharony
1 randommagneticfieldamongspinsbypushingpositivefield
0 [2] mapped the RFIM with a field of random sign and fixed “downhill”withrespecttoapotentialfield(the“height”)de-
5 magnitude with disordered bonds to an experimentally real- fined for each spin. This redistribution and coalescence of
0 izable system: the dilutedantiferromagnetin a field (DAFF)
positive field and negativefield “sinks” allows for the deter-
/ [3]. At low temperatures, the glassy behavior of the DAFF
t mination of same-spin domains. The data structuresused to
a seen in bothexperimentandtheoryleadsto non-equilibrium
organize the pushes and the heuristics used for updating the
m
effectssuchashistorydependenceandabroadrangeofrelax-
potentialfieldaredescribedinSec.III.Theresultsforthetim-
- ation times. These observations are qualitatively consistent
d ing using various heuristics, summarized in Sec. IV, should
withpredictionsbybothFisherandVillainofanexponential
n beusefulfordesigningfurtherextensivestudiesoftheRFIM.
slowingdownnearthecriticalpointandofalowtemperature
o Visualizingthe auxiliaryfieldsleadsto a clearerexplanation
c phasedescribedbythezero-temperaturecriticalpoint[4,5]. ofthetimingresultsandwasimportantin guidingourwork.
:
v Exponentialslowingdownalsoaffectsoptimizationmeth- InSec.V,weusethesevisualizationstopresentaqualitative
i odssuchassimulatedannealing[6], thataremodeledonthe overviewoftheoperationsofthedifferentimplementationsin
X
dynamicsofthe physicalsystem. The largebarriersto equi- thedistinctphasesoftheRFIM.Theprimaryresultsofthepa-
r
a libration make it very time consuming to sample configura- per, namelythe recommendedchoicesfor the PR algorithm,
tion space accuratelyat finite temperaturesor to find the ex- aresummarizedinSec.VI.
actzero-temperaturegroundstatesusingsuchmethods.Find-
ingthepartitionfunctionfortheRFIMatfinitetemperatureis
NP-hard[7,8]. However,thezero-temperatureFM-PMtran- II. RANDOM-FIELDISINGMODEL
sition is expected to be in the same universality class as the
finite temperaturetransition. So it is fortunate that there are Therandom-fieldIsingmodelhasanon-trivialgroundstate
alternatemethodsforquicklyfindingtheexactRFIMground duetothecompetitionbetweentheferromagneticinteraction
state. Thesemethodsarebasedonamappingfromtheprob- thattendstoalignneighboringspinsandtheinfluenceofran-
lem of finding the ground state of the RFIM to that of find- domfields, whichtendstoforcespinstopointinrandomdi-
ingthemaximumflow—or,equivalently,theminimumcut— rections.Takingthestrengthoftheferromagneticinteractions
2
betweenneighboringspins to be J anddesignatingthe local andcomputerscience texts[21]. Inthis paper,we limitour-
randomfieldsbyh ,theenergyofaspinconfigurationis[1] selves to a description of the algorithm, neglecting proofs
i
of correctness, butincludinga “physical”descriptionfor the
H=−JXsisj −Xhisi, (1) variantofthePRalgorithmwehaveused.
Intuitively,thePRalgorithmfindsthedomainsofuniform
hiji i
spininthegroundstatebyrearranging(“pushing”)themag-
where the spin s at a givensite i on a d-dimensionallattice netic field. If the bond between two spins is strong enough,
i
with n = Ld sites can take on valuess = ±1 and the sites thefieldononespincanberemovedfromonespinandadded
i
iandj inthesuminthefirsttermarenearestneighbors. We toitsneighbor,possiblyinfluencingthedirectionoftheneigh-
studytheGaussianRFIM,wheretheh areindependentvari- boringspin. Thisrearrangement(andchangeinthebonds,as
i
ables chosen from a Gaussian distribution with mean 0 and described below) is the push. Subsequent pushes can then
variance ∆2J2, with periodic boundaryconditions. The pa- affect distant spins. As a consequence of pushes, positive
rameter∆characterizesthestrengthofthedisorderrelativeto and negative fields originally located on separate spins can-
theferromagneticinteraction. cel. The cancellation leads to domains of uniform sign for
Asthedisorderdominatesoverthermalfluctuationsatlarge the excess fields. This domain growth by rearrangement of
length scales and low temperature in more than two dimen- fieldislimitedbythestrengthofthenearest-neighborbonds,
sions,thegroundstatesofthismodelareofinterest.Indimen- which“carry”therearrangement.Thepushoperationreduces
sions d > 2, there is a zero-temperature transition between orremovestheinteractionsbetweenneighboringspins. When
twophasesatthecriticaldisorder∆ = ∆ . When∆ < ∆ , thefieldhaslargevariationscomparedtothebondstrengths,
c c
theferromagneticinteractionbetweennearestneighborsdom- themagneticfieldcannotbepushedveryfarasbondsbecome
inatesandthespinstakeonameanvaluem=n−1 s with saturatedandblockfurtherrearrangement.Conversely,inthe
Pi i
|m|=6 0inthelimitn→∞.Inthecase∆>∆ ,randomness limit of weak fields, large domains form, as the bonds fa-
c
dominatesandthegroundstateis“paramagnetic”|m|=0,as vor alignment of nearest neighbor spins: in the language of
n→∞. PR, the large capacity of the bonds relative to the strength
Indimensionsd = 1and2,thegroundstateisinthepara- ofthefieldsallowsforlong-rangerearrangementsoftheran-
magnetic phase at any ∆ > 0 in the thermodynamic limit dom field. In the limit of very weak fields, the field can be
(i.e.,∆ =0).Butthecorrelationlengthξ,characterizingthe pushed anywhere (no bonds are saturated), so there is only
c
rangeofspin-spincorrelationsorthesizeofuniformspindo- one domain, whose orientation is simply a result of the sum
mains,divergesas∆−2whend=1(see,e.g.,[16]).Forsam- oftherandomfieldsonallthespins.Butwhenthefieldisvery
plesofsizeL≪ξ,i.e.,for∆<∆ ∼L−1/2,where∆ (L) strong(∆ ≫1),therearrangementoffieldislimited,and,in
x x
is the sample-size-dependent crossover disorder, the ground mostcases,thedomainshavesinglespinsandtheorientation
stateisessentiallyferromagnetic.Ithasbeenarguedthatln(ξ) ofaspinsisinthesamedirectionasitsmagneticfield.
scalesasaninversepowerof∆whend=2[17,18,19]. This Theotherbasicoperationistherelabeloperation. Thisop-
very rapidly growing correlation length gives rise to an ap- erationupdatesanauxiliaryfield,theheight,definedforeach
parent ferromagnetic phase in d = 2 at even moderate ∆ spin. Following the more detailed constraints described in
x
in finite samples [20]. When we take the limit ∆ → 0 in Sec. IIIA, this height guides the pushes. Most simply put,
d = 1andd = 2, we willbetaking∆ ≪ ∆ . Manyofthe pushesarealwaysin thedownhilldirection. Whena pushis
x
results for the algorithm, such as which algorithm is fastest, notpossiblefromagivensite,therelabeloperationincreases
hold independentlyof the existence of a true physicalphase themagnitudeoftheheightofthatsite. Thisrelabellingwill
transition. eitherallowapushtobeexecutedorwillidentifythespinas
havingaparticularsigninthegroundstate.
The basic operationscan be carried out using a variety of
orderingmethods. Thechoiceofmethodaffectsthe running
III. PUSH-RELABELALGORITHM:DATASTRUCTURES
timeofthealgorithm. Aspecificalgorithmisdefinedbytwo
ANDHEURISTICS
sets of choices, which are defined and described in detail in
thefollowingsubsections:
WenowdiscussthemotivationforusingthePRalgorithm
and outline its structure. In particular, we define the auxil-
1. the dynamically-determined order local operations
iaryfieldsandbasicoperationsthatoperateonthosefieldsto
(pushes and relabels), which is organized by a choice
determine the ground state. Picard and Ratliff showed that
ofdatastructure,and
any quadratic optimization problem, such as the RFIM, can
be mapped onto a min-cut/max-flowproblem [9]. There are 2. heuristicmanipulationsofthe auxiliaryfields, namely,
a number of algorithms for solving max-flow [21], though globalupdatesandgaprelabeling.
CherkasskyandGoldberg’sresultsshowthatthePRalgorithm
[12]isoftenthebestalgorithmforsolvinglargeproblemson WeimplementedthePRalgorithminJavaandC++. TheJava
avarietyofgraphs[13]. Fordetailedexplanationsofthemap- implementation can visualize the evolution of the auxiliary
pingoftheRFIMtoamax-flowproblemandtheproofsofthe fieldsusedbyPR.TheC++codereliedontheoriginalCcode
correctnessofthePRalgorithm,seereviewsofapplicationsof [25]developedbyCherkasskyandGoldberg[13]. Thecodes
combinatorialoptimizationto statistical physics [22, 23, 24] canbedownloadedfromourwebsite[26].
3
A. Auxiliaryfieldsandpush-relabeloperations priorityqueueimplementedasaheap(HPQ),alowestheight
priorityqueue(LPQ)alsoimplementedasaheap,andadou-
The PR algorithm uses three auxiliary fields to guide the bleFIFOqueue(DFIFO)thattreatspositiveandnegativeex-
rearrangementof the magnetic field and to enforce the con- cess symmetrically. (We also tried a stack or last-in first-out
straints given by the bond strengths. One field is the resid- (LIFO)structure,butrejectedit duetovastly longerrunning
ualbondstrength(morecommonlyreferredtoastheresidual timesatlarge∆.)
capacity in the literature [21]). This residual interaction be- TheFIFOstructureisalistofactivesitesinthelattice. The
tween sites is denoted r . Initially, r = r = J for all site at the front of the list, i, is removed. Excess is pushed
ij ij ji
nearestneighborpairs(i,j),butingeneralr neednotequal away from i if possible. If any inactive neighborsj of i are
ij
r duringtheexecutionofthealgorithm. Thisresidualbond made active through a push operation, they are added to the
ji
strength definesthe pathsalong which excessmagnetic field endofthelist. Iftherearenoaccessibleneighbors,iisrela-
canbepushed. Asitej issaidtobereachablefromasiteiif beled. IfiisstillactiveafterthisPRstep,itisalsoappended
thereisadirectedpathfromitoj withr > 0forallbonds totheendofthelist.
kl
(k,l)alongthatpath. TheHPQ structure[27] is morecomplexthan FIFO. Like
For each site i, an excess field e and a height field u the FIFO structure it contains a list of all active sites. The
i i
(whichis oftencalled“distance”or“rank”)are alsodefined. listisnotorganizedbythetemporalorderinwhichsiteshave
At the beginningof the algorithm, the excess, which can be beentreated,aswouldbethecaseinaFIFOqueue,butbythe
positiveornegative,issetequaltotherandomfieldstrength, heightof the sites. The first site in the HPQ is always a site
e = h , and the height field u is set equal to the dis- withmaximalheight. Whenwe removeasite fromthefront
i i i
tancetothenearestreachablesitewithnegativeexcess. Sites ofthequeue,thesitewiththehighestheightthatisstillinthe
fromwhichnositewithnegativeexcesscanbereachedhave queuemovestothefront.Ifweaddasitewithalargerheight
u =∞. Siteswithnegativeexcesshaveu =0. thananyofthesitesalreadyinthequeue,thenewsitemoves
i i
ThePRalgorithmmaintainstheheightfieldu inafashion tothefrontofthequeue.IfafteraPRstepasiteisstillactive,
i
designed to move the positive excess “downhill” (to smaller itisre-addedtothequeue. SinceaPRsteprelabelsanactive
heights), i.e., towards sites with negativeexcess, where pos- site by increasing its height by at least one before adding it
sible. If it becomes impossible to rearrange excess towards back to the queue, such a site is still of maximalheight and
a site with negative excess, the site is given a height label isaddedatthefront. Thusthesamesitemightbeactedupon
u = ∞ (in practice, ∞ is represented by n, the number of by the algorithmmany times in a row. In practice, the HPQ
i
spins). A site j issaidto beaccessiblefromi iftheresidual issimplyimplementedbysortingthesitesintobinsgivenby
bondstrengthr > 0andu = u +1. Anysitewithposi- theirheightvalues.
ij i j
tiveexcessandheightu <∞isactive. (Inthedoublequeue The LPQ structure [28] is exactly the same as the HPQ
i
methodmentionedinSec.IIIB,negativeexcesssitescanalso structure, except that the list is reversed, so that sites i with
beactive.) lowestheightuiaresubjecttoPRoperations.
Thepushoperationmovesexcessfromanactivesiteitoan FIFO,HPQ,andLPQenforceanasymmetrybetweenpos-
accessibleneighborj. Lettingδ = min(r ,e ),theexcesses itiveandnegativeexcess. Whyshouldonesign(positiveex-
ij i
areupdatedbye → e −δ,e → e +δ,whiletheresidual cess) be pushedaroundwhile the other (negativeexcess) re-
i i j j
bond strengths are updated according to r → r −δ and mainsstatic(exceptbyrearrangementofpositiveexcessonto
ij ij
r →r +δ. Ifanactivesitehasnoaccessibleneighbor,so anegativeexcessoflessermagnitude)? Certainly,thephysi-
ji ji
thatnopushesarepossible,itisrelabeled:theheightuiisset calproblemisunalteredbythereplacementhi → −hi. Our
to one greater than the height of the lowest neighbor j (i.e., fourth implementation treats positive and negative excesses
minimalu overneighborsj) forwhich r > 0. If nosuch onanequalfooting. Negativeexcessthenmovesfromlower
j ij
neighbor exists, the height of i may immediately be raised heights (here, ui < 0 is a possible height label) to higher
to u = ∞. We call the combination of all possible push heights. WeimplementedthisastwoFIFOqueues(DFIFO),
i
operationsfromasinglespinpossiblyfollowedbyarelabela oneforpositiveandonefornegativeexcesssites. Bothqueues
PRstep. areupdatedsimultaneously. Asnotedinthenextsection,we
The PRalgorithmterminateswhenno activesites remain. alsousedtwoheightfields,oneforeachsignofexcess,when
The total number of PR steps needed to complete the algo- implementingDFIFO.
rithm is denoted by N . The assignment of spin orienta-
PR
tionsinthegroundstateisfoundbyexecutingaglobalupdate
(Sec. IIIC) to finalize which sites have u = ∞. Sites with C. Heuristics
i
u = ∞areassigneds = 1,whiletheremainingspinshave
i i
si =−1. If these queues are adapted as described so far, the algo-
rithm,thoughpolynomialinn,istooslowtobepracticalfor
studyinglargersystems. Heuristicscanbeusedtomanipulate
B. Datastructures the height field to guide the push-relabel operations. Good
heuristicsarecrucialtothepracticalityofthealgorithm.
WeimplementedthePRalgorithmusingfourdifferentdata All our implementationsuse the “globalupdate” heuristic
structures: a first-in-first-out queue (FIFO), a highest height for initialization of the u . A global update [13] is used for
i
4
twopurposes. Itcreatesgradientsintheheightfieldfromthe ing of the data structure, nor do they reflect the time needed
siteswithpositiveexcesstothenearestsitewithnegativeex- for the heuristics, such as the global update. The CPU time
cess, allowing the excess fields to move efficiently towards tisthereforeanotherusefulmeasureoftheperformanceand
annihilation. Itaddition,itidentifiesregionsthatare discon- ultimately what we want to minimize. For the timing mea-
nectedfromtherestofthespinsanddonotcontainanyneg- surementweusedSUN’sJavaVirtualMachine1.4.2onDual
ativeexcess; suchregionsare labeledasspin-upregionsand 1GHzPIIImachineswith512MBofRAMrunningLinux.
are removedfromfurtherconsideration. Withoutthis identi- Herewesummarizetheresultsofourtimingruns. Wefirst
ficationitwouldtakeoftheorderofnr operationstomarka describeourresultsfortwo-dimensionallattices. Inthiscase
regionofr spinsasup. Theglobalupdateisimplementedas thereisnophasetransitionfortheRFIMgroundstatesatfinite
a breadthfirstsearchstartingfromthesetof negativeexcess ∆. However,asthecorrelationlengthξdivergesveryrapidly
sites and takes of the order of n operationson a hypercubic as∆decreases, thereis acrossoverdisordervalue∆ (L)at
x
lattice. Duringthissearch,theheightofeachvertexissetto whichξ exceedsthe linear system size L. For∆ < ∆ (L),
x
thedistancetothenearestreachablesitewithnegativeexcess the 2D system is effectivelyferromagnetic. The dependence
(sink). One important result of global updates is that local of∆ onLisveryslow[17,18,19]. Inoursimulations,we
x
minimaoftheheightfield withpositiveheight,oftenresidu- find∆ (8)≈1.7and∆ (4096)≈0.55.Inthreedimensions,
x x
alsofformernegativeexcesssites,areeliminatedatsiteswith thereisatruetransitionat∆=∆ ≈2.27[14].
c
nonnegativeexcess. Followingcommonpractice, we choose Inthissection,wecomparethetimingsforFIFOandHPQ
theintervalbetweenglobalupdatestobefixed,withaglobal structures and then present results for DFIFO. (Results for
updateexecutedaftereveryΓpush-relabel(PR)steps. LPQ are includedin the discussion of the results for d = 3,
The global update needs to be modified for the double near the end of thissection.) When seeking to minimize the
queue(DFIFO)approach. Inparallelwiththeheightfieldfor runningtime, thereisa competitionbetweenfrequentglobal
positive excess, based on a search from nodes with negative updates (small Γ), which improve the efficiency of the PR
excess, it is naturalto construct the height field for negative steps and infrequent updates (large Γ), which save the cost
excesssitesusingpositiveexcesssitesas“sinks”.Thiscreates ofperformingaglobalupdate. We firstdeterminethe global
anambiguityforsiteswithzeroexcess.Shouldtheirheightbe updateintervalthatminimizestheCPUtimeoverallchoices
determinedinrelationtothepositiveornegativeexcesssites? of∆.Generally,foralldatastructures,wefindaminimumin
Wedidnotdirectlyresolvethisambiguity,butinsteaduseda tforΓ ≈ n(oratleastthatperformanceisnotimproved
min
schemewith separateheightfieldsforthe positiveandnega- forotherΓ). TheFIFOandHPQstructuresthengivecompa-
tiveexcesses,withseparaterelabelingandglobalupdates. rableresultsawayfromthecrossover∆ ≈ ∆ ,thoughFIFO
x
Inaddition,whenusingtheHPQdatastructure,weusegap ispreferableatsmall∆andHPQisfasteratlarge∆. DFIFO
relabeling [13]. If there is a heightu , such that no site has takessignificantlymoretimethaneitherHPQorFIFOtofind
g
height u = ug, but there are active sites with u > ug, i.e., thegroundstate for∆ontheorderof∆x. We alsofindthat
thereisagapinthesetofheights,allsiteswithu>u areas- the discretenessofglobalupdatescan leadto plateausin the
g
signedthemaximumheight. Thisreducestheneedforglobal numberofglobalupdatesN asafunctionofΓ. The3Dre-
G
updatesasactivesiteswithu>u aredisconnectedfromthe sultsaregenerallyconsistentwithourresultsfor2Dlattices.
g
sinks. When using HPQ, regions with larger height tend to WefindthatLPQalwaysperformsatleastaswellasHPQand
have their height raised uniformly, leading to the possibility isfasterfor∆>∆ .
c
ofcreatingagapquickly. Asthegapisaglobaltally,itisnot
efficient at detecting regions that have separated themselves
fromtheirsurroundingslocally,sogaprelabelingislessuseful A. Choosingtheglobalupdateintervalintwodimensions
in, e.g., FIFO, whereupdatesare carriedoutwithoutrespect
toheight.Tofacilitategaprelabeling,theHPQalgorithmuses We first examine the effect of varying the update interval
an array that containsthe numberof sites (whetheractive or Γfortwodimensionalsamples. Globalupdatesareonlyuse-
not)ateachpossibleheight. Agapiscreatedwhentheoccu- fulif they reducethe numberof PR steps needed to find the
pationnumberatanon-maximalheightisreducedtozero. solution. As we mentioned in Sec. IIIC, global updates are
expensive, taking of the order of n operations, and should
thereforenotbeperformedtoofrequently.Tofindtheoptimal
IV. COMPARISONOFIMPLEMENTATIONS updateintervalΓ ,wevariedΓforeachdatastructureover
min
severalordersofmagnitudeforvariousvaluesoftherandom
Tofindthebestavailablecombinationofdatastructureand field strength ∆ and examined how both NPR, the number
heuristics for solving the Gaussian RFIM, we measured the of push-relabelcycles used to find the ground state, and the
runningtime of the algorithm for a variety of combinations. runningtimetvaried.
We measure the running time for the algorithm to find the
ground state using both CPU time t in seconds and the to-
talnumberofPRstepsN . PRstepsarethecoreoperations 1. GlobalupdateintervalforFIFO
PR
ofthealgorithmandgiveusamachineindependentmeasure
oftheimprovementduetotheheuristics. However,PRsteps Immediately after a global update, push operations are
don’taccountforthe timeneededforthe internalbookkeep- guaranteed to move excess towards the nearest sinks. After
5
aglobalupdate,however,somesinksaresoonannihilatedby 1e+07
excesspushedintothem. Ifasinkatsitesisannihilated,the L=25
heightfieldaroundsnolongerindicatestheshortestdistance L=50
L=100
toasink;theheightfieldstillslopestowardssalthoughthere 1e+06 L=200
is no longer a sink. The extent of this “misleading” height L=400
field depends on the density of sinks and is limited by the
distancebetweensinks. WhenΓissmall, globalupdatesare PR1e+05
N
frequentenoughthatthe heightfield generallyleadspositive
excess towards a sink. As long as this is the case, increas-
ingthefrequencyofglobalupdatesdoeslittletoreduceN .
PR 10000
Therefore,forsmallΓ,weexpectN tobealmostindepen-
PR
dentofΓforFIFO.
Fig. 1 displays algorithm costs N and t for the case
PR 1000
∆ = 2.2, which exceeds ∆ . The expectation that N is 1 10 100 1000 10000 1e+05 1e+06 1e+07
x PR Γ
independent of Γ at small Γ is confirmed by the plot of the
meanofN asafunctionofΓdisplayedinFig.1(a).Aswe
PR
increaseΓ, alargerandlargerfractionoftheheightfieldbe- 1000
comesincorrectbetweenglobalupdates: moreminimaofthe L=25
L=50
heightfieldnolongercontainsinks. TheplotshowsthatNPR 100 L=100
startstoincreasesignificantlyaboveΓ ≈ 0.1n. Tominimize L=200
L=400
therunningtimet, we needtobalancethecostofadditional
10
global updates with the reduction in N . A global update
PR
takesoftheorderofnoperationsand,infact,Fig.1(b)shows
s] 1
aminimumintvs.ΓatΓmin ≈ n. Toverifyourassumption t[
that Γ ≈ n, we rescaled the curvesin both Fig. 1(a) and
min
0.1
(b)bydividingΓ,N ,andtbyn. Thecollapseofthedata
PR
verifiesthe minimumin t vs. Γ atΓ ≈ n andshowsthat
min
0.01
the runningtime of the algorithmscales nearlylinearly with
nin2D(Fig.2),atfixed∆.
0.001
10 100 1000 10000 1e+05 1e+06 1e+07
Γ
2. DetailedlookattheFIFOdata
Figure 1: Mean running time for finding the ground state of a 2D
While determining the optimal global update interval for RFIMsystemvs.ΓusingFIFOand∆=2.2,wherethecorrelation
theFIFOstructure,wenoticedtwofeaturesofinterestinthe length is much smaller than the sample sizes used, ξ ≪ L. For
data.ThesefeaturesarethetendencyforN tobeaninteger clarityinthisfigureandotherfiguresinthissection,statisticalerror
PR
multiple of Γ and, at small ∆, a separation of the mean run barsarenot shown, but areconsistent withtheapparent deviations
fromsmooth curves. Thefigureshows(a)thenumber of PRsteps
timesbetweenpositivelyandnegativelymagnetizedsamples.
A close look at the data reveals piecewise linear behavior NPRusedtofindthegroundstateand(b)CPUtimet[s]asafunction
oftheglobalupdateintervalΓforL = 25, 50, 100, 200,and400.
intheaveragerunningtimeofthealgorithm,whenmeasured
ForsmallΓ,NPR isnearlyindependentofΓshowingthatfrequent
byN ,asdisplayedinFig.3. Theselinearregionsarecon-
PR global updates are unnecessary. At Γ somewhat less than 0.1n =
sistentwithNPR =kΓ,forintegerk,asshownbythedashed 0.1L2,NPR startstoincrease. TheminimumintheCPUtimetis
linesinthefigure.Theselinearregionscoincidewithplateaus atΓmin ≈n,wherethechangeintimeneededfortheadditionalPR
inthenumberofglobalupdatesexecutedduringthesolution, stepsbalancesthechangeintimeneededfortheglobalupdate.
N ,whenplottedvs.Γ(seetheinsetinFig.3). Theplateaus
G
in N vs. Γ reflect the effect of the global update, which
G
causeslargechangesintheheightfieldandcanbringtheaux- operationsidentifytheup-spinregions. WhenΓisfinite,but
iliary fields close to a solution. Note that a ground state is largerthanthenumberofPRstepsneededtofindallofthedo-
foundbythealgorithmwhenallofthepositiveexcessiscon- mainboundariesandsmallerthanthetotalnumberofPRsteps
finedtoregionsthatareisolatedfromsinks. Thisisolationis neededto terminatethe algorithm, the first globalupdateef-
due to saturated bonds that block the rearrangementof flow fectivelyterminatestherunningofthealgorithmandN =1.
G
(and to the cancellation of positive excess with negative ex- There is another interval at smaller values of Γ where one
cess in regions accessible to the sinks). The algorithm will global update is executed before the domain boundaries are
notterminateuntiltheblocked-offregionshavetheirheights determinedandthesecondglobalupdateterminatesthealgo-
raised to u = n. Without global updates (Γ = ∞), there rithm,givingN =2.Thispatterncontinuestohighervalues
i G
is a separation between the times when the push-relabelop- ofm. Weemphasizethatthedatapresentedisaveragedover
erationseffectivelydeterminethe domainboundariesbysat- 100 samples at each value of Γ; the data therefore indicates
urating the appropriate bonds and the time when the relabel thatthe fluctuationsin the locationsof the linear regionsare
6
1e+07 1e+07
L=50 10
L=100
L=200
L=400 NG
1e+06
n
/R PR1e+06 1
P N 10000 1e+05 1e+06 1e+07
N Γ
1e+05
10000 1e+05
0.01 0.1 1 10 100 10000 1e+05 1e+06 1e+07
Γ/n Γ
10 Figure 3: Plot showing the piecewise linear behavior of NPR, the
meannumberofPRstepsneededtofindthegroundstate,whenplot-
L=25
L=50 tedasafunctionofΓ,theintervalbetweenglobalupdates,forFIFO
L=100 in d = 2. The parameter values are L = 200 and ∆ = 2.2 and
L=200
L=400 thedataisaveragedover100samples. NPR increaseslinearlyover
intervalsinΓ, NPR = kΓ, withthecurvesfork = 1,2,3,...in-
1 dicatedbythedashedlines. Thelargechangescausedbytheglobal
n
]/ updateoperationcanleadtoplateausinthenumberofglobalupdates
s
t[ NGvs.Γ,asshownintheinset.
globalupdateisincluded,theheightlabelofallsitesissetto
0.1 its maximum value once no more sites with negative excess
exist. As the global relabeling takes only of the order of n
steps,thismakestherunningtimesforup-anddown-magne-
0.01 0.1 1 10 100
Γ/n tizedsamplesmoresimilar,towithinO(n)intotalmagnitude,
forΓ=O(n).
Figure 2: Scaling of PR steps NPR and CPU time t with system
size n for FIFO (∆ = 2.2). (a) The average number of PR steps
3. GlobalupdateintervalforHPQ
persite,NPR/n, asafunctionofnumber ofpush-relabel stepsper
siteperglobalupdate,Γ/n. Thedatacollapseisconsistentwiththe
number ofoperationsscalingwithn(thebestthatcanbeexpected For ∆ < ∆ , the behavior of the timing for HPQ vs. Γ
x
forξ≪L).(b)CPUtimetpersite(t/n)vs.Γ/n.Thevaluesalong is very similar to the behavior of the timing for FIFO, with
bothaxeshavebeendividedbynthenumberofsites.Afaircollapse Γ ≈ n. In contrast, in the paramagnetic regime, gap re-
opt
developsforL≥100.TheminimumisatΓmin/n≈1,independent labelingdetectsencloseddomainsefficientlyandquicklyre-
ofsystemsize.Theminimumisshallow,however,andevenmissing
ducesthenumberofactivesites. Forlarge∆, then,frequent
Γminbyafactoroftenincreasestherunningtimeonlybyafactorof
globalupdates(Γ < n)quicklydominatetherunningtimet.
abouttwo.
ForΓ>n,tbecomesindependentofΓ(Fig.4),asglobalup-
datesarenotexecutedbeforethesolutionisfound.Consistent
with theoptimalbehaviorforsmall∆ andthe independence
quitesmall.
of t from Γ at large ∆, we will use Γ ≈ n for the HPQ
Wealsonotedastrongup-downasymmetryintherunning opt
structure. AsΓ = ∞isareasonablechoiceforlarger∆,we
time at small ∆. For ∆ ≪ ∆ and Γ = n, it takes about
x willsometimesincludethischoiceforcomparison.
50%longertofindthesolutionifthe groundstate isspin-up
than if it is spin-down. For Γ → ∞ the ratio of mean run-
ningtimesbecomesverylarge. Atsmall∆, thegroundstate
is determined by the sum over all random fields. If the sum B. FIFOvs.HPQintwodimensions
is negative, the algorithm is done as soon as all the positive
excessfieldshavebeenannihilated,butifthesumispositive Now that we have found an optimal value for Γ for each
allthesiteswithremainingexcessfieldsmustbemovedupto data structure, we can directly compare the performance of
height n before they are labeled inaccessible. This local re- FIFO and HPQ for various ∆. We compare timings using
labelingcantakeoftheorderofn2 steps,asessentiallyeach Γ = n for FIFO and both Γ = n and Γ = ∞ for HPQ.
site mustbemovedstepwisetothemaximalheight. Whena We will referto the latter two choicesas HPQ andHPQ ,
n ∞
7
1000 100
L=25 L=8
100 L=50 FIFO, d = 2 L=16
L=100 L=32
L=200 L=64
10 L=128
L=256
L=512
t[s] 1 / nR 10 LL==12002448
NP L=4096
0.1
0.01
0.001
10 100 1000 10000 1e+05 1e+06 1e+07
Γ 1
0 1 2 3 4
∆
Figure 4: CPU time t[s] for HPQ with gap relabeling vs. Γ for
∆=2.2>∆x.Inthisregimegaprelabelingiseffectiveandglobal Figure5: [Coloronline]Numberofpush-relabeloperationspersite
updates are almost unnecessary. For Γ < n, the running time in- neededtosolveforthegroundstate, NPR/n, plottedasafunction
creasesduetofrequentglobalupdates. ForΓ>n,therunningtime of∆,fora2DsquarelatticewithGaussiandisorder. Thenumberof
staysconstant.Thereisnodiscernibleminimuminthetvs.Γcurve, spinsn = L2 rangesfrom82 to40962. Duetothelargerangeof
butΓ = ndoesnothurtperformanceatlarge∆andimprovesper- system sizesand hence running times, we haveused alogarithmic
formanceforsmaller∆. axisforplottingNPR/n. Statisticalerrorbarsandindividualpoints
are not displayed, as the point spacing is small and the statistical
uncertaintiesareverysmall,exceptforL=4096.Theglobalupdate
respectively. WecomputedthemeanvaluesoftandNPR for intervalusedwasΓ = n. PlateausareseenwhenNPR isclosetoa
alargerangeofsystemsizesanddisorders.Therunningtimes multipleofΓ,especiallyatlarger∆.Eachcurveshowapronounced
wereaveragedoveratleast100(sometimesasmanyas105) peak, whichcanbeusedtodefine∆x(L). Atthisvalue of∆,the
correlationlengthisoftheorder ofthesystemsizeandthesystem
samples,chosenindependentlyforeachvalueofLandΓ.
crossesoverfromaferromagneticregimewithonelargedomainto
Fig. 5 displays the running time per spin for FIFO with
aparamagnetic regimewherethecorrelationlengthissmallerthan
Γ = n and Fig. 6 presents a comparison of all three com-
thesystemsize. Thecrossoverfield∆x approacheszeroasLgoes
binationsforL=100. Thefirstnotablefeatureinthedatais
toinfinity,buttheapproachislogarithmicallyslow.
that all three combinationsshow a pronouncedpeak in N
PR
atavalueof∆ (L)between0.5and2.0overtherangeL=8
x
to4096.Thispeakisreminiscentofthecriticalslowingdown
near a phase transition. It moves slowly to the left with in-
creasing system size (Fig. 5), consistent with ∆x(L) decay- 1e+06 HHeeaapp,, ΓΓ==n ∞
ingwith increasingL. As notedin Sec. II, thereis nophase Queue, Γ=n
transitionin2D,butthereisacrossoverforanyfinitesizesys-
temfromaferromagneticregimeatlow∆toaparamagnetic
L=100
regime at large ∆ as the correlation length becomes smaller
R
thanthesystemsize. NP
1e+05
For small ∆, HPQ and FIFO perform very similarly.
n
HPQ ,ontheotherhand,is2to4timesslowerinthisregime.
∞
Near the crossover, the running times for HPQ and FIFO
n
start to deviate and the difference is largest at ∆ , where
x
HPQ is about 3 times slower than FIFO and HPQ is yet
n ∞
anotherfactorof3slowerthanHPQn. 100000.01 0.1 1 10 100
For ∆ > ∆x, FIFO starts to lose its advantageasdomain ∆
sizes become smaller and gap relabeling becomes more and
moreeffectivewithincreasing∆. When∆≈3,HPQ starts
n
tooutperformFIFO;when∆≈20,FIFOtakestwiceasmany Figure 6: [Color online] Number of PR steps NPR in d = 2 for
pushoperationsasHPQ . L = 100. Thisplot shows thetimingforHPQ with(Γ = n) and
n
without (Γ = ∞) global updates, compared with FIFO. Note the
significant improvement the global update gives to HPQ for field
strengths below and at the finite-size crossover.FIFO needs fewer
C. TimingsforthedoubleQueue(DFIFO)
push operations than HPQΓ=∞ over almost 3 orders of magnitude
in∆. Thisdifferencebecomesmorepronouncedasthesystemsize
Our method for simultaneously rearranging both positive increases.
and negative excess, DFIFO, performs well for small ∆ in
8
d = 2 but never better than FIFO. For large values of ∆, it 50
is much slower than FIFO and HPQ. DFIFO performs even L=8
worsein3D.Thereasonbecameclearwhenwevisualizedthe L=16
40 FIFO, d = 3 L=32
rearrangementsof excessfield. The visualizationshowsthat L=64
positive and negative excesses tend to miss each other. The L=128
L=256
potentiallandscape,givenbythetwoheightfields,isnotup- 30
n
dated when excess is pushed and is not consistentlyupdated / R
byrelabels. Therearrangementofexcessisgenerallyguided P
N
by the last global update: the negative excess is pushed to 20
where the positive excess was at the time of the last global
update,andviceversa. Thislackofcoordinationbetweenthe
10
twoheightfieldsbecomesmorepronouncedinhigherdimen-
sions,wheretherearemorepossiblepathsbetweensites. Itis
possiblethatusingsomeothercombinationoftheheightfields 0
2 3
fora heuristicwouldproducebetterresults. Onealternative, ∆
for example, would be to combine the separate height fields
forpositiveandnegativeexcessesintoasinglefield. Positive
excesswouldmovedownthe heightfield gradientandnega- Figure 7: [Color online] Push-relabel operations (cycles) per site,
tiveexcesseswouldmoveupthesamegradient. NPR/n, plottedvs. disorder strength ∆ind = 3, for L = 8, 16,
32,64,128,and256,forthepush-relabelalgorithmusingtheFIFO
datastructure. TheglobalupdateintervalisfixedatΓ = n. There
isaclearpeakintherunningtimenear∆ ≈ 2.3,whichisingood
D. Resultsforthreedimensions(d=3)
agreementwiththecriticalfield∆c = 2.270(5)found[14]forthe
ferro-toparamagneticphasetransitionforaGaussiandistribution.
Inthreedimensions,Γ=nremainsagoodchoiceforboth
FIFO and HPQ. As in the case d = 2, due to a very broad
120
minimuminthedependenceoftonΓ,theexactvalueofΓis
Priority queues
notcritical. Thisisinagreementwithpreviousresultswhere
d = 3
100
Γ = nhasbeenfoundtobeagoodchoiceforavarietyof
min
sparseanddensegraphs[13].
Thepeakinthetimingofthealgorithmisquitepronounced, 80 L=8, LPQ
L=16, LPQ
as in the 2D data, but in this case the location of the peak n L=32, LPQ
convergesto a fixed value at large L. We gathered data for / R 60 L=64, LPQ
thenumberofpush-relabeloperationsN persite,N /n, NP L=128, LPQ
PR PR L=8, HPQ
vs. ∆, for sizes L ranging from 8 to 128 for HPQ and LPQ 40 L=16, HPQ
anduptosize256fortheFIFOdatastructure.Theplotofthe L=32, HPQ
FIFOdata(Fig.7)showsagrowthinN /nwithLandcon- L=64, HPQ
PR 20 L=128, HPQ
vergenceto a well-defined curve for the running time for ∆
intheparamagneticrange.Intheferromagneticregime,there
0
is a slow growthof the runningtime with sample dimension 2 3
∆
L. Both the LPQ and HPQ data structures are significantly
slowerthantheFIFOdatastructure,atmoderatevaluesof∆
(witharatioof≈3forthepeakrunningtimesatL=128),as
Figure 8: [Color online] Push-relabel operations (cycles) per site
seen in Fig. 8. The LPQ data structureis significantlyfaster
NPR/n vs. ∆ in d = 3, for L = 8, 16, 32, 64, 128, and 256,
thantheHPQ structureonthe paramagneticsideofthe peak
whenusingtheHPQandLPQdatastructures.Thereisagainaclear
in the running time, but is very similar in speed on the fer- peakintherunningtimenear∆≈2.3,butthepeakissignificantly
romagneticside. Therunningtime forLPQ convergesmuch broaderandhigherforHPQ,thanforLPQ.Notethatthenumberof
morequicklythantheHPQ versiontoanL-dependentvalue cyclesgrowsmorequicklywithLthanfortheFIFOstructure.
asLisincreasedatfixed∆>∆ .
c
i.e., where r = 0. A sample snapshotfrom the simulation
ij
V. QUALITATIVEDESCRIPTION withbothoftheseoptionsactivatedisshowninFig.9.
When∆issomewhatlargerthan∆ ind = 2, the differ-
x
In order to better understand the timing results, we have encesbetweenthetemporalprogressusingtheHPQandFIFO
visualizedtheevolutionoftheheightfieldsandtherearrange- data structuresare clearly seen in the dynamicvisualization.
mentofexcess. Thevisualizationcode[26]usesacolormap As FIFO cycles through all active sites, all the positive ex-
todisplaytheheightfield. Theprogramhastheoptiontodis- cessesmoveataroughlyuniformspeeddowntheheightgra-
play the location of sites with excess(white for positiveand dients. TheutilityoftheglobalupdateatlatetimesforFIFO
blackfornegative)andtoindicatewheretheflowissaturated, isapparent:whenΓislarge(infrequentglobalupdates),afew
9
For samples with weak disorder, the algorithm with the
HPQ data structure doeslead to a few positive excessesrac-
ing around the network while the others sit idly. The time
to raise a large positive domain high enoughto create a gap
is large. Again, active sites tend to movearoundquite a bit,
while other regions remain unchanged. Towards the middle
of the algorithmexecution,the lattice oftendisplayslines or
rings of positive excess sitting on an equipotential line near
asink,oftenimmediatelyadjacenttoit. Thesestructuresare
at low height and cannot move until all excesses of greater
height have been removed. This should slow the algorithm
down. Thepositiveexcessesdonot“screen”thesinks, since
the global update doesn’t differentiate between sinks of dif-
ferentcapacity.Thisexpectationisconsistentwithourresults
fortheLPQ,whichperformssomewhatbetterthanHPQ,es-
peciallyathigherfieldstrengths.
Atlargefieldstrengths,theproportionofactivesiteswhose
random field strength is greater than the strength of their
bonds to neighbors (i.e., ∆ > 4 in 2D) is large. The algo-
rithmactsonmostsitesonlytwice. Onecanclearlyseehow
HPQsweepsthelatticeandraises(nearly)allthesitesofini-
tialheighttwo(thosewithnoadjacentsinks)tothemaximum
height, then doesthe same to sites of heightone. FIFO also
Figure 9: [Color online] Image of the state of the auxiliary fields generallymakestwoorthreepasses, thoughitsfirstpassex-
duringthesolutionofa2DRFIMofsizeL=50andwithdisorder ecutespushesonthesiteswherepushingispossibleandrela-
strength∆=2.Thispictureshowsasnapshotoftheheightfieldui, belsotherwise.
the location of non-zero excess e(i), and bonds with zero residual
strength(rij =0),whileexecutingthepush-relabelalgorithm(HPQ
data structure). The red (darkest non-black) regions have already
VI. SUMMARY
been identified as up; the heights in this region have the maximal
valuen = 2500. Thecolorsthatarenotblackorwhitecorrespond
totheheightfieldintheremainderofthesample:blue(darkercolors) Ashasbeendemonstratedbefore[11,14,29,30,31,32,33,
indicatesalow,whilegreen(lightercolors)indicatesahighheight. 34],thePRalgorithmisanefficientmethodtofindtheexact
The white squares represent spins with a positive excess field (the ground state for the RFIM at T=0 in any dimension. Here
heightatthesesitesisnotindicated); blacksquaresrepresent spins weinvestigatedhowtoimplementPRtofindthegroundstate
withanegative excess (and haveheight zero). Theblack linesare
mostquickly. We compareda numberofdatastructuresand
dualtothebondsthataresaturated(i.e.,havezeroresidualstrength).
theeffectofglobalandgaprelabelingontheperformanceof
Thesesaturateddualbonds(rij = 0orrji = 0)“block”therear- thealgorithm.
rangementofexcessei. Moreup-spinregionsmaybeidentifiedas
Inagreementwithpreviouswork[34], ourdetailedresults
morebondsaresaturatedandthedualbondslinkuptoisolateare-
gion. Thealgorithmterminateswhenallpositiveexcessisconfined recommend the FIFO-queue combined with global updates
toup-spinregions. every Γ ≈ kn steps, with k a number near unity. FIFO
performs much better near the crossover disorder (∆ for
x
d = 1,2) and the critical disorder (∆ in d = 3) and never
c
positiveexcessesareseentoskatearoundinregionscontained performssignificantlyworsethanHPQ.TheexactvalueofΓ
bysaturatedbonds. Thepusheshavefoundtheminimumcut isnotcrucialsincetheminimumintvs.Γisverybroad.
bysaturatingthebondsthatseparatethepositiveandnegative We also tried an implementation that treats positive and
spinregions,butthealgorithmhasn’tconfirmedthatfactyet. negativeexcessesequally.Thisimplementation,however,suf-
Thelocalrelabels,whichtendtoraiseasite’sheightonlyby fersfromthelackofcoordinationbetweentheheightfieldsfor
one step at a time, are inefficient in raising a large up-spin the two sets of excesses and positive and negative excesses
domainto maximalheight. The remainingpositive excesses tendto miss eachother. Itis likely thatthis algorithmcould
areshuffledaroundwithinthedomain,slowlyincreasingthe be improved by using a single height field to coordinate the
heightbysmallrelabels. Whenglobalupdatesareinfrequent, motionofthepositiveandnegativeexcess.
the algorithm may terminate by a final globalupdate, which A more flexible approach to the global update interval
findsthatasetofpositiveexcessesisisolatedfromallsinks. might also be useful in speeding up simulations. It is likely
This confirmsthe picture discussed in Sec IVA2. HPQ, on thatadaptivelymodifyingtheglobalupdatessothattheyare
theotherhand,tendsto actrepeatedlyonthesame site. Iso- executedwhenthesinkdensitychangesbyadefinedfraction
latedregionswithpositiveexcessareraiseduniformlyabove orpacketsofexcesshavetravelledagivendistancewouldop-
therestofthesample.Thisallowsthegapheuristictoquickly timizethealgorithmduringeachstageofthesolution. There
identifyisolatedpositivespindomains,evenwhenΓislarge. is still a lot of room for other modification, e.g., cutting off
10
the breadthfirst search to reflectsaturation and non-uniform Acknowledgments
intervalsforglobalupdatethatdependonthenumberofsinks
thathavebeenannihilatedsince thelast globalupdaterather
thanthenumberofPRsteps.
We foundthatthevisualizationoftheoperations(see[26]
for the source code) greatly improved our understanding of This work has been supported by the National Science
thealgorithm. Thiscodemaybeusefulinsuggestingfurther Foundation under grants ITR DMR-0219292 and DMR-
improvementstothealgorithm. 0109164.
[1] A.P.Young,SpinGlassesandRandomFields,vol.12ofSeries [19] E.T.SeppäläandM.J.Alava,PhysicalReviewE63,066109
onDirectionsinCondensedMatterPhysics(WorldScientific, (2001).
1998). [20] M.AlavaandH.Rieger,PhysicalReviewE58,4284(1998).
[2] S.FishmanandA.Aharony,J.ofPhys.C:SolidStatePhys.12, [21] T.H.Cormen, C.E.Leiserson,andR.L.Rivest,Introduction
L729(1979). ToAlgorithms(MITPress,Cambridge,Massachusetts,1990).
[3] D. P. Belanger, Spin Glasses and Random Fields(World Sci- [22] M.J.Alava,P.M.Duxbury,C.F.Moukarzel,andH.Rieger,in
entific, 1998), vol. 12 of Series on Directions in Condensed PhaseTransitionsandCriticalPhenomena,editedbyC.Domb
MatterPhysics,chap.ExperimentsontheRandomFieldIsing andJ.L.Lebowitz(AcademicPress,2001),vol.18.
Model,pp.251–275. [23] A. K. Hartmann and H. Rieger, Optimization Problems in
[4] J.Villain,PhysicalReviewLetters52,1543(1984). Physics(Wiley-VCH,2002).
[5] D.S.Fisher,PRL56,416(1986). [24] A. A. Middleton, New Optimization Algorithms in Physics
[6] S.Kirkpatrick,C.D.Gelatt,andM.P.Vecchi,Science220,671 (Wiley-VCH,2004), chap. CountingStatesandCountingOp-
(1983). erations,pp.71–100.
[7] J.C.Anglesd’Auriac,M.Preissmann,andR.Rammal,Journal [25] A. W. Goldberg, Andrew gold-
dePhysique-Lettres46,173(1985). berg’s network optimization library,
[8] SorinIstrail,inProceedingsofthethirty-secondannual ACM AndrewGoldberg’sNetworkOptimizationLibrary, URL
symposium on Theory of computing (ACM Press, 1999), pp. http://www.avglab.com/andrew/soft.html.
87–96. [26] D. C. Hambrick, RFIM applet, RFIMApplet (2003), URL
[9] J.C.PicardandH.D.Ratliff,Networks5,357(1975). http://www.phy.syr.edu/research/cmt/RFIM/.
[10] F.Barahona,J.Phys.A18,L673(1985). [27] B. V. Cherkassky and A. V. Goldberg, Algorithmica 19, 390
[11] A.T.Ogielski,PRL57,1251(1986). (1997).
[12] A.GoldbergandR.Tarjan,JournaloftheAssociationforCom- [28] A. V. Goldberg and R. Kennedy, SIAM Journal of Discrete
putingMachinery35,921(1988). Mathematics10,390(1997).
[13] B. V. Cherkassky and A. V. Goldberg, Algorithmica 19, 390 [29] S.BasteaandP.M.Duxbury,cond-mat/9801108v3(1998).
(1997). [30] P.M.DuxburyandJ.H.Meinke,PhysicalReviewE(2001).
[14] A.A.MiddletonandD.S.Fisher,PRB65,134411(2002). [31] J. Esser, U. Nowak, and K. Usadel, Physical Review B 55
[15] A.A.Middleton,PRL88,017202(2002). (1997).
[16] G.Schröder,T.Knetter,M.J.Alava,andH.Rieger,Eur.Phys. [32] A.Hartmann,Ph.D.thesis,Ruprects-Karls-Universität,Heidel-
B24,101(2002). berg(1998).
[17] K.Binder,ZeitschriftfürPhysikB-CondensedMatter50,343 [33] A.K.HartmannandU.Nowak,cond-mat/9807131(1998).
(1983). [34] E.Seppälä,Master’sthesis,HelsinkiUniversityofTechnology,
[18] G. Grinstein and S. k. Ma, Physical Review Letters 49, 685 FIN-02150Espoo,Finland(1996).
(1982).
