Belief propagation for graph partitioning

Belief propagation for graph partitioning P etr ˇ Sulc 1 , 2 , 3 , Lenk a Zdeb oro v´ a 1 1 Theoretical Division and Cen ter for Nonlinear Studies, Los Alamos Nation al Laboratory , Los Alamos, NM 87545, USA 2 New Mexico Consortium, Los Alamos, NM 87544, USA 3 F acult y of Nuclear Sciences and Physical Engineering, Czech T ec hnical Unive rsi t y , B ˇ reho v´ a 7, CZ - 115 19 Prague, Czec h Republic E-mail: sulcpetr@ gmail.com, lenka.zd eborova@gmai l.com Abstract. W e study the beli ef propagation algorithm for the graph bi- partitioning problem, i .e. the ground state of the f erromagnetic Ising mo del at a fixed m agnet ization. Application of a message passing s c heme to a mo del with a fixed global parameter i s not banal and we show that the magnetization can i n fact b e fixed i n a lo cal wa y w i thin the b elief propagation equations. Our metho d prov ides the full phase di agram of the bi-partitioning pr oblem on random graphs, as we ll as an efficien t heuristic solver that we ant icipate to b e useful i n a wide range of application of the partitioning problem. P A CS num b ers: 75.10.Nr, 05.70.Fh, 05.70.Ce, 02.70.-c Keywords: g raph partitioning, belie f pro pa gation, global co nstraint, ra ndom g raphs, graph bisec tio n, r eplica symmetry bre a king 1. In tro duction Graph par titio ning problem was one of the firs t optimization problems treated with methods of statistical mechanics o f diso rdered systems [1, 2]. Since then other applications of the theo ry of spin gla sses in har d optimizatio n and cons train t satisfaction problems attracted a lot of int eres t and ma n y re mark able res ults were obtained. As anticipated in the early works [2, 3], understa nding of the e nergy landscap e and the phase transitions in the s pace o f solutions leads to under standing of algorithmic har dness of the problems [4, 5], and even more r emark ably it leads to a developmen t of a new class of heuristic a lgorithmic techniques [4]. Now a da ys, the cavit y metho d [6] ser v es as a sta te of art technique for understanding random optimization problems, a nd its applicatio n on a given ins tances of the pro ble m is a base for a class o f one o f the mos t promising he ur istic so lv ers, kno w n a s mes sage passing a lg orithms in computer sc ience. Despite all this activ ity in the field, nor the phase diagr a m neither a message passing algo rithm for pa rtitioning a gr aph into t wo gro ups of a given s ize has b een work ed o ut. The main reas on that makes the graph partitioning a tricky pro blem to tr e a t is the ex istence of a global constraint that fixes the size of the tw o gr oups. The aim of this ar ticle is to fill this gap, and give the phase diagram of the g r aph bi-partitioning on spar se random graphs and asso ciated b elief propaga tio n algo rithm. Belief pr op agation for gr aph p artitioning 2 1.1. Partitioning pr oblem: Setting and applic ations A gr aph G ( V , E ) is given by the set o f vertices V a nd e dges E . If an element ( i, j ) belo ngs to the s et of edges we say that vertices i a nd j a re connected. The g raph bi-partitioning pro blem consis ts of dividing vertices o f the gr aph into tw o disjoint sets of a given siz e , so as to minimize the num b er of connectio ns b et ween vertices from different gr oups. The problem is known to b e NP-complete [7], a nd hence there is a go o d r eason to b elieve that no exa ct p olynomia l a lgorithm exists. The gr a ph partitioning pr oblem is encountered in ma n y imp ortant applications. T o give few examples: In an electric circuit design one needs to know o n which b oard to place the different c omponents to minimize the num b er of links be tw een different bo ards [8]. In pa rallel computing o ne has to partition data and tasks among several pro cessors in or der to minimize the communication b etw e en them [9]. Partitioning is also closely related to da ta clustering and communit y detection [1 0]. The list could contin ue fo r lo ng, and it is hence crucia l to develop efficient heuristic a lgorithms that give g oo d solutio ns to the pr oblem. A lar ge volume of litera ture on heuristic metho ds for gr aph par titioning exists. One of the ear ly fundamental works in the field is [1 1], its r unning time is, how ever, O ( N 2 ) so it is no longer used in prac tice. Simulated annealing techniques can b e used, see e.g . [12, 13 ]. A lo cal search based metho ds such as the extremal optimiza tion of [14] were suggested. There is a whole class of sp ectral partitioning metho ds that use the eigenv ectors of the Laplacian of the co nnec tiv ity graph, see e.g. [1 5]. How ever, the current state of art metho d for partitioning, that is used in most practica l applications, is based on the multi-lev el progra mming: The no des a re gr ouped into sup er-no des and the sup er-no des group ed again, at the end the system size is very s mall and the problem is so lv ed exactly and the gro uping o f no des is then unwrapped. The multi- level pr ograms use elements from many other appro ac hes, see [1 6] fo r an exc ellen t review. W e do no t anticipate that b elief propa gation de velop ed in this pap er, will be by itself comp etitiv e with the highly tuned implementations of the multi-level metho ds. How ever, we do anticipate that it ca n b e used as a comp onent of these implemen tations. F or example, in the multi-lev el algo rithms o ne nee ds to estimate the pro babilit y that t wo no des ca n be gr ouped in the same sup er-no de — this is exactly what b elief propag ation is designed to compute very fast and efficiently . The gra ph bi-partitioning pr oblem is equiv a len t to finding the ground state of the Ising mo del with fixed magnetiza tion. The energy in the Is ing mo del is given by the following Ha miltonian: H = − X ( ij ) ∈ E S i S j , (1) where S i is the Ising spin (either +1 of − 1) on the i -th vertex of the gra ph. The magnetization m , − 1 ≤ m ≤ 1, is given by 1 N X i S i = m , (2) where N is the num b er of vertices. Therefore , the pr oblem of finding a configuration of spins that minimizes (6) while demanding magnetiza tio n m to b e fix ed is equiv a len t to dividing vertices into tw o gr o ups of size N (1 + m ) / 2 a nd N (1 − m ) / 2 such that the num b er of links b et ween them is minimal. F or m = 0, the g raph is divided in to t wo gr o ups of equal size, i.e. the g raph bisection. The co s t of a g raph partitioning Belief pr op agation for gr aph p artitioning 3 at a given mag netization, that we call b ( m ), is g iven as the num b er of edge s b et ween different groups divided b y the tota l num b er of vertices. The r elation b et ween b ( m ) and gr ound sta te energ y E ( m ) of the Ising mo del at mag ne tiza tion fixed to m is b ( m ) = E ( m ) + M 2 N , (3) where N is the num b er of no des, and M the num b er of edges. 1.2. Pr evious r esult s on bi-p artitioning ra ndom gr aphs In gr aph theory es timating the asymptotic size o f the bisection width in r andom regular gr aphs, i.e. graphs of a fixed degree chosen uniformly at rando m fro m all the po ssible ones, is a classical question. Many upp er and low er b o unds w ere derived. The currently b e st known upp er and low er bo unds on bisection width in random regular graphs are by [17, 18, 19, 20, 2 1] and we summariz e their numerical v alues in T a ble 1 and Fig. 4. F or Erd˝ os-R´ enyi ra ndom graphs with N → ∞ vertices a nd mean degree α (every edge is present with probability α/ ( N − 1)), the siz e o f the larges t comp onent is g N + o ( N ), where g sa tisfies the following equation: g = 1 − e − αg . (4) In or der to div ide the gr aph into t wo parts of size N (1 + m ) / 2 and N (1 − m/ 2) such that the num b er of edges betw een the t wo is zero, the size o f the largest comp onen t g m ust be at max im um (1 + m ) / 2. That is p ossible for average degree α < α s where α s = − 2 1 + m log 1 − m 2 . (5) F or α > α s an extensive num b er of edg es needs to b e c ut in the minimal bipartitio n. The v a lue α s is hence in a sense the sa tisfiabilit y threshold fo r gra ph par titioning o f Erd˝ os-R´ en yi rando m g raphs. This is further discussed in [22 ], where the a uthors a lso obtain a n interesting upp er b ound on the bisec tion width ( m = 0). In statistical physics many articles a ddressed the ra ndom graph bi-partitioning problem, see e.g . [2, 12, 23, 24, 2 5, 26, 27, 2 8], but as far as we can tell they a ddr ess only cases wher e (A) the ma gnetization is fixed to z e ro, (B) the fluctuations in the degree of the ra ndo m graph are negligible , i.e. the graphs are either dense or regular . The computational techniques used in the ab ov e mentioned pap ers do not g e ne r alize to the non- zero magnetization case nor to gra phs with fluctuating degree, as e.g. to the E rd˝ os-R´ enyi ra ndom gr aphs. W e will give a mo re detailed ex planation of why the techn iques do not generalize in section 4.3 . This also justifies nov elt y of the a pproach developed in this ar ticle. 1.3. Contribution of t his article If the gr ound state energy of the ferro magnetic Ising mo del (1) was a con vex function of the mag netization m then a n external magnetic field (playing the ro le of the chemical po ten tial from the gra nd-canonical ensemble) co uld b e used to compute E ( m ) with a standard cavit y metho d [2 9]. How ever, r andom graphs are mean field top olog ies and the energ y at fixed magnetiza tion E ( m ) do es not hav e to b e and in this case is no t a conv ex function, similarly as in the fully co nnected Curie-W eiss mo del. The pr oblem of imp osing the v alue of the ma gnetization is hence more challenging. Belief pr op agation for gr aph p artitioning 4 A metho d to explor e the non-conv ex pa rts o f thermo dynamical p otentials within the Bethe-Peierls (Belief P ropagation) approximation was suggested in [3 0], and used later e.g. in [31, 3 2]. The main idea is to introduce a n unifor m externa l magnetic field (or c hemical potential) a nd adjust its v alue after every up date of the lo cal cavit y fields. W e use this metho d for partitioning graphs, and we arg ue that it (or its generaliz a tion to the replica s y mmetry breaking scheme) is asymptotically ex act on s parse random graphs. The ma in pra c tica l contribution of this article is the b elief propa gation algor ithm for graph pa r titioning pr oblem that w e b elieve to b e of use in the v ar ious applications of the pro blem. W e study the b eha vior a nd per formance of the algor ithm on ra ndom graphs but we anticipate it will b e meaning ful and useful fo r o ther families of graphs, complex netw orks for e x ample. W e also c o mpute the phase diagram of (1) a t fix ed magnetization. In [22] it was argued that in the E rd˝ os-R´ enyi gr aphs at zero magnetizatio n the glassy trans ition happ ens at some av erage degree strictly larger than the satisfia bilit y thresho ld, α c > α s , we indeed confirm this conjecture , we compute α c and several other quantities of interest. An interesting side re ma rk, discussed in s ection 4.3, co ncerns the cas e treated in the prev io us works: the reg ula r random graphs at zero magnetiza tio n. There the av erage prop erties of the graph bi-partitioning a re equiv alent to those o f the spin g la ss problem. W e argue why this equiv a lence do es not g eneralize to non-zero magnetization or no n-regular gra phs. More detaile d discuss ion a bout the equiv alence can b e fo und in [33]. 2. Ca vi t y metho d at fixed magneti zation As w e expla ined in the in tro duction, the graph partitioning is equiv a len t to the ferromag netic Ising mo del at fixed magnetization m ∗ . The magnetizatio n will be fixed via an externa l magnetic field h which app e ars in the Ha miltonian as H h = − X ( ij ) ∈ E S i S j − h X i S i . (6) The g round state energy density of (1) and (6) a r e related via the Legendre transformatio n e ( h ) = e ( m ) − hm , so that the parameter h has to b e chosen such that ∂ e ( h ) ∂ h    h ∗ = − m ∗ . (7) If e ( h ) is the gr ound energy density of (6) w ith field h cor resp onding to magnetization m , the corres ponding par tition cost (3) of the gr aph is b = e ( h ) + hm + α 2 2 , (8) where α is the mean degr e e of the g raph. 2.1. Belief pr op agation e quations The Be the- P eierls a ppr o ximation, or the Belief-Propag ation equa tions, aims to describ e the Boltzma nn measure of (6) µ ( { S i } ) = e − β H h ( { S i } ) Z , (9) Belief pr op agation for gr aph p artitioning 5 where β is the inverse tempe r ature. The graph pa rtitioning problem corres p onds to β → ∞ . In this section we summar ize the well known be lie f propaga tion equa tions for this problem. F or a detailed deriv a tion see [34, 3 5]. In the mo s t standa rd form o f b elief propaga tion e quations [34] one introduces ψ i → j S i to b e the proba bilit y that v ariable i ta k es v alue S i given the interaction on ( ij ) is abs en t. On a tree (cycle fr e e) graph then ψ i → j S i = 1 Z i → j e β hS i Y k ∈ ∂ i \ j  X S k e β S i S k ψ k → i S k  , (10) where Z i → j is normalization ensuring ψ i → j +1 + ψ i → j − 1 = 1. After a fix ed p oint of eqs. (10) is found the Bethe fr ee energy (or the log -partition function) is given as [34] − β F ( h ) = log Z = X i log Z i − X ( ij ) log Z ij , (11) where Z i = X S i e β hS i Y k ∈ ∂ i  X S k e β S i S k ψ k → i S k  , (12) Z ij = X S i ,S j e β S i S j ψ i → j S i ψ j → i S j . (13) A t a given v alue o f the exter nal magnetic field h the av era ge magnetization is computed as m = − [ ∂ F ( h ) /∂ h ] / N , using (11) one gets m = 1 N X i P S i S i e β hS i Q k ∈ ∂ i  P S k e β S i S k ψ k → i S k  P S i e β hS i Q k ∈ ∂ i  P S k e β S i S k ψ k → i S k  . (14) In order to write the zero temp erature limit, β → ∞ , of the ab ov e equations we int ro duce mor e suitable mes sages (usually called cavity fields) h i → j e 2 β h i → j ≡ ψ i → j +1 ψ i → j − 1 (15) One then obtains eq uations equiv a len t to the replica symmetric equations in [29]. The self-consistent e quations for messa ges (10) b ecome h i → j = h + X k ∈ ∂ i \ j  max (1 + h k → i , 0) − max ( h k → i , 1)  ≡ F ( { h k → i } ) . (16) Note that the term in the s um is − 1 for h k → i ≤ − 1, + 1 for h k → i ≥ 1, and h k → i for − 1 < h k → i < 1. The Bethe estimate of the gr o und s tate energy is E ( h ) = X i E i − X ( ij ) E ij , (17) where from (12-13) we obtain E i = h + α + 2 X k ∈ i max (0 , h k → i ) − 2 max [ h + X k ∈ i max (1 + h k → i , 0) , X k ∈ i max ( h k → i , 1)] (18) E ij = 1 + 2 ma x (0 , h i → j ) + 2 max (0 , h j → i ) − 2 max (1 + h i → j + h j → i , 1 , h j → i , h i → j ) . (19) Belief pr op agation for gr aph p artitioning 6 And finally the contribution to the sum P i in the express ion for the magnetization (14) is in the zero temp erature limit equal to +1 if h + X k ∈ ∂ i  max (1 + h k → i , 0) − max ( h k → i , 1)  > 0 , (20) and − 1 otherwise (if the t wo terms ar e equal the contribution is 0). F o r notatio n let us ca ll this function m = M ( h, { h i → j } ) . (21) 2.2. Population dynamics at fixe d magnetization In order to calculate the average ground state ener gy (17), and thus the partitioning cost b , for a given ense m ble of r andom graphs one implements the popula tio n dynamics metho d [6, 3 5]. In the standard p opulation dynamics one would up date eq s. (1 6) with a given v alue o f external ma gnetic field h till conv er gence or till maximum num b er of iterations and then o ne w ould compute the g round energy and the co r resp onding ma gnetization. If this is done with the a bov e equations for gra ph bi-par titioning then the res ulting magnetization will alwa ys b e either +1 for h > 0 or − 1 for h < 0. W e how ever wan t to find the ground s tate energy at magnetization v alues − 1 < m ∗ < 1. In or der to do that we will not keep the external field h cons tan t. Instead after every itera tion of (16) we use the cur r en t v alued of fields h i → j and update the v alue of h in such a wa y that m ∗ = M ( h new , { h i → j } ) where m ∗ is the desir e d v a lue of the magne tiza tion. The resulting p opulation dynamics co de is sketc hed in algo r ithm 1. Note that the alg orithm 1 uses bisection metho d in each itera tion in orde r to fix the mag ne tiza tion. F or d - r egular graphs , M ( h, { h i → j } ) as a function of h for given v alues of { h i → j } is contin uous monotonic function and therefor e the a lgorithm always manages to fix the desir ed magnetization. F or general g r aphs M ( h, { h i → j } ) may hav e less well b eha ved form and we will discuss that in s ection 3. 2.3. 1RSB at fi xe d magnetization As may be a nticipated from the relation b et ween gra ph bisection and the spin g lass [2] the b elief pr opagation equa tions (replica symmetric appr oach) a r e not alwa ys asymptotically exac t for the g raph bi-pa rtitioning. Instead in so me regions of parameters the problem is glassy a nd the replica symmetry brea king a pproach is needed for a n exact solution, just like in the Sher rington-Kirk pa tric k mo del [36, 3 7]. The replic a symmetry br e aking a pproach for spa rse random gra phs w as developed in [6, 29] and is w ell established. Hence we o nly p oint out the difference in the equations that leads to fixing a non-triv ial v alue of the mag netization. In order to write the 1 RSB equations we follow clo sely the approach o f [29]. W e int ro duce a complexity function Σ( E ), i.e. n umber o f thermo dynamical s ta tes at a given ener gy , a nd its Legendr e transfo rm Φ( y ) also calle d the replica ted p otential − y Φ( y ) = − y E + Σ( E ) , ∂ y Φ( y ) ∂ y = E . (22) Every thermo dynamical sta te has a cor resp o nding v a lue o f the cavit y field h i → j and according to [29] the self-co nsisten t equation for the distribution o f cavit y fields over Belief pr op agation for gr aph p artitioning 7 Algorithm 1 Population dyna mics algor ithm for BP on d -r e g ular random g r aphs with fixed magnetization m ∗ h ← 0 Randomly initialize message s h ( i ), i = 1 , 2 . . . N for j = 1 to max do for i = 1 to N do Randomly select d indices in k = 1 , 2 . . . N Calculate h ( i ) from { h ( k ) } us ing eq. (16) end for h 1 ← h − 1 h 2 ← h + 1 while | h 1 − h 2 | < criterion do Calculate mag netiza tion m with external field h using eq . 21 if m < m ∗ then h 1 ← h end if if m > m ∗ then h 2 ← h end if h ← ( h 1 + h 2 ) / 2 end whi le end for Calculate E using (17) (av erag ed ov er so me num b er o f itera tions) return E , h states is P i → j ( h i → j ) = 1 Z i → j Z Y k ∈ ∂ i \ j d P k → i ( h k → i ) e − y E i → j δ [ h i → j − F ( { h k → i } )] , (23) where F ( { h k → i } ) is defined by e q . (1 6). The reweigh ting factor is defined by E i → j = − lim β →∞ 1 β log Z i → j where Z i → j is the normalizatio n cons tan t in (10) a nd is g iven by an equation ana lo gous to (18). Once a fixed po in t of (23) is found the po ten tial Φ( y ) is computed as follows Φ( y ) = P i Φ i − P ij Φ ij with e − y Φ i = Z P OP e − y E i , e − y Φ ij = Z P OP e − y E ij , (24) where the notation R P OP = R Q k ∈ ∂ i \ j d P k → i ( h k → i ) and the energy contributions ar e given by (1 8-19). The ener gy of the system is then computed accor ding to (22 ) a s E = P E i − P ij E ij with E i = R P OP E i e − y E i R P OP e − y E i , E ij = R P OP E ij e − y E ij R P OP e − y E ij . (2 5) And the ma gnetization m = P i m i / N , wher e m i = − ∂ E i ∂ h = − R P OP ∂ E i ∂ h e − y E i R P OP e − y E i + y R P OP ∂ E i ∂ h E i e − y E i R P OP e − y E i − y R P OP E i e − y E i R P OP ∂ E i ∂ h e − y E i ( R P OP e − y E i ) 2 . (26) Belief pr op agation for gr aph p artitioning 8 Note that ∂ E i ∂ h = ± 1 dep ending on the sign in eq. (21). Again the only difference betw ee n the usual 1RSB and 1RSB at fixed magnetization is that after every itera tion the externa l magnetic field is chosen a new v a lue such that mag netization computed from (26) is equal to the desired v a lue m ∗ . Solving the 1 RSB equa tions is often tedio us and to obtain the phase diagr am it is often sufficient to investigate the conv ergenc e of the b elief propaga tion iterations. This is equiv a len t to analy zing the lo cal stability of the replica symmetric solution tow a rds r eplica symmetry br eaking, a s done or ig inally b y de Almeida and Thoules s [38]. Within the po pulation dynamics we use the noise propa gation metho d (for a deriv ation see app endix C of [39 ]). In the p opulation dy na mics algorithm to gether with cavit y fields h i → j , one keeps track o f the qua ntit y v i → j = X k ∈ ∂ i \ j ∂ h i → j ∂ h k → i v k → i (27) after every sweep of B P iteratio n we no r malize the v alues v i → j by λ in such a wa y that P ( v i → j /λ ) 2 = 1. Parameter λ then plays a role of a certain Lyapunov exp onent and the b elief pr opagation do es no t conv erge if and only if o n av erage λ > 1. W e hav e found that BP never conv erges o n r egular gra phs for a n y v a lue of magnetizatio n − 1 < m ∗ < 1. Nevertheless, the v alue of the ener g y calculated with B P gives a go o d low e r bo und on the a ctual energy o f the mo del [40]. F or the Erd˝ o s -R ´ e nyi gr aphs with given magnetization, we found a pha se transition fro m a re plica sy mmetr ic region where BP is asymptotically exact to a glassy region where RSB solution would be need to obtain the asymptotically exact v alue of the gro und state energy (this phase transition is shown in Fig. 1). 3. BP as a heuristi c solver Equations for the b elief pr opagation derived in the previous sectio n can b e used on a given g r aph a s we sketc h in algorithm 2. The parameter memory , which we set to 0 . 7 in our simulations, is intro duced in order to pr ev ent mess ages fro m oscilla ting. If the algor ithm do es not conv erge a fter a given maximum num b er of iteratio ns, it is terminated. Howev er , ev en if the alg orithm do es not conv er ge, the ca lculated E still provides a reasona ble estimate of the bisection cost that is on av e r age a low er b ound of the true av erag e cost. In the presented algorithm, we introduced a s ligh tly differe n t metho d to fix the magnetization by manipulating h . In the algor ithm 2, we sort all the lo cal cavit y fields and se t h so that N (1 − m ) / 2 of them are nega tiv e (or zer o ) a nd the rest p ositive (or zero ). It follows from the definition of messag es (15) tha t the p ositiv e v alue of lo cal c avity field mea ns tha t s pin o n this no de is more likely to b e equa l to 1, nega tiv e means that the spin is more likely to b e − 1. If the lo cal cavit y field is exactly equa l to zero, the spin in a given no de is unbiased (free). This can b e use d to actually obtain a graph partition. Howev er a decimation technique, alg o rithm 3, achieves m uch b e tter results in particular when man y free or almost free spin a re pr esen t, rep orted in Figs. 3 and 4. Note that in the present form the decimation so lv er has running time quadratic in the s ize of the system. Ho wev er, linear r unning time ca n b e achiev ed without significant lo ss of p e r formance by decimating a finite fra ction of s pins after every iteration, as in the survey pro pagation algor ithm of [4]. Belief pr op agation for gr aph p artitioning 9 Algorithm 2 BP algor ithm for partitioning o f a given g raph ∀ i, j Initializ e messages h i → j and field h rando mly iter ← 0 rep eat for all i ∈ V do conv ergence ← 0 lo cal field[i] ← P k ∈ ∂ i  max (1 + h k → i , 0) − max ( h k → i , 1)  h i → j new ← h + lo cal field(i) - h j → i conv ergence ← conv erg ence +   h i → j new − h i → j   h i → j new ← memory ∗ h j → i + (1 − memory) ∗ h i → j new end for sort(lo cal field) h = - lo cal field[ N (1 − m ) / 2 ] iter ← iter +1 un til conv ergenc e < ǫ OR iter > maximum iterations compute E using equa tion (17) return E , h Algorithm 3 Decimation algo rithm rep eat Run algo rithm 2 Cho ose a vertex i such that lo cal field( i ) is the hig hest (or low est if this is a n even iteration) and fix a ll outgoing message s fro m this no de to + ∞ ( - ∞ fo r an even iteration). Fix spin in vertex i to +1 ( − 1 in even itera tio n) un til Num b er of fixed spins to +1 or − 1 reaches the v alue r equired to fix de s ired magnetization m 4. Beha vior of the metho d and results In this sectio n we dis cuss the be havior of the b elief pro pagation a lgorithm o f r andom regular and Erd˝ o s-R ´ enyi random g raphs. W e, howev er, a nticipate that qualitatively similar b ehavior as on the Erd˝ os-R´ enyi rando m graphs will b e seen on o ther gr a ph families. 4.1. Phase diagr am of Er d˝ os-R ´ enyi gr aphs bi-p artitioning The mo st interesting to discuss is the b ehavior of the algo rithm on a given gr aph and the decima tion. In particular : Is the function M ( h, { h i → j } ) (21) contin uous in h such that any v alue o f the magnetization can b e fixed? Do es the exter nal field h converge in the iterations ? Do es the decimation a c hieve low ener gy states ? W e choose t ypical Erd˝ os -R ´ enyi rando m gra phs to illustrate the b eha vior and answer these questions . An Erd˝ o s-R ´ enyi random gr aph of av erag e degr ee sma ller that o ne, α < 1, basically lo oks lik e a collection of small disconnected trees. Let us hence first discuss how do es the algo rithm b ehav e on a tr ee. O n a tree the b elief propag a tion equations (1 6) hav e only one p ossible fixed p oint for every v alue of h . F or h > 0 all h i → j = h + d i − 1, where d i is the degree of node i and magnetization m = 1, for h < 0 all h i → j = h − d i + 1 and m = − 1, and for h = 0 all h i → j = 0 and m = 0. When fixing ma g netization to some v alue − 1 < m ∗ < 1 the third fixed p oint is Belief pr op agation for gr aph p artitioning 10 attractive. B ut provides a wrong v alue of magnetiza tio n m = 0, a s in this case the function M ( h, { h i → j } ) is a s tep function. The fixed po in t also doe s not provide muc h of a useful information ab out the cos t of splitting a tre e on tw o groups of a g iven size. So in a sense our algorithm do es not b ehav e very well on tr ees, that is a kind of un usual situation for b elief propagatio n. The decima tio n algo rithm, how ever, works well and is able to o btain r e asonably go o d partitions even on a tree. This is b ecause once a spin is fixed the infor mation propaga tes and is taken care of co rrectly . But back to Er d˝ os-R´ enyi random graphs, for mean connectivities ab ov e the per colation thres hold but low er than the satisfia bilit y thresho ld 1 < α < α s (given by (5), a nd depicted in Fig. 1) one finds that the algor ithm 2 conv erge s to a c o nfiguration such that o n the g ian t comp onent of the graph all lo cal fields ar e po sitiv e (o r neg a tiv e). Thu s all s pins on the gia n t comp onent will b e chosen to be +1 (-1). In the r est of the graph (that is a ll the small comp onen ts) the lo cal fields as well as the external field h are negative (p ositive) but very close to zero . In s uch a case again, the BP algorithm without decimation is not very efficient in ac tua lly dividing the comp onent s into t wo prop erly s ized g roups. H ow ever, the decimation alg orithm achieves this task quite well (note that if one spin on a small co mponent is fixed, than all the other vertices orient in the sa me direction). 0 0.2 0.4 0.6 0.8 1 1 2 3 4 magnetization, m average degree, α b=0 b>0 glassy phase α s α c RSB threshold satisfiability threshold 0 0.02 0.04 0.06 0.08 0.1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 bisection width, b average degree, α EO , N=2k BP , N=100k 1.35 1.4 1.45 1.5 0 0.002 0.004 Figure 1. Left: The plot shows tw o phase transitions i n partitioning of Er d˝ os- R´ en yi random graphs. The satisfiability threshold α s , eq. (5), ab o ve which the giant component hav e to b e cut to bipartition the graph. And the glass transition, α c , at whic h the belief propagation equations stop to con ve rge and replica symmetry br eaking is needed to descri be correctly the system. Note that α s < α c ( m = 0) = 1 . 472 as anticipat ed i n [22]. Righ t: Bisection width b on Erd˝ os-R´ en yi graph as a function of the mean connec tivity computed by a veraging o ver 2 graphs of size N = 100000 with algorithm 2. The data are compared with the exact av erage bi section w i dth b calculated with the extremal optimization heuristics f or N = 2000, data from [ 22] . As the repli ca symmetri c result provides a low er bound on the energy and the exact ground states on systems of finite sizes are i n this case l arger that the asymptotic v alues, the asymptotic v alue mu st lie betw een the t wo curves. The inset zooms i n to the phase transition r egion. After the satisfiability thresho ld (5), the giant c omponent is bigger than the nu mber of vertices tha t are in the lar ger of the tw o groups, so inevitably one will hav e neig h b ors with opp osite spins in the ground state. There are tw o p ossibilities: (A) BP conv erg es or (B) B P do es not co n verge. If BP do es conv erg e , i.e. b ellow α c , then it converges to a configuration where the giant comp onent is divided into tw o groups (p ositive a nd negative lo c al fields) and all the other comp onents of the gra ph are orie nted in one dir ection (the one that ha s smaller num b er of vertices on the g ian t Belief pr op agation for gr aph p artitioning 11 comp onen t). In o r der to fix the prop er magnetiza tion o n the g ian t c omponent the external field is nonzer o even when the total ma gnetization m ∗ = 0. BP do es not conv erge ab ove the replica sy mmetr y breaking threshold α c depicted in Fig. 1. But even in such ca s es the snapshots of fields are meaningful and the decimation a lgorithm a c hieves g oo d energies e v en w hen the non-co n vergence is ignored, a s illustrated in s ection 4.2. In fact on the Erd˝ os -R ´ enyi ra ndo m graphs there is a firs t or der phase tr ansition at zer o mag netization. At the transition the der iv ative of the energy with r espect to ma g netization has a discontin uity . On b oth sides of the transitio n a meta-s table state exists with spino dal p oints at v alues of ma gnetization corresp onding to the half size of the giant comp onent. This phase transition and lines corre s ponding to the meta-stable state and the spino dal p oin t are illustra ted in the fig ure 2. 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 -1 -0.5 0 0.5 1 b m d=1.6 d=1.44 Figure 2. The figure shows the partitioning cost b as a function of m f or tw o differen t Erd˝ os-R´ enyi graphs with mean connectivities 1 . 44 and 1 . 6 (and of si zes N = 100000). In the simulation, the messages were randomly ini tialized for m = − 1 and then b was calculated with algorithm 2. Magnetization was then slightly increased to m + ∆ m and messages h i → j we re i nitialized with their v alues from s im ulation with previous m . The dashed curv es corresp ond to the case when the system orients the spins on the small comp onen ts in the less fav orable wa y (that is, +1 for m > 0 and vice versa). The dashed curv es end at a spinodal point where the giant component is di vided in half . How to understand this phase transition: Consider large p ositive magnetization, in the lowest cos t s olution the giant comp onent a nd larg e part o f the small co mponents are p ositive and a small part of the comp o nen ts are nega tiv e. As the magnetization is decreased the small comp onen ts are all turning nega tiv e, a nd also parts of the gia n t comp onen t turn negative. The externa l field is negative in that reg ion in order to keep the sma ll comp onents neg ative. Ev en after half of the spins b ecome nega tiv e the system do e s not rea lize that it is les s costly to turn everyb o dy , instead if the magnetization is s lowly decreased further the b elief pr opagation equa tions indicate that a larger frac tio n of the giant co mp onent s hould b e nega tive. As the magnetization is decreased the nega tiv e exter nal field b ecomes clo s er to zero, at the p oin t the exter na l field flips to p ositive v alues the small comp onent turn to p o sitiv e direction and the system realizes this gives muc h low er cost. This po in t c orresp onds to a spino dal po in t. Of co urse this discussio n could b e rep eated by changing the works po s itiv e for negative and vic e versa. The phase transition, meta-sta ble state, and spino dal po in t are illustra ted in figure 2. If the mag ne tiza tion is no t changed gra dua lly , dep ending on the initial conditions Belief pr op agation for gr aph p artitioning 12 the alg orithm do es conv erge to one or the o ther of the br anc hes, more likely to the low e r one. This is a nice prop ert y as if more stable divisions a r e pr e sen t in real net work our a lgorithm might b e able to find them (or a t le a st those o f them with cons iderably large basin of attrac tion). 4.2. Performanc e of the BP de cimation In this section we illustrate acc ur acy of the decimation BP solver on random 3- r egular graphs. Regular graphs are in some sense the ha rdest ca se for gra ph bi-partitioning as they lo ok lo cally alike from every p oint a nd no apparent s tr ucture ca n b e explor ed to decide if tw o no des sho uld be in the same g roup or not. 0.1 0.11 0.12 0.13 0.14 0.15 0.16 100 1000 10000 100000 bipartition cost, b system size, N decimation solver b 1RSB =0.1138 EO, exact average b(N) 0 0.02 0.04 0.06 0.08 0.1 0.12 -1 -0.5 0 0.5 1 bipartition cost, b magnetization, m decimation solver 1RSB BP Figure 3. Left: Decimation results for 3-regular random graphs of differen t sizes, compared to presumably exact a verage ground state energies as computed f rom the extremal optimization heuristics by [14]. Als o shown is the asymptotic cost b = 0 . 1138 calculated by 1RSB method. Note that the decimation algorithm is far better than the b est known algorithmic b ound b = 0 . 1 6 [19]. Right: The plot shows replica symmetric (BP) results, 1RSB results and p erformance of the decimation algori thm for the partition cost b as a function of the magnetization m for 3-regular random graphs. The BP population dynamics algorithm was with N = 10000, 1RSB s olutions were obtained from a simulation with N = 30000. The decimation results were were a veraged ov er 10 different gr aphs, each with N = 2000. If Fig. 3 we show the average bisection co s t achieved b y the decimation so lvers on graphs of different size. W e co mpa re to the asymptotic v a lue of the cost a nd to the av erage v alues obtained fro m extremal optimizatio n heuristic of [14, 33] that are exac t (or at least very clo se to exact), we s ee that our decimation solver achieves energie s very closed to the gro und states. I n pa rticular note that the best prov able a lgorithmic bo und for 3- regular graph bisectio n is b = 0 . 1 6 [19] whic h is far above wha t decima tion achiev es. In the r igh t part of Fig. 3 we compar e the pa rtition cost a s a function of the magnetization m a s obtained fr om (a) the p opulation dynamics solving the BP equations, (b) reso lutio n of the 1RSB equatio ns fro m sec. 2.3 under the assumption that for every edge the distribution of fields P i → j ( h i → j ) is the same — this being called the facto r ized s olution in [29], and (c) decimation so lv er run on gr aphs o f size N = 2000 . Belief pr op agation for gr aph p artitioning 13 4.3. R andom r e gular gr aphs at zer o magnetization In this subsection we want to discuss the bisectio n (zer o magnetiza tion) of ra ndom regular g raphs. This case has b een treated in [2, 23, 24, 25, 2 6, 28] using analo gy with spin gla sses, i.e. Hamiltonian H SG = − X ( ij ) ∈ E J ij S i S j , (28) with r andom J ij = ± 1 has b een so lv ed instead of fixing magnetizatio n to zero via an external field. Indeed, note that in r a ndom regular graph it is mor e than reasona ble to assume that the tw o groups in graph bise c tion hav e exactly the s ame prop erties and hence the first order phas e tra ns ition that w e hav e seen at m = 0 in the Erd˝ os -R ´ e nyi gr aph is exp ected to disapp ear. Conseq uen tly , the s lope of the g round state e ( m ) at m = 0 is exp ected to b e zero , a nd hence als o the v alue of external field to which o ur alg orithm conv erges is zer o h = 0 . W e remind that cavit y fie lds h i → j can b e interpreted as a change in the ground state ener gy of (6) w he n link ( ij ) is remov ed fro m the g raph. If h is an integer then also all h i → j hav e to be integers in the the final solution of the problem. The cavit y equations can then be parameteriz ed by fr action o f neg a tiv e, po sitiv e and z e ro c avity fields h i → j . The only wa y to achiev e zero magnetiza tion is then to set the fraction of negative and p ositive cavit y fields equal. And this lea ds exactly to the sa me equations as M´ e z ard and Parisi obtained in [29] and justifies the a pproach of [2, 23, 24, 25, 26, 28]. Consequences a nd g eneralization of this equiv alence will b e descr ibed in [33]. W e wan t to s tr ess that at non-zer o magnetizatio n the co rresp onding externa l field h do es not take an in teger v alue and hence no str aight forward relatio n to the s pin glas s problem exis ts. Also a s lo ng as the degr e e o f the graph is not constant there mig h t be a ro om for a first orde r phas e transitio n at m = 0 due to a symmetries b et ween the tw o groups in the bise c tion - a s illustr ated in the Erd˝ os-R´ enyi g raphs. If the first order phase transition is present that at m = 0 the exter na l field h 6 = 0 and hence aga in no straightforward analo gy with the spin glass pr oblem exists. Th us the appro ach developed in this pap er is the o nly one know that is able to trea t no n-regular gra phs or non-z e r o v alues of the mag netization. In T able 1 and Fig. 4 we summarize the known rigor ous b ounds for bisection widths in ra ndom r egular gr aphs. W e als o summarize results of b elief propaga tion obtained from our po pulation dyna mics, and the results from 1RSB calculatio n using int eger v alues o f the cavity fields. Bo th the latter are only a pproximation to the full- step r e plica sy mmetry breaking res ult that would pre s umably be exa ct in this case. Finally we compare with p erforma nc e o f o ur decimation BP solver. In particular Fig. 4 illustra tes how accur ate the decimation solver is. Not that the true v alue of the bisection width must lie be tween the decimation a nd 1RSB data p oints. 5. Discussion The main pr actical contribution of this a rticle is the b elief pr opagation algorithm for graph partitioning problem that we anticipate to be useful in the v ario us applications of the pa rtitioning problem. W e s tudied the b ehavior and p erformance of the a lgorithm on r andom gr aphs but we anticipate it will be meaningful a lso for o ther families of graphs, complex netw orks in particula r. Compar ed to other par titioning a lgorithms Belief pr op agation for gr aph p artitioning 14 d b low b up b RS b 1RSB b BPdec 3 0.101 0.1666 0.1125(2) 0.113846 0.1180 (3) 4 0.22 0.3333 0 .2579(2) 0 .263527 0.27 2 (1) 5 0.319 2 0 .5028 0.40 72(3) 0.41 2398 0.422(2) 6 0.480 3 0 .6674 0.57 56(3) 0.58 5414 0.5975 (9 ) 7 0.648 6 0 .8502 0.74 30(4) 0.75 2171 0.766(2) 8 0.822 6 1 .0386 0.92 32(4) 0.93 6595 0.955(2) 9 1.001 2 1 .2317 1.10 22(4) 1.11 453 1.133(1 ) T a ble 1. This table summ arizes the b est known low er b ound (2nd column, [18, 17], the Bollobas’s b ound d/ 4 − p ( d ln 2) / 2 is the best known f or d ≥ 5) and the upper b ound (3rd column, [ 19, 20, 21]) b ounds for random regular graph bisection. In the 4th column we give results for the bisection from the p opulation dynamics for b elief propagation, these n umbers are i den tical to the ones obtained in [ 26] with non-integer cavit y fields. The 5th column gives results of the 1RSB calculation wi th inte ger fields as deve lop ed in [29]. And the final column shown perf ormance of our implementat ion of the BP decimation algorithm for graphs of size N = 2000. -0.85 -0.8 -0.75 -0.7 -0.65 -0.6 -0.55 3 4 5 6 7 8 9 rescaled bisection , (2b - d/2)/d 1/2 degree , d SK upper bound decimation solver 1RSB BP lower bound Figure 4. W e pl ot data fr om T able 1 rescaled as (2 b − d 2 ) / √ d as a f unction of the degree d . According to [2] for large d the true v alues should con ve rge to the ground state energy of the Sherr i ngton-Kirkpatric k mo del, E = − 0 . 763219. BP has the adv a n tage that is provides infor mation ab out pro babilit y with which a certain no de is in a certain groups. It is also able to s ee different lo cally stable divis ions of the gr aph - as illustrated by the first order phase tra ns ition in Erd˝ os- R ´ enyi g raphs at zero magnetizatio n. In real world netw o rks the pa r titioning cost at different v alues of magnetization m may lead to a non- trivial information ab out communities in the net work and information ab out their s ignificance. No te also that our a pproach is straightforwardly genera lizable to k -partitioning the graph into k g roups of a fixed size. Ac knowledgmen t W e thank Stefan Bo ettcher for s haring with us his da ta from the extr emal o ptimization algorithm that we used for compar ison in Figs. 1, 3 . W e thank Cris Mo ore for Belief pr op agation for gr aph p artitioning 15 po in ting to us the mea ning o f the first order phase tr a nsition a t zero ma gnetization and the existence of the s pinodal lines illus trated in Fig. 2. W e also thank Flore n t Krzak ala, Ma rk Newman, Allo n Percus, and F ederico Ricci-T ersenghi for v ar ious very useful dis c ussions ab out this work. P ˇ S acknowledges partial supp ort of New Mexico Consortium via NSF collab orative grant CCF- 0829945 o n “Harnes s ing Statistical Physics for Computing and Communications” References [1] S. Ki rkpatric k, C. D. Gelatt Jr ., and M . P . V ecc hi. Optimization by simulated annealing. Scienc e , 220:671–680, 1983. [2] Y. F u and P . W. Anderson. Appli cation of statistical mecha nics to NP-complete problems i n com binatorial optimization. J. Phys. A , 19:1605–1620, 1986. [3] M. M´ ezard and G. Parisi. Replicas and optimization. J. Physique , 46:L771–L778, 1985. [4] M. M´ ezard, G. Parisi, and R. Zecch ina. Analytic and al gori thmic solution of r andom satisfiability problems. Scie nce , 297:812–81 5, 2002. [5] F. Krzak ala, A. Montanari, F. Ricci- T ersenghi, G. Semerj ian, and L. Zdeborov´ a. Gibbs states and the set of solutions of random constraint satisfaction problems. Pr o c . Natl. A c ad. Sc i. U.S.A , 104:10318, 2007. [6] M. M´ ezard and G. Pa risi . The b ethe lattice spin glass revisited. Eur. Phys. J. B , 20:217, 2001. [7] M.R. Garey and D.S. Johnson. Computers and intr actability: A guide to the the ory of NP- c ompleteness . F reeman, San F rancisco, 1979. [8] G. Karypis , R. Aggarwal, V. Kumar, and S. Shekhar. M ul tilev el hypergraph partitioning: applications in vlsi domain. IEEE T r ansactions on V ery L ar g e Sc ale Inte g r ation (VLSI) Systems , 7(1):6979, 1999. [9] A. Pothen. Graph partitioning algorithms with applications to scientific computing. T ec hnical repor t, Norfolk, V A , USA, 1997. [10] M . Girv an and M. E. J. Newman. Commu nity structure in so cial and biological net works. Pr o c e e dings of the National A c ademy of Scienc es of the Unit e d St ate s of Americ a , 99(12):782 1–7826, 2002. [11] B. Kernighan and S. Lin. An efficient heuristic pro cedure for partitioning graphs. Bel l System T e c hnic al Journal , 49:291 – 307, 1970. [12] J. R. Banav ar, D. Sherr ington, and N. Sourlas. Gr aph bipartitioning and statistical mech anics. J. Phys. A Lett. , 20:L1–L8, 1987. [13] G. R. Sc hreib er and O. C. Martin. Cut si ze statistics of graph bisection heuristics. SIAM Journal on Optimization , 10(1):231– 251, 1999. [14] S. Bo ettc her and A. G. Percus. Extremal optimization for graph partitioning. Phys. R ev. E , 64(2):0261 14, Jul 2001. [15] L. Hagen and A.B. Kahng. New spectral metho ds for ratio cut partitioning and clustering. IEEE T r ansactions on Computer-Aide d Design of Inte gr ate d Cir cuits and Syste ms , 11(9):10741085 , 1992. [16] G. Kar ypis and V. K umar. A fast and high quality multilev el scheme for partitioning i rregular graphs. SIAM Journal on Scienti fic Computing , 20(1):359 – 392, 1999. [17] B. Bollobas. The isop erimetric num ber of random r egular graphs. Eur op e an J. Combin , 9:241– 244, 1988. [18] A. V. Kosto c hk a and L. S. Melni k o v. On b ounds of the bisection width of cubic graphs. In J. Nesetril and M. Fiedler, editors, F ourth Czecho slovakian Symp osium on Combinatorics, Gr aphs and Complexity , pages 151–154. Elsevier Science Publisher, 1992. [19] B. M onien and R. Preis. Upper b ounds on the bisection wi dth of 3- and 4-regular graphs. In MF CS ’01: Pr o c e e dings of the 26th International Symp osium on Mathematic al F oundations of Computer Scienc e , pages 524–536, London, UK, 2001. Springer-V erlag. [20] J. D ´ ıaz, N. Do, M. J. Serna, and N. C. W or mald. Boun ds on the max and mi n bisection of random cubic and random 4-r egular graphs. The or. Comput. Sci. , 307(3):531–5 47, 2003. [21] J. Diaz, M.J. Serna, and N.C. W ormald. Bounds on the bisection width for random d - regular graphs. The or etic al Computer Scienc e , 382(2):120 – 130, 2007. Latin A m erican Theoretical Informatics. [22] A. G. Percus, G. Istrate, B. Goncalv es, R. Z. Sumi, and S. Bo ettc her. The p eculiar phase structure of random graph bisection. Journ al of Mathematic al Physics , 49(12), 2008. arXiv:0808.1549v2 [ cond-mat.stat -mech]. Belief pr op agation for gr aph p artitioning 16 [23] D. Sherrington and K. Y. M. W ong. Graph bipartitioning and the bethe spin glass. Journal of Physics A: Mathematic al and Gener al , 20(12):L785–L791 , 1987. [24] K. Y. M. W ong and D. Sherrington. Graph bipartitioning and spi n glasses on a random net work of fixed finite v alence. Journal of Physics A: Mathematic al and Gener al , 20(12):L793–L799 , 1987. [25] M . M´ ezard and G. Parisi. Mean-field theory of r andomly fr ustrated systems wi th finite connect ivity . Eur ophys. L et t . , 3:1067–107 4, 1987. [26] K. Y. M. W ong, D. Sherrington, P . Mottishaw, R. Dewar, and C De Dominicis. Graph partitioning and dilute spin glasses: the minimum cost solution. J. Phys. A: Math. Gen. , 21:L99–L105, 1988. [27] W. Liao. Replica symmetric solutions of the graph-bipartitioning pr oblem wi th fixed, finite v alence. J. Phys. A: Math. Gen. , 21:427–440, 1988. [28] Y. Y. Goldschmidt and C. De Dominicis. Replica symmetry breaking i n the s pin-glass m o del on lattices with finite connectivit y: Appli cation to graph partitioning. Phys. R ev. B , 41(4):2184– 2197, F eb 1990. [29] M . M´ ezard and G. Parisi. The ca vity metho d at zero temp erature. J. Stat. Phys. , 111:1–34, 2003. [30] C. Di, A . M on tanari, and R. Urbank e. W eigh t distributions of LDPC co de ensembles: Comb inatorics meets statistical physics. In International Symp osium on Information The ory . IEEE, 2004. [31] T. Mora and M. M´ ezard. Geometrical organization of solutions to random linear B o olean equations. J. Stat. Me ch. , 10:P10007, October 2006. [32] F. Krzak ala, F. Ricci-T ersenghi, and L. Zdeborov´ a. Elusive glassy phase i n the random field ising mo del. arXi v:0911.155 1v1 [cond-mat.dis-nn], 2009. [33] S. Bo ettc her and L. Zdeb oro v´ a. A conjecture about max-cut size and bisection width in random regular graphs. in preparation, 2009. [34] J.S. Y edidia, W.T. F reeman, and Y. W eiss. Understanding b elief pr opaga tion and its generalizations. In Exploring Artificial Intel ligenc e in t he Ne w Mil lennium , pages 239–236. Science & T ec hnology Bo oks, 2003. [35] M . M´ ezard and A. M on tanari. Physics, Information, Computation . Oxford Pr ess, Oxf ord, 2009. [36] D. Sherr ington and S. K irkpatric k. Solv able m odel of a spin-glass. Phys. R ev . L ett. , 35:1792– 1796, 1975. [37] G. Parisi. The order parameter for spin-glasses: A function on the interv al 0–1. J. Phys. A , 13:1101–111 2, 1980. [38] J. R. L. de Almeida and D. J. Thouless. Stability of the Sherr ington-Kirkpatric k solution of a spin-glass mo del. J. Phys. A , 11:983–990, 1978. [39] L. Zdeb oro v´ a. Statistical physics of hard optimization problems. A cta Physic a Slovac a , 59(3):169– 303, 2009. [40] S. F ranz, M. Leone, F. R i cci-T ersenghi, and R. Zecc hina. Exact solutions for diluted spin glasses and optimization problems. Phys. R ev. L et t. , 87(12):127 209, Aug 2001.