A tree-decomposed transfer matrix for computing exact Potts model partition functions for arbitrary graphs, with applications to planar graph colourings

A tree-de comp osed transfer matrix for computi ng exact P otts mo del p artition functions for arbitrary graphs, with applications to planar graph colouri ngs Andrea Bedini 1 , 2 ‡ and Jesp er Lykk e Jacobsen 3 1 Dipartimento di Fisica, Univ er sit` a deg li Studi di Milano , I-2 0133, Milano , Italy 2 INFN, Sezio ne di Mila no, I- 20133 , Milano, Ita ly 3 LPTENS, ´ Ecole Nor male Sup ´ erieure, 2 4 r ue Lho mond, 7 5231 Paris, F rance E-mail: andrea. bedin i@mi.infn.it , jesper .jaco bsen@ens.fr Abstract. Combining tree dec o mp o sition and tra nsfer matrix techniques provides a very general a lg orithm for computing exact partition functions of statistica l mo dels defined on arbitrary graphs. The algorithm is particularly efficient in t he c ase of planar graphs. W e illustrate it b y computing the Potts mo del pa r tition functions and chromatic po lynomials (the n umber of prop er v ertex colourings using Q colour s) for large samples of ra ndom planar graphs with up to N = 100 vertices. In the la tter case, our alg orithm yields a sub-expo nential a verage running t ime of ∼ exp(1 . 516 √ N ), a substantial impr ov ement ov er the ex po nential r unning time ∼ exp(0 . 245 N ) pr ovided by the hitherto b est known a lgorithm. W e s tudy the statistics of chromatic ro ots of random pla na r gr aphs in some deta il, compa ring the findings with results f o r finite pieces o f a r egular lattice. P ACS n umbers : 02.70.-c, 05 .50+q, 89 .70.Eg Keywor ds : T ree dec omp osition, transfer matrix, c hromatic p olynomial, random planar graphs 1. In tro duction A t ypical problem in statistical ph ysics is to compute the partition function Z of some mo del with short-ranged in teractions b etw een discrete degrees of freedom defined on the vertices o f some graph. The partit io n function is a w eighted sum ov er the states of these degrees of freedom. If the g r a ph is a regular latt ice (i.e., consists of a large n um b er M of iden tical la y ers) one usually rewrites Z in terms of matrix elemen ts of T M , where T is a t ransfer matrix (or imaginary- t ime ev olution op erator in an equiv alen t path in tegral form ulation) that corresponds to the addition of a single la y er to the ‡ Pr esent address: MASCOS, Department o f Mathematics and Statistics, The Univ ers it y of Melb ourne, VIC, 3010, Australia T r e e-d e c omp ose d tr ansfer matrix 2 lattice. This suggests solving the problem by diagonalising T , so that, in some sense, one has substituted algebraic complexit y fo r com binatoria l complexit y . When the lattice is plana r, the exact diagonalisatio n o f T in the limit of an infinitely large system (or thermo dynamic limit) can in many cases b e achiev ed b y using the p ow erful to ols of in tegrability [1]. Alternativ ely , there exists man y efficie nt means o f diagonalising T n umerically for rather large sys tems. In man y a pplications one is how ev er in terested in mo dels defined on graphs that are not regular, but incorp ora t e some elemen t of randomness. One m ust then distinguish whether t his randomness is of the annealed or t he quenc hed type. F or the annealed disorder, Z is a double sum o v er the g raphs and t he v ertex degrees of freedom. In t he case of planar graphs, this double sum can in many cases b e ev aluated analytically b y random matrix tec hniques [2]. But arg ua bly , the ph ysically most relev an t scenario is that of quenc hed disorder, where one w ould t ypically need to a verage the free energy log Z , not just Z , o v er the disorder realisations. Clearly , a maximal wa y to deal with quenc hed randomness would b e to ev aluate Z indep enden tly for eac h graph in the sample. This is generally not p ossible analytically , ev en for pla nar graphs. The purp ose of this pap er is to exhibit an efficien t algorithm that ev aluat es Z exactly for an arbitrarily giv en graph. This algorithm applies to a ra ther large class of mo dels of statistical phys ics, as outlined ab ov e, and it do es not require the graph to b e planar. Rather than delving in to this generality —we shall br iefly come bac k to the issue in the discussion—w e ha ve here c hosen to fo cus on one significan t sp ecific application, namely the ev aluation of the Potts mo del [3] partitio n function Z G ( Q, v ) (whic h, a fter a c hange of v ariables, equals the T utte p olynomial [4] known in graph theory) for a giv en planar graph G . W e are particularly intere sted in the sp ecial case v = − 1, where Z G ( Q, v ) ≡ χ G ( Q ) equals the n um b er of Q -colourings (c hromat ic p olynomial) of G . W e no w outline the ph ysical motiv ation for studying this problem, b efore ending t he introduction by giving a brief accoun t on t he w orking principle of our algorithm and its p erfo rmance (computational complexit y). The Q -colouring pro blem consists in assigning to eac h of the v ertices of a gra ph G any one of Q differen t colours, in suc h a w ay that a dj a cen t ve rtices carr y different colours. Eac h suc h assignmen t is kno wn as a prop er v ertex colouring. The Q -colouring problem arises within phy sics in studies of frustrated a ntiferromagnetic spin mo dels, and in spin glass theory in particular [5]. The n um b er of po ssible colourings (p ossibly zero) can b e sho wn [6 , 7] to b e a p olynomial in Q , kno wn as the c hromatic p olynomial χ G ( Q ) of the g raph. In particular it mak es sense—and is useful—to generalise the original coun ting problem and to study the pro p erties of χ G ( Q ), considered as a p olynomial o f a f o rmal v ariable Q . The history of the Q -colouring problem is long a nd in teresting, and we refer the reader to [8] for an extens ive list of references . The case where G is a planar gra ph has attracted particular in terest. A w ell-kno wn result is t hen the fo ur- colouring theorem [9] whic h states that χ G (4) > 0 for an y planar G . Other results exploit the ex tension T r e e-d e c omp ose d tr ansfer matrix 3 of χ G ( Q ) to non-in teger v alues of Q . One in teresting ques tio n is whether there exists, in the planar case, some Q c so that χ G ( Q ) > 0 for all Q ∈ [ Q c , ∞ ). This statemen t has b een established as a theorem [1 0 ] for Q c = 5, a nd the so-called Birkhoff-Lewis conjecture [10] t ha t this extends to Q c = 4 is widely b eliev ed to b e true. The case where G is a r e gular pla na r lattice has also b een in tensiv ely studied. In particular, the lo cation and prop erties of c hromatic ro ots, χ G ( Q ) = 0, in the complex Q -plane has b een studied fo r a v ariet y of lattices and b oundary conditions (see [11 , 8, 12, 13, 14 , 15 ] and references therein). Some of the mec hanisms resp onsible for the generation of real c hromatic ro ots in the region Q ∈ [0 , 4], close to or exactly at the so-called Beraha n umbers B k = (2 cos( π /k )) 2 with k ≥ 2 in teger, ha v e eve n b een understo o d ana lytically [16, 17]. One should also men tion that there e xists a fa mily of planar gra phs—planar triangulations actually—with real c hromatic ro ots con v erging t o Q = 4 from b elo w [18], meaning that the v alue o f Q c in the Birkhoff-Lewis conjecture cannot b e lo w ered further. One issue whic h has b een left largely unansw ered b y these studies is the distribution of the lo cation of (real and complex) c hromatic ro ots for ensem bles of randomly chose n planar graphs. This situation app ears to b e the most interes ting from the p o int of view of the physic s of disordered frustrated systems. F ro m the mathematical side, it has b een pro v en [19] that c hromatic root s are dense in the complex plane (exce pt ma yb e in the disc | Q − 1 | < 1) for a sp ecial class of planar graphs (g eneralised theta graphs). Ho w ev er, this might not hav e muc h to do with the ph ysical ques tion of the c hromatic ro ots of a typic a l planar graph. One w ould a lso w an t to kno w whether typical ro ots accum ulate a t the Beraha n um b ers, and whic h graph c haracteristics (connectivit y , av erage co ordination n um b er,. . . ) might b e resp onsible for the lo cation of ro ots. In order to elucidate the amazingly in tricate b ehaviour of c hromatic ro ots in the limit of large graphs, it is clearly desirable to hav e efficien t means of computing χ G ( Q ). Let us recall that in theoretical computer science, an imp ortant outstanding question is whether the class P of decision problems that can b e answer e d in p olynomial time coincides with the class NP of problems fo r whic h a pr o p osed answ er can b e veri fi e d in p olynomial time. NP-complete problems are those to whic h an y problem in NP can b e reduced in p olynomial time. A t presen t, no p olynomial-time algo rithm has b een found for an y of the thousands of kno wn NP-complete problems and it is hence widely b eliev ed that P 6 = NP . Like wise, one can define a counting analo gue o f NP , denoted b y #P , a s t he class of en umeration pro blems in which the ob jects b eing coun ted are the p ossible answ ers a ccepted by an NP machine. If an NP problem is often of the form “Are there an y solutions that satisfy certain constrain ts?”, the corresponding #P problems ask “ho w man y” rather than “are there any”. The class # P-complete of the hardest problems in #P is defined likew ise. Clearly , #P-complete problems a r e at least as hard as NP-complete problems. It is know n [20] that the coun ting of prop er v ertex colourings is a #P-complete pro blem for Q = 3 , and the same is then true a fortiori f or the computation of χ G ( Q ), let a lone Z G ( Q, v ) [21], for generic v alues of Q a nd v . T r e e-d e c omp ose d tr ansfer matrix 4 In practice this means that an y algo rithm computing χ G ( Q ) can b e exp ected to hav e a running time that increase s exp o nen tially with the num b er of v ertices N . Ho w ev er, lo w ering the co efficien t o f the exp onen t can still mak e a h uge difference f o r studying the issues outlined ab o ve . One cen tral goa l of this pap er is to mak e a substantial step in that direction, b y lo w ering the a ve rag e running time ∼ exp(0 . 245 N ) of the previously b est kno wn algo rithm [22] to ∼ exp(1 . 51 6 √ N ), as illustrated in Fig . 1 b elow. The impro v emen t from exp onential to sub-exp onen tial asymptotics is not only imp ortant from a t heoretical p oint of view. T o get a roug h idea of what this means in practice, in ten seconds the algorithm [22] will compute χ G ( Q ) for a t ypical planar graph o f N = 40 v ertices, whereas our a lgorithm can deal with N = 100 in the same time. Algorithmic progress on #P-complete problems related to graph theory and net w ork design has b een made b y sev eral, usually widely separated, comm unities. On one hand, statistical physic ists hav e shown that the relev an t pa r t it ion functions can b e constructed in analogy with the path integral formulation o f quantum mec hanics. T o t his end, the configuration of a par t ially elab orated gra ph are enco ded a s suitable quan tum states, and the constan t-time surface is sw ept o v er the gra ph by means of a time ev olution op erator known as the transfer matr ix. Although rarely stat ed, this approac h is v alid not only for regular latt ices but also for arbitrary graphs. On the other hand, gr a ph theorists hav e used that graphs can b e divided into “w eakly in teracting” subgraphs through a so-called tree decomp osition [23, 24], and solutions o btained for the subgraphs can b e recursiv ely c ombined in to a complete solution. The principle underlying the algorit hmic prog ress to b e describ ed in this pap er is that tree decomp osition and transfer matrix metho ds can b e combined in a v ery natural w a y . The main idea is tha t the tree decomp osition is compatible with a recursiv e generalisation of the t ime ev olution concept. Borro wing ideas from quan tum field t heory , the com bination of partial solutions is obtained by the fusion of suitable state space s. The resulting algo rithm w orks on an y graph, and can readily b e adapted t o many other problems of statistical phys ics, by suitable mo difications of the state spaces and the fusion pro cedure. In particular, w e will apply this tec hnique to the problem of computing the c hromatic p o lynomial on planar graphs, obtaining exact solutions, for graphs with N ≃ 100 v ertices, in only a few seconds. The la y out of the pap er is as follows. I n section 2 w e briefly recall the relatio n b et w een the P otts model partition function Z G ( Q, v ) and the c hromatic p olynomial χ G ( Q ). In section 3 w e presen t our algorithm and discuss its p erformance. The results for c hromatic ro ots of planar graphs are g iv en and discussed in section 4. Finally , we giv e our conclusions in section 5. T r e e-d e c omp ose d tr ansfer matrix 5 2. P ott s model and vertex colou rings W e first recall the relation b et w een the vertex colouring problem and t he P o t ts mo del. Consider a graph G = ( V , E ) with ve rtices V and edges E , and let σ i = 1 , 2 , . . . , Q b e the colour of v ertex i ∈ V . Then Z G = X σ Y ( ij ) ∈ E e K δ ( σ i ,σ j ) (1) is the partition function of the Potts mo del on G . The K r onec k er delta δ ( σ i , σ j ) = 1 if σ i = σ j , and 0 otherwise. Inserting the ob vious iden tity e K δ ( σ i ,σ j ) = 1 + v δ ( σ i , σ j ), with v = e K − 1, and expanding out the pro duct w e obtain the F ortuin- Kasteleyn expansion [25] Z G ( Q, v ) = X A ⊆ E v | A | Q k ( A ) , (2) where k ( A ) is the n um b er of connected compo nen ts in the subgraph G ′ = ( V , A ). Ob viously , in the an tiferromagnetic limit K → −∞ (or equiv alen tly v → − 1) t he only surviving configurations are prop er Q -colourings of the graph G . Indeed, the sp ecial case χ G ( Q ) = Z G ( Q, − 1) is a p olynomial in Q , known as the chr omatic p olynomial , and equals the num b er of v ertex colourings. 3. Algorithm 3.1. T r ansfer matrix The computation of the partition function Z G ( Q, v ) using the expansion (2) has been described in [26, 8] for the case where G is a finite piece of a r e gular lattice (notably , but not exclusiv ely , a planar one). W e first describ e ho w t his traditional tra nsfer matrix metho d extends to the case where G is an arbitrary graph—i.e., not necessarily par t of a regular lattice, nor necessarily planar . In short, the com bined action of linear op erators builds a superp osition of all configurations app earing in the part it io n funciton with their correct statistical (Boltzmann) w eigh t en tering as coefficien ts. T o better illus tra t e this pro cedure, consider the follo wing example g r a ph G 1 2 3 4 5 6 7 8 9 (3) W e first ha v e to define the order { v t } in whic h ve rtices will b e pro cessed. This order is the bas is for the construction of a “time slicing” of the graph. With eac h time step is asso ciated a b ag (a v ertex subset) of a ctiv e vertices . A ve rtex b ecomes activ e as so o n as T r e e-d e c omp ose d tr ansfer matrix 6 one of its neighbours is pro cessed and it sta ys activ e un til it is pro cessed itself. T aking the v ertices in lexicographic order w e obta in the following decomp osition: 1 2 3 2 3 4 3 4 5 6 4 5 6 8 5 6 7 8 9 6 7 8 9 7 8 9 8 9 9 (4) where w e wrote in b old face the v ertex b eing pro cessed at each time step. Eac h bag has its o wn set of basis states consisting of the p artitions § of t he currently activ e ve rtices. F or instance, in the first time step, the basis states are the fiv e pa rtitions of the three-elemen t se t { 1 , 2 , 3 } : | 1 2 3 i , | 1 2 3 i , | 1 2 3 i , | 1 2 3 i , | 1 2 3 i . These partitio ns describ e how the activ e v ertices are in terconnected through A ∩ E t , where E t ⊆ E is the subset of edges having b een pro cessed at time t . A state is a linear sup erp osition of basis states. Pro cessing a vertex consists in pro cessing edges connecting it to unpro cessed v ertices and then deleting it. Since each edge e ∈ E may or ma y not b e presen t in A w e pro cess an edge ( i, j ) b y acting on the state with an op erator of the form 1 + v J ij where 1 is t he iden tit y op erator and J ij a join op er ator . A join op era t o r acts on a basis state b y amalgamating the blo c ks con taining v ertices i and j . J ij | i j i = | i j i , J 2 ij = J ij . (5) V ertex deletion is defined in terms o f a del e tion op er ator D i that remo v es i f r o m the partiton and a pplies a fa ctor Q (resp. 1) if i was (resp. w as not) a singleton. D i | i j · · · i = Q | j · · · i , (6) D i | i j · · · i = | j · · · i . (7) F or example, in (4), pro cessing the first bag means computing the fo llo wing comp osition D 1 ( 1 + v J 12 ) ( 1 + v J 13 ) | 1 2 3 i , whic h giv es ( Q + 2 v ) | 2 3 i + v 2 | 2 3 i , concluding the first time step. When a new activ e v ertex is encoun tered it is inserted, as a singleton, in each partition comp osing the curren t state. After pro cessing the last bag, the complete partition function (2) is obtained as the co efficien t of the empt y par tition resulting from the deletion of the last activ e v ertex. A t eac h step, the time and memory requirem ents are determined b y the b ag size n whic h is the num b er of v ertices sim ultaneously active . If the graph is plana r, the n umber § Recall that a partition of a finite set S is a collectio n of nonempty , pairwise disjoint subsets of S whose union is S . The subsets a re called blo cks . T r e e-d e c omp ose d tr ansfer matrix 7 of par t itions to b e considered is at most the Catalan n umber C n , i.e., the n um b er of non-crossing partitions. The corresp onding generating function is C ( z ) = ∞ X n =0 C n z n = 1 − √ 1 − 4 z 2 z . (8) If the g raph is not planar, the n um b er of par t it ions is at most the Bell n um b er B n , with generating function B ( z ) = ∞ X n =0 B n z n n ! = exp(e z − 1) . (9) W e hav e C n = 4 n n − 3 / 2 π − 1 / 2 [1 + O (1 / n )], whereas B n gro ws sup er-exp onen tially . | 3.2. T r e e de c omp o s ition It turns out that the decomp osition (4) of G is a sp ecial case of a more general construction. By definition, a tr e e de c omp osition [23, 24] of a graph G = ( V , E ) is a collection of bags, orga nised as a tr ee (a connected graph with no cycles), and satisfying the follo wing requiremen ts: (i) F or eac h i ∈ V , there exis ts a bag con taining i ; (ii) F or eac h ( ij ) ∈ E , there exists a bag containing b o th i and j ; (iii) F or an y i ∈ V , the set of bags con taining i is connecte d in the tree. The previous decomp osition (4) is just a sp ecial case o f a tree decomp osition (a p ath de c omp osition ). As an example of the general construction, applied to (3), consider 1 2 3 2 3 4 3 4 5 3 5 6 5 6 7 4 5 8 5 8 9 (10) where the arro ws form the unique path that connects eac h ba g to the cen tral one (the r o ot of the tree). The adv an tage of w orking with tree instead o f path decompo sitions relies o n the fact that in the former case a decomp osition with smaller bags can b e obtained (the latter b eing just a sp ecial case). Therefore, the num b er of states one has to k eep tr a c k of is exp onen tia lly smaller, and the gain is significant. The transfer matrix approa ch can b e adapted naturally to t his new general setting: Prop erties (i)–(ii) gua r a n tee that each edge and v ertex are processed within a definite bag. Prop erty (iii) implie s that each ve rtex has a definite life time in the recursion, its insertion and deletion b eing separated by t he pro cessing of all edges inciden t on it.   Note that the restriction to non-crossing partitions holds true even if the ordering of the vertices, such as in (4), do es not resp ect the planar embedding (3) of the whole graph (whic h will in general be unknown). Because we choose to s tore par titions independently of their cro ssing pr op erty , it has not b een r e lev an t to b e able to explicitly keep track of the planar embedding. All that matters is that the n umber of par titions do es not exceed C n , and this is true since such a non-cro ssing repres ent a tio n could be obtained “by hand” b y insp ecting the plana r embedding at each stag e of the tra ns fer pro c e ss. T r e e-d e c omp ose d tr ansfer matrix 8 In this new ve rsion the algorithm starts from the ro ot of the tree, whic h can b e c hosen a rbitrarily , and runs throug h the tree recursiv ely . Alternativ ely , this can b e though t of as starting from the lea v es and building up the t ree inductiv ely . Consider first the case o f a paren t bag P with only o ne da ugh ter bag D . Going up from D t o P in the recursiv e step inv olv es deleting vertices D \ P , inserting v ertices P \ D and finally pro cessing edges in P . A tree dec omp osition do es not sp ecify when an edge e = ( ij ) m ust b e processed. A simp le recip e w ould b e to pro cess e as so on as one encoun ters a bag containing b oth i and j . Ho w ev er w e note that this f r eedom of c hoice can b e exploited to optimise the algorithm (see b elo w). 3.3. F usion W e now discuss the case when a paren t bag P has sev eral daugh ters D ℓ with ℓ = 1 , 2 , . . . , d . In this case , v ertex deletions and insertions a r e follo w ed by a sp ecial fusion pro cedure. Supp ose first d = 2, and let P 1 b e a partitio n of D 1 ∩ P with w eight w 1 , a nd P 2 a partition of D 2 ∩ P with w eight w 2 . The fused state is then P 1 ⊗ P 2 = P 1 ∨ P 2 , and it o ccurs with w eigh t w 1 w 2 . Here ∨ denotes the j o in o p eration in the pa r t ition lattice. ¶ F or definiteness, let us explain ho w this can b e computed in practice. First c ho ose some set E 1 so that P 1 = Y e ∈ E 1 J e ! S 1 , (11) where S 1 is the all-singleton partition o f D 1 ∩ P . Note that it is not necessary to choo se E 1 as a subset of E , although this can alw a ys b e done. Since J e J e ′ = J e ′ J e the order in the ab o v e pro duct is irrelev ant. Similarly c ho ose E 2 for P 2 . The fused state c an then b e constructed explicitly as P 1 ⊗ P 2 = Y e ∈ E 1 ∪ E 2 J e ! S 12 , where S 12 is the all-singleton partition of ( D 1 ∪ D 2 ) ∩ P . F or d > 2 daugh ters, the complete fusion can b e accomplished b y fusing D 1 with D 2 , then fusing the result with D 3 , and so on. ¶ W e recall that a partia l order of the partitions follows from defining a partition P a to b e a r efinement of P b , and we write P a  P b , provided that each blo ck in P a is a subset o f so me blo ck in P b . With this partial order , the set of par titio ns fo r m a lattic e , in the sense that any tw o partitions P 1 and P 2 po ssess a la rgest low er b ound called me et a nd denoted P 1 ∧ P 2 and a sma llest upp er b ound called join and denoted P 1 ∨ P 2 . The meet P 1 ∧ P 2 has a s blo cks all nonempty intersections of a blo ck fr om P 1 with a blo ck from P 2 . The blo cks of the join P 1 ∨ P 2 are the s ma llest subsets which are exactly a union of blo cks from b oth P 1 and P 2 . In the partition lattice, the b ottom (o r finest) element is the unique partition in which each blo ck has size 1 (the all-single ton partition), and the top (or coars est) e le ment is the unique pa rtition in w hich there is a s ingle blo ck (the a ll-connected partition). T r e e-d e c omp ose d tr ansfer matrix 9 Let us illustrate the fusion pro cedure for G with the tree decomp osition (10). Aft er pro cessing the t w o left-most bags and deleting vertex 2, the propagating state is ω 1 | 3 4 i + ω 2 | 3 4 i = ( ω 1 + ω 2 J 34 ) | 3 4 i , (12) where ω 1 = Q 2 + 3 v ( Q + v ) and ω 2 = Q 2 v + 3 Qv 2 + 4 v 3 + v 4 . By symmetry , the same result is obtained for the tw o top righ t bags (replacing 4 b y 5) and for the t wo b ot t o m righ t bags (replacing 3 by 5). + The fused state arriving in the cen tral bag is then ( ω 1 + ω 2 J 34 )( ω 1 + ω 2 J 35 )( ω 1 + ω 2 J 45 ) | 3 4 5 i = ω 3 1 | 3 4 5 i + (3 ω 1 ω 2 2 + ω 3 2 ) | 3 4 5 i + ω 2 1 ω 2 ( | 3 4 5 i + | 3 4 5 i + | 3 4 5 i ) , (13) from whic h the result Z G ( Q, v ) follows up on deleting v ertices 3,4,5. 3.4. Pruning Problem sp ecific features can be exploited to reduce further the n um b er of basis states to b e cons idered. As an example of this, note that in the c olour ing case ( v = − 1) the op erator O ij = 1 + v J ij asso ciated with an edge ( ij ) is a pro jector, O 2 ij = O ij , and it annihilates the s ubspace of partitions where i and j are connecte d. It follow s that one can discard basis states in whic h t w o v ertices are connected, as so on as one discov ers that an edge b et w een them is g oing to b e pro cessed la t er within the same bag or within the pa ren t ba g. Especially b efo re fusions this simple tric k reduce s substan tially the n um b er of basis states a nd th us leads to a big sp eed up. 3.5. Perf o rmanc e F or a p lana r gra ph, the state of a ba g of siz e n is spanned b y C n basis states. (F or a non-planar graph, replace C n b y B n .) The memory needed by the algorithm is therefore prop ortional to C n max , where n max is the size of the largest bag. The time needed to pro cess one edge in a bag o f size n is prop ortional to C n . Ho w ev er, most of the time is sp en t fusing states. F or a paren t P with d daugh ters D ℓ , the num b er of basis state pairs to b e fused is d X ℓ =1 C |D ℓ − 1 ∩ P | C | D ℓ ∩ P | , (14) where we ha v e set D k = ∪ k ℓ =1 D ℓ . Eac h of these elemen tary fusions (i.e., the joins P 1 ∨ P 2 ) can b e done in time linear ∗ in the n um b er of participating v ertices. Note that w e can c ho ose the order of successiv e fusions so as to minimise the quan tit y (14 ). + Needless to say , the fusion pro cedure will work also in cases where the g raph p osse s ses no such symmetries. In our example, we hav e chosen G as a r ather symmetric g raph in o rder to keep the actua l computations as simple as po ssible. ∗ W e repr esent partitions as lists of n umbers linking ea ch pa rticipating vertex to the blo ck it belo ngs to. The s et E 1 in (11) can b e obtained in a single scan of this list w her e one keeps tr ack of the last T r e e-d e c omp ose d tr ansfer matrix 10 It is therefore essen tial for t he algorithm t ha t one knows ho w to obtain a g o o d tree decomp osition. The minim um of n max − 1 ov er all tree decomp ositions is kno wn as the tree width k , but obtaining this is another NP-hard pro blem. Ho w ev er, the simple algorithm GreedyFillI n [27] giv es an upp er b ound k 0 on k and a tree decomposition of width k 0 in time line ar in the num ber of v ertices N . F or uniformly g enerated pla nar graphs we find that for N = 40—a v alue enabling comparison with algo rithms that determine k exactly—that k 0 = k nearly alw a ys ( k 0 = k + 1 w ith probabilit y ≃ 10 − 3 ). When k 0 is small enough that a ll the pa r titions can fit in to the computer’s memory , the algorithm has prov en to b e ve ry fa st with execution time in the or der of seconds. T o summarise, supp ose that we wish to compute Z G ( Q, v ) or χ G ( Q ) for a planar graph of N v ertices and M edges, and that w e are giv en a tree decompo sition of B bags, with n max b eing the size of the largest bag. The num ber of (binary) fusions to b e p erformed is F = P i ( d i − 1 ), where d i is the num ber of daughters of bag i . F or simplicit y we a ssume a computational mo del where t he o p erations on the w eights can b e done in unit time. ♯ F or the v ersion of the a lgorithm without pruning, we hav e therefore sho wn that the worst-case running time is asymptotically prop ortio nal to ( N + M + B ) C n max n max + F C 2 n max n max . (15) Note that this is linear in N a nd M for a fixed n max (cf. [28]). W e c ho ose to test our algorithm aga inst the one presen ted b y Haggard et al. in [22]. W e first generated a uniform sample of 100 planar graphs for eac h size N = | V | b et w een 20 and 100 using F usy’s algorithm [31]. W e then ran four differen t algo rithms o v er this sample: the algo rithm of [22], our first path-based transfer matrix algorit hm, the new tree-based v ersion algorithm and a tree-based v ersion using the ab o v e pruning optimisation. P ath decompo sitions were obta ined with a v arian t of the G reedyFillIn algorithm in whic h the resulting tr ee decomp o sition w as forced not to branch. Av erage running times are presen ted in Fig. 1, in the form of effectiv e expo nential fits. While these fits clearly confirm the efficiency of the tree decomp osed approac h, they are actually misleading a nd hide an imp ortant fact. Namely , the tree width of a planar g raph is b o unded b y α √ N , as follo ws from the famous planar separator theorem [29]. The currently b est upper b ound on the constan t is α < 3 . 1 82 [30]. These facts—com bined with (15) and (8), and the heuristic observ ation of the near-ideality of the tree decomp ositions obtained by t he GreedyFill In heuristics—immediately imply a n upper bound on the w orst- case running time on the tree-based algorithms of 16 3 . 182 √ N = exp( 8 . 822 √ N ). It follo ws that the mean running time of these algorithms m ust b e o f the f o rm exp( β √ N ), with the constants β satisfying β tree+pruning ≤ β tree . The fits for β , using t he same data as in Fig. 1, a re sho wn in Fig. 2. seen vertex for ea ch blo ck. In this repres e ntation, applying J e is cons tant time, apart from a p o s sible linear-time pro cedure to bring the re s ulting list in a cano nical form. This last pro cedure can b e safely po stp oned after all the join op er a tions, making the whole elemen tar y fusion linear in the nu mber of participating vertices. [Note that in [28] the join op eration is claimed to b e quadr atic in n .] ♯ The w eig h ts are num be rs for an ev aluation of Z G ( Q, v ) or of χ G ( Q ), p olynomia ls with in tege r co efficients for χ G ( Q ), and p oly nomials in tw o v aria bles with integer co efficients for Z G ( Q, v ). T r e e-d e c omp ose d tr ansfer matrix 11 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 1 0 0 N 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 1 0 1 1 0 2 average running time (seconds) avg (path) ~ exp(0.406 N) avg (tuttepoly) ~ exp(0.245 N) avg (tree) ~ exp(0.141 N) avg (tree + pruning) ~ exp(0.103 N) Figure 1. Average running time in seco nds on a random planar gr aph of size N vertices. E ach p oint is averaged over 1 00 g r aphs. 3.6. T e chni c al asp e cts F or completeness w e give a few implemen tational details. P artitions relev ant to the propagating state are k ept in a ha sh table for fast (amortised constant time) a ccess. The corres p onding w eigh ts ar e p olynomials in Q whose co efficien ts rapidly exceed the mac hine integer size M = 2 32 when running on big gra phs. T o solv e this issue without requiring a dditional memory w e used mo dular arithmetics: t he algorit hm w as run man y times with co efficien ts computed mo dulo primes p < M , and the original co efficien ts w ere reconstructed b y the Chinese remainder theorem. 4. Results As a simple application of our algorithm w e o btained the distribution of the c hromatic ro ots { Q | χ G ( Q ) = 0 } for ensem bles of planar graphs of up to N = 1 0 0 vertice s. This problem is in teres ting b oth in statistical mec ha nics and in graph the ory [4, 16, 17, 8]. As known f rom Lee-Y ang theory , the ro ots signal phase tra nsitions. Before turning to the actual results, it is useful to give some results fo r r e gular planar graphs of a few hundred v ertices. T o this end, w e ha ve used a particular T r e e-d e c omp ose d tr ansfer matrix 12 4 5 6 7 8 9 1 0  N 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 1 0 1 1 0 2 average running time (seconds) ¯ t t u t te p o l y  e 0 . 2 4 5 N ¯ t t r e e  e 1 . 8 4 2  N ¯ t t r e e + p r u n i n g  e 1 . 5 1 6  N Figure 2. Average running time in seco nds on a random planar gr aph of size N vertices, now fitted to the form exp( β √ N ). The data p oints a re identical to those of Fig. 1 . adaptation of the traditiona l (path-decomp osed) tr a nsfer matrix algorithm that ta k es sp ecial adv an tage of the regular lattice structure [32] to compute χ G ( Q ) for L × L sections of the square and triangular lat tice (the latter b eing considered as a square latt ice with added diagonals). The b oundary conditions are free–free, in the terminology of [8]. The c hromatic p olynomials for L = 10 , 12 ha v e b een c hec k ed against [12], while those for L = 14 , 16 , 18 a re new. The resulting lo cations of the chromatic ro ots in the complex Q -plane are shown in Fig. 3 for the square lattice, and in Fig. 4 for the triang ular lattice. As in earlier w ork [8, 12, 13] w e no tice that most o f the ro ots fall on conne cted curv es which, v ery roughly speaking, t end to form a e gg- shap ed figure with a pair of prongs on the righ t side. Without going in to details (whic h are num erous, see [33, 11, 16, 17, 8, 12, 13, 34] for more exact statemen ts) in the limit L → ∞ the egg- shap ed pa rt is exp ected to enclose a segmen t Q ∈ [0 , Q 0 ] of the real axis, with Q 0 = 3 for the square lattice , and Q 0 ≈ 3 . 81967 for the t r ia ngular lattice. In the in terior of the egg one observ es additional isolated real zeros right at, or extremely close to, the so-called Beraha num b ers B k = (2 cos( π /k )) 2 ( k = 2 , 3 , 4 , . . . ) (16) T r e e - d e c o m p o s e d t r a n s f e r m a t r i x 1 3 0 0.5 1 1.5 2 2.5 3 Re(Q) -2 -1 0 1 2 Im(Q) 10 x 10 12 x 12 14 x 14 16 x 16 18 x 18 Figure 3. Chromatic ro ots for L × L sections of the squa re lattice with free b oundar y conditions, a nd L = 10 , 12 , 14 , 16 , 18. Beraha num b ers B k with k = 2 , 3 , 4 , 5 are shown as small vertical line s e g ments. The reason wh y this is so is link ed to par t icularit ies o f the represe ntation theory of the quan tum algebra that underlies the tw o-dimensional P otts mo del [16, 17]. The existence of a finite-size equiv alen t of Q 0 is clearly visible fro m Figs. 3 – 4. One could pro p ose defining Q 0 ( L ) a s the la rgest r eal ro ot, but one should b e careful since in some cases the egg-shap ed curv e do es no t p ossess a ro ot righ t on the real axis. Barring these (and other) difficulties, one could speculate that if the result for an arbitrary planar g raph w as qualitatively similar to that of t hese regular la ttices, the probability distribution of complex ro ots w ould lo o k lik e some broadened-out egg. In additio n, the distribution of real ro o ts would b e a sup erp osition of sharp p eaks centere d at the Beraha n um b ers, and a broad bac kground distribution correspo nding to the Q 0 ( L ). This is indeed more-or-less what w e observ e. F or a sample of 2500 uniformly draw n [31] planar graphs o f N = 100 v ertices w e show the distributions of complex c hromatic ro ots in Fig. 5, and the distribution of the subset of real ro ots in Fig. 6. R egarding the complex root s, w e note that although c hromatic roo ts of planar graphs ha v e b een T r e e - d e c o m p o s e d t r a n s f e r m a t r i x 1 4 0 1 2 3 4 Re(Q) -2 -1 0 1 2 Im(Q) 10 x 10 12 x 12 14 x 14 16 x 16 18 x 18 Figure 4. Chromatic ro ots for L × L sectio ns o f the triangular lattice with free bo undary conditions, a nd L = 10 , 12 , 14 , 1 6 , 18 . Beraha num be rs B k with k = 2 , 3 , 4 , 5 , 6 , 7 , 8 are shown as small vertical line segments. sho wn to b e dense in the complex plane [1 9] (except ma yb e in the disc | Q − 1 | < 1), t ypical ro ots are clearly quite close to the origin. As to the re al c hromatic ro ots, t he absense of ro ots on the negative real axis, a nd the in terv als (0 , 1) and (1 , 3 2 / 27] fo llo w from a theorem [35 ] (se e also [36]). The ro ots found here resp ect this theorem a s w ell as t he Birkhoff-Lewis conjecture [10]. Apart from that, Fig. 6 show s as exp ected a sup erp osition of a broad bac kground distribution and sharp p eaks cen tered at Q = B k with k = 2 , 3 , 4 , 5 , 6 (we ha ve B 5 ≃ 2 . 6180 3). It s eems likely t hat for N → ∞ this bac kground will ex tend to the in terv al (32 / 27 , 4) and p eaks will o ccur a t all Beraha n um b ers. W e also exp ect that the maxim um of the bac kground d istribution ma y b e shifted furt her to the righ t by requiring the graphs to b e t wo- o r three-connected, and so presumably the p eaks a t the first few B k should stand out more clearly . W e plan to in ve stigate this and further issues in more detail elsewhere [37]. T r e e-d e c omp ose d tr ansfer matrix 15 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 2 . 5 3 . 0 3 . 5 x  3  2  1 0 1 2 3 y 1e-03 1e-02 1e-01 1 10 p ( x , y ) Figure 5. Distribution for complex chromatic r o ots Q = x + iy obtained from a sample of 250 0 plana r graphs of size 1 0 0. The nor ma lisation is such that p ( x, y ) is the exp ected n umber of ro o ts in [ x − 0 . 05 , x + 0 . 05] × [ y − 0 . 05 , y + 0 . 05] counted with their m ultiplicity . 5. Conclusion In this pap er w e ha ve sho wn ho w to o btain an efficien t alg orithm fo r solving instances of graph-related #P-complete problems. The case o f the P otts partition f unction Z G ( Q, v ), and in particular its sp ecialisation the c hromatic p olynomial χ G ( Q ) = Z G ( Q, − 1), w ere studied in detail, and some results on the distribution of real and complex chromatic ro ots for unifo rmly dra wn random planar graphs were given. Our algorithmic appro a c h has some ov erlap with a n algorithm describ ed—but apparen tly not implemen ted—b y Noble [28]. This author also com bines partial ev aluat io ns of Z G ( Q, v )—treated by him in the T utte p olynomial formulation [4 ]—in the basis of par titions with t he notion of tree decomp osition. There are how ev er a n um b er of imp orta nt differences. F ir st, our use of tra nsfer matrices—and sp ecifically of its factorisatio n within eac h bag, i.e., o ur c hoice to pro cess a single edge at a time within eac h bag — leads to a b etter time complexit y . In par t icular, treating a daugh ter bag with no sisters has w orst-case running time C 2 n max n 2 max + C n max 2 n max ( n max − 1) / 2 in Noble’s T r e e-d e c omp ose d tr ansfer matrix 16 1 . 0 1 . 5 2 . 0 2 . 5 3 . 0 3 . 5 4 . 0 x 1 0 - 4 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 1 0 1 1 0 2 p ( x ) Figure 6. Pro ba bility distribution for the subset of chromatic r o ots in Fig. 5 that lie on the r eal axis. Normalised to unit total area . approac h ve rsus C n max n max in o urs. Second, w e treat binar y fusions more efficien tly , exploiting the fact that a given pair P 1 , P 2 of input partitions giv es rise to a unique output partition P 1 ∨ P 2 : this reduces the time complexit y for binary fusions from w orst- case time C 3 n max n 2 max in [28] to C 2 n max n max in our approac h. Third, our algo rithm app ears easier to describ e and implemen t—as w e hav e indeed done—a nd it av oids having to deal with the Q → 0 limit ( x = 1 in the notations of [28]) as a particular case. It is th us clear that what matters the most is to obtain a g o o d tree decomp osition, and in particular the running time is ess entially determined by n max . In Fig. 7 we sho w the distribution of n max obtained fro m applying the algorit hm GreedyFillIn to random planar graphs as a function of N . It is plainly visible from Fig . 7 that b o t h the mean and w orst- case n max exhibit a slo w er than linear g r owth with N . Assuming that GreedyFillIn yields alw a ys an n max whic h is of the order of the t r ue tree width, w e w ould hav e n max ≤ cst √ N , and the data of Fig. 7 app ear to b e compatible with suc h a scaling. W e defer the further inv estigation of the asymptotic b ehaviour of GreedyFillIn to future work [37]. W e should also mention that the use of tree decomp osition in the contex t of de cision problems has previously b een adv o cated b y a num b er of authors (see [38] a nd references T r e e-d e c omp ose d tr ansfer matrix 17 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 N 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0 n max 1e-05 1e-04 1e-03 p ( N , n max ) Figure 7. Proba bility distr ibution p ( N , n max ) of obtaining a tree decomposition of maximum bag size n max (i.e., width n max − 1) when applying the Gre edyFi llIn algorithm to random planar graphs, a s a function of the gra ph size N . The data is computed o ver a sample of 10 0 planar gra phs for each size b etw een 20 and 500, binned in blo cks of 1 0 o n the x a xis and nor malised to unit total ar e a. therein). In these w orks, the solution is found b y doing dy namic programming on the tree-decomp osed graph. Let us conclude b y commen ting on the generality o f o ur algorithm. G oing through the list of kno wn NP-complete problems related to graphs and net w ork design, one easily realises that most of them can b e promoted to coun ting problems (whic h cannot b e easier), and a coun ting p olynomial ana lo gous to χ G ( Q ) can b e defined. Adapting our algorithm to t ho se cases usually requires just a problem-sp ecific definition of the states, of the single-edge transfer pro cess, and of the fusion pro cedure, whereas t he bulk of the metho d can b e tak en o v er without c hanges. It is worth underlining that our ability to solv e efficien tly an enu meration problem can b e exploited to obtain a sp ecific instance that solv es the corresp onding decision problem with only an additional linear time factor. In particular, w e ha v e adapted (but not y et implemen ted) our algorithm to coun ting v ersions o f the following problems: Hamiltonian w alks and cycles, lo ng est path, ve rtex co v er, dominating set, feedbac k v ertex set, and minim um maximal indep endent set. W e plan to rep ort o n these problems elsewhere [37]. T r e e-d e c omp ose d tr ansfer matrix 18 Ac knowled gmen ts This w ork w as supp orted b y the Europ ean Comm unity Net work ENRAGE (grant MR TN-CT-200 4-00561 6) and b y the Agence Natio nale de la Rec herc he ( g ran t ANR- 06-BLAN-012 4-03). The autho r s thank A.D. Sok al for discussions and for commen ting on a n earlier v ersion of the man uscript. References [1] R.J. Baxter , Exactly solve d mo dels in statistic al me chanics (Academic P ress, Lo ndon, 1 982). [2] M.L. Mehta, R andom matric es , 3rd edition (Els evier, Amsterdam, 2004). [3] R.B. Potts, Pr o c. Ca mb. Phil. So c. 48 , 106 (1952 ). [4] W.T. T utte, Cana d. J. Math. 6 , 8 0 –91 (195 4). [5] R. Mulet, A. Pagnani, M. W eight and R. Zecchina, Phys. Rev. Lett. 89 , 2687 01 (2002 ). [6] G.D. Birk hoff, Ann. Math. 14 , 42 (191 2). [7] H. Whitney , Bull. Amer. Ma th. So c. 38 , 57 2 (193 2). [8] J. Salas and A.D. Sok al, J. Sta t. Phys. 104 , 609– 699 (2001). [9] K. App el and W. Haken, Bull. Amer. Math. So c. 82 , 711 (1976). [10] G.D. Bir khoff a nd D.C. Lewis, T rans. Amer. Math. So c. 60 , 355 (1946 ). [11] R.J . Ba xter, J. P hys. A 2 0 , 5 241 (1 987). [12] J .L. J a cobsen a nd J. Sa la s, J . Stat. P hys. 104 , 701 (20 01). [13] J .L. J a cobsen, J . Salas and A.D. Sok al, J. Stat. Phys. 112 , 921 (2 003). [14] J .L. J a cobsen a nd J. Sa la s, J . Stat. P hys. 122 , 705 (20 06). [15] J .L. J a cobsen a nd J. Sa la s, Nucl. P hys. B 7 83 , 2 38 (20 07). [16] H. Sa leur, C o mm. Math. P h ys . 132 , 65 7 (199 0). [17] H. Sa leur, Nucl. P hys. B 3 60 , 2 19 (19 91). [18] G.F. Royle, Anna ls Combin. 12 , 1 95–21 0 (20 08). [19] A.D. Sok al, Combin. Pr obab. Comput. 13 , 221 –261 (2004 ). [20] D.J.A. W elsh, The c omputational c omplexity of some classic al pr oblems fr om st atistic al physics , in G.R. Grimmett and D.J.A. W elsh (eds.), Disor der in physic al systems , Clarendon Press, Oxford, 1990, pp. 30 7–321 . [21] F. J aeger, D. V ertiga n a nd D. W e ls h, Math. Pro c. Camb. Phil. So c. 1 08 , 3 5–53 (1990 ). [22] G. Hagga rd, D.J. Pearce and G. Royle, Computing T utt e p olynomials . T echnical r ep ort NI0 9024- CSM, Isaac Newto n Ins titute for Ma thema tical Scienc e s , 2009 . [23] N. Ro ber tson and P .D. Seymour , J. Comb. Theory B 36 , 49-64 (1984 ). [24] H.L. B o dlaender, Theo r . Co mp. Science 2 09 , 1 –45 (1 998). [25] C.M. F ortuin a nd P .W. K asteleyn, Physica 57 , 536– 5 64 (1972). [26] H.W.J. B l¨ ote a nd M.P . Nigh tinga le, P hysica A 1 12 , 4 05 (19 82). [27] H.L. B o dlaender and A.M.C.A. Ko ster, Infor m. Co mput. 208 , 2 59–2 75 (2 010). [28] S.D. Noble, C o mbin. P robab. Co mput. 7 , 3 0 7–32 1 (1998). [29] N. Alon, P . Sey mour and R. Tho ma s, J . Am. Math. So c. 3 , 80 1–808 (199 0). [30] F.V. F o min and D.M. Thilikos, J . Graph Theory 5 1 , 5 3 –81 (2006). [31] E . F us y , Rando m Struct. Alg. 35 , 4 64–5 22 (2009 ). [32] J .L. J acobsen, Bulk, surfac e and c orner fr e e ener gy series for the chr omatic p olynomia l on the squar e and triangular lattic es , arXiv:100 5.3609 . [33] R.J . Ba xter, J. P hys. A 1 9 , 2 821 (1 986). [34] R.J . Ba xter, Pro c. Roy . So c. London A 383 , 43 (198 2). [35] B . Ja ckson, Combin . Pro bab. Comput. 2 , 32 5–33 6 (199 3). [36] C. Thomas sen, Co mbin . Pro bab. Comput. 6 , 49 7–506 (199 7). [37] A. Be dini and J.L . Ja cobsen, in prepa ration. T r e e-d e c omp ose d tr ansfer matrix 19 [38] F.V. F o min and D.M. Thilikos, Le c ture Notes in Comput. Sci. 299 6 , 5 6–67 (20 04).