Topological invariant tensor renormalization group method for spin glasses
aa r X i v : . [ c ond - m a t . d i s - nn ] O c t Topological Invariant Tensor Renormalization Group Method for Edwards-AndersonSpin Glasses Model
Chuang Wang,
1, 2, ∗ Shao-Meng Qin, and Hai-Jun Zhou State Key Laboratory of Statistical Physics, Institute of Theoretical Physics,Chinese Academy of Sciences, Zhong-Guan-Cun East Road 55, Beijing 100190, China University of Chinese Academy of Sciences, Beijing (Dated: July 9, 2018)Tensor renormalization group method (TRG) is a real space renormalization group approach. Ithas been successfully applied to both classical and quantum systems. In this paper, we study a dis-ordered and frustrated system, the two-dimensional Edward-Anderson model, by a new topologicalinvariant TRG scheme. We propose an approach to calculate the local magnetizations and nearestpair correlations simultaneously. The Nishimori multi-critical point predicted by the topologicalinvariant TRG agrees well with the recent Monte-Carlo results. The TRG schemes outperform themean field methods on the calculation of the partition function. We notice that it maybe obtaina negative partition function at sufficiently low temperatures. However, the negative contributioncan be neglected if the systems is large enough. This topological invariant TRG can also be used tostudy three-dimensional spin glass systems.
PACS numbers: 75.10.Nr, 05.10.Cc, 64.60.DeKeywords: Edward-Anderson model, spin glasses, tensor renormalization group
I. INTRODUCTION
Exploring the Edwards-Anderson (EA) model is sig-nificant but extremely difficult. The nature of spinglasses on three-dimension is still a heat debate betweenthe mean field picture and droplet picture . over thepast 30 years. For the two-dimensional (2D) model, be-sides its interests in statistical physics, it has wild appli-cations on image processing , computer vision , whichis usually referred to as the Markov random field inthe computer scientists community. In this paper, wepropose a coarse-graining method for EA model on 2Dsquare lattice and calculate local physical quantities si-multaneously by the tensor renormalization group (TRG)method.TRG is a real space renormalization group approachinitially introduced by Levin and Nave for classicalferromagnetic Ising spin systems on 2D regular lattices.This method is an extension of the density matrix renor-malization group method for one-dimensional quantumsystems . The basic idea is to perform a coarse-grainingprocess on a tensor network. Matrix low rank approxi-mation is used to cut the degree of freedom of tensorindices up to a maximum value D through the singularvalue decomposition.Shortly after the introduction of the initial TRGmethod, an improvement was made by Xiang and co-authors , who proposed a backward iteration to cal-culate the environment tensor and improved the resultsby considering the effect of the environment. The TRGmethod has excellent performance on the classical ferro-magnetic Ising model, the Potts model and the dilutedferromagnetic model , etc. It also becomes a crucialtool to handle 2D quantum systems . Very recentlya further improvement, namely the topological invariantTRG method, was proposed in the papers to extend the TRG to three-dimensional (3D) ferromagnetic Isingcases.Unlike the ferromagnetic Ising model, the EA spinglasses model is heterogeneous, disordered and frus-trated. It is intrinsic hard. The problems of finding aground state of the 2D EA model with external fieldand the general 3D EA model are proved to belong toNP-hard class , which is commonly believed that noalgorithm can solve them within polynomial time. Inthe previous study, the mean field approximation and Monte Carlo Sampling , transfer matrix method ,numerical exact algorithm for 2D without the externalfield are used to calculate local properties for indi-vidual finite size instances. These methods are combinedwith finite size scaling to investigate the thermodynam-ics limit properties. The duality relationship andreal space renormalization methods are also employedto study the phase diagram and universality. TRG canbe exploited in both of two roles. It can be served as anapproximate calculator of physics quantities for a singleinstance, and it may also be used as a new renormaliza-tion method to directly investigate critical phenomenon.We, here, focus on the former role. To our best knowl-edges, there is no work on applying TRG on spin glassesuntil now.In this paper, we proposed two main approaches.Firstly, we show a new topological invariant coarse-graining scheme based on the work . It avoids two prob-lems when the method is directly applied on EA model:cutting extra freedom of indices and inversing singularmatrices. In the ferromagnetic Ising model, these twoproblems don’t exist. Secondly, we propose an approachto compute local physical quantities simultaneously. Forexample, all single-spin magnetizations can be calculatedby a single sweep of coarse-graining procedure and back-ward procedure. These two approaches are also usefulfor other heterogeneous systems. In the numerical cal-culation, TRG may get a negative value of the partitionfunction at very low temperature, which is the majordifference between the spin glasses model and other het-erogeneous systems . We show that the contributionof the negative part is comparable to the error flunctu-ation for a large enough system, and therefore it can beneglected. In the high temperature region, TRG out-performs the mean field method, belief propagation andgeneralized belief propagation , while the mean fieldmethods are failed in the lower temperature because ofthe convergence problems. The Nishimori multi-criticalpoint is calculated by our TRG scheme. The re-sults agree well with the recent Monte-Carlo results .We emphasis that the original TRG method can alsobe directly applied to any heterogeneous systems, includ-ing EA model, similar to the works on diluted ferromag-netic model . The advantage of topological invariantscheme is that it can be extended to 3D cases .The paper is arranged as follows. In the remainder ofSection 1, we introduce the EA model and show how toconvert it to a tensor network. In section 2, we demon-strate our topological invariant TRG procedure. In sec-tion 3, we show how to calculate local physical quantitiesby backward iteration. In section 4, we list some nu-merical results to test the validation of this method. Insection 5, we discuss the further improvement and appli-cations. II. THE MODELSA. The Edward-Anderson Model
We consider the classical 2D EA model on a pe-riodic square lattice with discrete coupling constants.The system consists of N spins { σ i } , M coupling con-stants { J ij } and N local external fields { h i } . Each spin σ i takes value from { +1 , − } . The overall spins state σ = ( σ , σ , . . . , σ N ) is referred to as a configuration.The energy function is defined as H ( σ ) = − X ( ij ) ∈ E J ij σ i σ j − X i ∈ V h i σ i , (1)where E and V denote the edge set and vertex set of thesystem, respectively.For a single instance of the EA model, the couplingconstants and local external fields are fixed accordingto predefined distributions. In this paper, the value of J ij is randomly chosen from the binomial distribution P ( J ij ) = pδ ( J ij ,
1) + (1 − p ) δ ( J ij , − δ ( x, y ) isthe Kronecker delta symbol, which is 1 if x = y , oth-erwise is 0 . The model parameter 0 . ≤ p ≤ p = 0 .
5) to thepure ferromagnetic system ( p = 1). The configuration σ is supposed to follow the Gibbs-Boltzmann distribution: p ( σ ) = 1 Z exp (cid:2) − βH ( σ ) (cid:3) , s j s k Φ ( ij ) Φ ( im ) Φ ( ik ) Φ ( il ) s m s l s i (a) ˜ V ( im ) s m s l ˜ U ( im ) s ik s k s im ˜ V ( ik ) ˜ U ( ik ) ˜ V ( ij ) ˜ U ( ij ) s j s ik ˜ V ( il ) s i s il ˜ U ( il ) (b) s ik T k s im s ij T m T l s il T i T j (c) FIG. 1. Construction of a tensor network: (a) The neighbor-hood of a vertex i . (b) Each matrix φ ( ij ) is split into twomatrices by the singular value decomposition, so that eachvertex i is now surrounded by four matrices which share thecommon index s i . (c) Summing over the index s i , the neigh-bor four matrices contract to be a tensor T i . where Z = P σ exp [ − βH ( σ )] is the partition function.It is useful to rewrite the distribution as a production ofa set of non-negative weight factors p ( σ ) = 1 Z Y ( ij ) ∈ E ψ ij ( σ i , σ j ) Y i ∈ V ψ i ( σ i ) , (2)where the weight factors have the form ψ ij ( σ i , σ j ) =exp[ βJ ij σ i σ j ], ψ i ( σ i ) = exp[ βh i σ i ]. If all the externalrandom fields are zero, the partition function and thepair-spin correlations can be calculated exactly in poly-nomial time . However, for general external fields { h i } , the problem is proved to be in the NP-hard class . B. Tensor Networks
Any two-body interaction system can be transformedinto a tensor network, in which the partition function ofthe system is equal to the trace of all the tensors. Thetransformation is not unique. Here we show a symmetricmethod. The transformation of the EA model on 2Dsquare lattice at a site i is illustrated in Fig. 1. Firstlyeach Ising spin σ i is mapped to a Boolean variable s i =(1 − σ i ) / ∈ { , } , so that each weight factor ψ ij ( σ i , σ j )can be expressed as a matrix Φ ( ij ) , where the elementin s i -th row and s j -th column is Φ ( ij ) s i s j = ψ ij (1 − s i , − s j ). Note that the C-programming-language conventionis used, in which the index starts from 0. Meanwhile,each external weight factor ψ i ( σ i ) of field h i is mappedto a vector Φ ( i ) , of which the s i -th element is Φ ( i ) s i = ψ i (1 − s i ). Next step, we perform the singular valuedecomposition on the matrix Φ ( ij ) , such thatΦ ( ij ) s i s j = X s ij U ( ij ) s i s ij d s ij V ( ij ) s j s ij , (3)where the matrices U ( ij ) , V ( ij ) are real orthogonal ma-trices and the vector d = ( d , d ) stores singular valuesin descending order. Each element in the vector d isnon-negative. The new variable s ij ∈ { , } is the in-dex of the singular vector d . Let ˜ U ( ij ) s i s ij = U ( ij ) s i s ij d s ij ,˜ V ( ij ) s j s ij = V ( ij ) s j s ij d s ij . Then we haveΦ ( ij ) s i s j = X s ij ˜ U ( ij ) s i s ij ˜ V ( ij ) s j s ij . Now, each variable i is surrounded by four matrices˜ U ( ij ) s i s ij , ˜ U ( ik ) s i s ik , ˜ V ( il ) s i s il , ˜ V ( i,m ) s i s im , where j , k , l , m are labels ofthe neighbor spins of the spin i . Finally, we sum over s i and get a tensor T is ij s ik s il s im : T is ij s ik s il s im = X s i ˜ U ( ij ) s i s ij ˜ U ( ik ) s i s ik ˜ V ( il ) s i s il ˜ V ( i,m ) s i s im Φ ( i ) s i . (4)The partition function of the original system is equal tothe result obtained by tracing over all the indices of thetensors defined on lattice sites: Z = X { s } Y i T is ij s ik s il s im . (5)We refer to the network of tensors constructed by theabove procedure as a tensor network. On the originallattice, each vertex i is associated with a tensor T i , andeach edge ( ij ) is associated with a tensor index s ij . Ingraphical language, the tensor network is similar to a fac-tor graph model with weight factors defined on verticesand state variables defined on edges, but a key differenceis that the elements of tensors are not necessarily non-negative. In the following discussions, we rewrite thetensor indices as i , i , i , i , i.e., T ii i i i for notationalsimplicity. III. TENSOR COARSE-GRAININGPROCEDURE
There are several ways to implement the tensorcoarse-graining procedure . Generally, each coarse-graining iteration consists of two steps. First is contract-ing two neighbor tensors into a new tensor with biggerindices freedom degree. It is an exact procedure. If thereis no computation limitation, the exact partition function j i j i ′ i j ′ j ′ i T ′ Tk j i ′ i ′ (a) ˆ R ˆ kR i ′ j j ′ ˆ i ′ ˆ i i (b) ˆ i ˜ T ′ ˜ Tk ′ j j ′ ˆ i ′ i ′ i (c) FIG. 2. (Color online) Demonstration of TRG: The top figureis the microscope of the circled region in the bottom figure.The two vertical tensors T and T ′ in (a) are contracted intoone tensor R in (b), and the associated two indices i and j of(a) are combined into one index ˆ i . If the degree of freedomsof the index ˆ i is larger than the cut-off parameter D , we usethe singular value decomposition to truncate this index andobtain the approximate tensors ˜ T and ˜ T ′ in (c). Bold linesindicates the freedom of associated indices are greater thenthe others, when the freedom exceeds the cut-off parameterD. could be got by the iteration of these contractions. Sec-ond is cutting the indices freedom degree approximately,so that the computation is tractable.We introduce our method for the tensor network de-fined on a 2D square lattice with the periodic boundarycondition expressed as Eq. (5). At the first step, eachtwo vertical neighbor-tensor pair T, T ′ are contracted asshowed in Fig. 2a. We sum over the common index k ,and the pair T, T ′ is unified into one tensor R : R ( i ,j ) ,j , ( i ,j ) ,i = X k T i ,k,i ,i T ′ j ,j ,j ,k . (6)The new tensor R has 6 indices i , i , i , j , j , j . Wecombine two indices in the same direction i , j as a unionindex ˆ i and i , j as another union index ˆ i , so that thenumber of indices of R is still 4, i.e., ˆ i , j , ˆ i , i . Afterthe contraction, the topological structure of the squarelattice is preserved, and the y-direction length shrinks tohalf, while the degrees of indices freedoms associated withthe edges along the x-direction increases to the square ofthe previous one.At the second step, the union indices ˆ i and ˆ i will betruncated alternatively along x-direction if their freedomdegrees are greater then a given cut-off parameter D .Specifically, let us consider the two horizontal neighbortensors R ˆ k,j , ˆ i ,i and R ′ ˆ i ′ ,j ′ , ˆ k,i ′ in Fig. 2b, which sharea same index ˆ k . We think R and R ′ as a sub-system inthe tensor network with the internal variable ˆ k and theboundary variables { ˆ k, j , ˆ i , i } and { ˆ i ′ , j ′ , ˆ k, i ′ } . Theboundary variables interact with other tensors, which canbe considered as the environment of the sub-system. Weare going to approximate the sub-system by another onewith a fewer freedom degree of internal variable such thatthe interaction with environment is as similar as possi-ble. Mathematically, it is done by the lower rank matrixapproximation. We rearrange the indices order of thetensor R as j , ˆ i , i , ˆ k , and group the first three indicesas an unique index i = ( j , ˆ i , i ). Then tensor R be-comes a matrix R i, ˆ k . In the same way, we get the matrix R ′ ˆ k,i ′ from the tensor R ′ , where i ′ = (ˆ i ′ , j ′ , i ′ ). We sumover the common index ˆ k to get a new matrix A : A i,i ′ = X ˆ k R i, ˆ k R ′ ˆ k,i ′ . (7)The sub-system is now expressed by the matrix A . Toexactly represent the boundary interaction, the minimumfreedom degree of the internal variable is the rank of A .A lower rank approximation is made by the singular valuedecomposition. The matrix A is decomposed in the re-duced form by A i,i ′ = rank( A ) − X k ′ =0 U i,k ′ d k ′ V i ′ ,k ′ . (8)The reduced singular value decomposition discards thezero elements of the singular vector d , which has no con-tribution to the sub-system. In the numerical compu-tation, singular values less than the criterion d i < ǫ =10 − are considered to be zero. If the rank of A is greaterthan the cut-off parameter D , we only keep the largest D singular values. Let a ′ = min { rank( A ) , D } . The ap-proximation of A is expressed as A i,i ′ ≈ ˜ A i,i ′ = a ′ − X k ′ =0 ˜ T i,k ′ ˜ T ′ i ′ ,k ′ , (9)where ˜ T i,k ′ = U i ′ ,k ′ d k ′ and ˜ T ′ i ′ ,k ′ = V i ′ ,k ′ d k . The matri-ces ˜ T and ˜ T ′ are non-singular, and therefore their inversematrices always exists. This property will be used in thenext section. Next, we expand the grouped indices i and i ′ and rearrange the order of indices to recover the tensor˜ T ˆ k,j , ˆ i ,i and ˜ T ′ ˆ i ′ ,j ′ , ˆ k,i ′ .By this way, the tensors R and R ′ in Fig. 2b are re-placed by ˜ T and ˜ T ′ in Fig. 2c, and the common index ˆ k isreplaced by k ′ whose degree of freedom is no greater than D . The above procedure of cutting off variable ˆ k by thesingular value decomposition guarantees that the matrix˜ A is a best approximation of A among all matrices withthe rank no greater than D , if the measure of the error isthe Frobenius norm k A − ˜ A k F . Note that A is a D × D matrix. The complexity of directly decomposition of A is O ( D ) . Considering that the rank of A is at most D , we could reduce the complexity into O ( D ). Detailsare illustrated in the appendix A. We now rotate the present tensor network 90 ◦ inFig. 2c, and then it has the same local structure of thetensor network as the one at the first step in Fig. 2a,while the length along x direction is reduced by half. Werepeat the step 1 and step 2 once more. The size of thetensor network shrinks half both in x and y directions.This is the complete step of a coarse-graining proce-dure. We repeat it, until the tensor network is reducedsmall enough to be tractable by brute-force summationto get the partition function. In this paper, the final sizeis 2 × A ( i ) in Eq. (9) inthe present layer tensor network to be a fixed value S m ,and we save the logarithm of the scale factor for the i ’thmatrix at the l ’th step as φ ( i ) l = ln( d ( i )0 ) − ln( S m ) , (10)where d ( i )0 is the maximum singular value of matrix ˜ A ( i ) .The total free energy density is f ( β ) = − N β X l X i φ ( i ) l + log Z r ! , (11)where Z r is the remaining scaled partition function cal-culated by contracting the final 2 × O ( D ), whilethe original method and the higher order TRG are O ( D ) and O ( D ), respectively. Practically, the pre-cision in calculating the free energy is better than theoriginal one for the same cut-off parameter D . Ourtensor coarse-graining method is based on the higher or-der TRG , where the exact contraction step is same, butthe approximate truncation is different. The higher or-der TRG method truncates all the indices associated withx-direction edges by the higher order singular value de-composition. However, we found that such a truncationscheme can’t report a sufficiently precise free energy den-sity value for the EA spin glass model. We only truncatethe x-direction indices alternatively, while the remaininghalf of the x-direction indices will be contracted in thenext step, so they are not necessary to be truncated. IV. MARGINAL PROBABILITY ANDBACKWARD PROCEDURE
The EA model has no translational symmetry, andtherefore the local magnetization depends on vertex po-sition. The marginal probability distribution of a vertex i is given by: P i ( s i ) = 1 Z X all indices T i ( s i ) Y j ∈ V \{ i } T j , (12)where s i is related to the spin σ i by σ i = 1 − s i andthe term, all indices, under the summation is referred toas all the indices of every tensor in the tensor network { T i | i ∈ V } , and T i ( s i ) is a tensor at vertex i when itsspin σ i is fixed to 1 − s i : T ii i i i ( s i ) = ˜ U s i i ˜ U ′ s i i ˜ V s i i ˜ V ′ s i i Φ ( i ) s i . (13)As showed in Eq. (12), P i ( s i ) can be computed by ordi-nary TRG method for any i in the tensor network witha special tensor T i ( s i ). However, it is impractical to cal-culate the marginal probabilities for all the vertices inthis way. In this work we use the backward iterationmethod to compute the marginal spin probability dis-tribution functions for all the vertices simultaneously.We define the environment tensor, or just called theenvironment, of a local tensor T i as M ii i i i = X all indicesexcept i ,i ,i ,i Y j ∈ V \ i T j , (14)where the summation is taken over all the indices of thetensor network expect the indices of the tensor T i . Anenvironment M i has the same indices as its correspon-dent tensor T i . The partition function can be re-writtenas Z = X i i i i T ii i i i M ii i i i . (15)And the marginal probability distribution is expressed as P i ( s i ) = 1 Z X i i i i T ii i i i ( s i ) M ii i i i . (16)Similarly, the nearest-neighbor pair-wise marginal distri-bution P ij ( s i , s j ) can be also expressed as a summationbetween a pair of neighbor tensors and the correspondenvironment: P ij ( s i , s j ) = 1 Z X i j j i j i k T ii ,k,i ,i ( s i ) T jj ,j ,j ,k ( s j ) × ˆ M ij ( i ,j ) ,j , ( i ,j ) ,i , (17)where ˆ M ij ( i ,j ) ,j , ( i ,j ) ,i is the environment of the tensor R ( i ,j ) ,j , ( i ,j ) ,i = P k T ii ,k,i ,i T jj ,j ,j ,k .We calculate environments of a tensor network at amore detailed level based on knowing the environmentsat a coarse-grained level, which we called the backwarditeration. We start from the final coarse-grained 2 × M i of a tensor T i at this level can be calculated directly by tracing theother three tensors. Given the environments M ˜ T , M ˜ T ′ of the tensors ˜ T ,˜ T ′ in Fig. 2c, we now show how to calculate the envi-ronments M T , M T ′ of the tensors T , T ′ at the detailedlevel in Fig. 2a. The definition of indices is the same asdescribed in the previous section and shown in Fig. 2.We start from the relation equation of the tensor ˜ T andits environment M ˜ T in Eq. (15) Z = X k ′ ,j , ˆ i ,i ˜ T k ′ ,j , ˆ i ,i M k ′ ,j , ˆ i ,i (18)= X k ′ ,k ′′ ,i ˜ T k ′ ,i δ k ′ k ′′ M k ′′ ,i , (19)where in the second line we group the indices ( j , ˆ i , i )as i , and insert a Kronecker delta function δ k ′ k ′′ .The tensor ˜ T ′ ˆ i ′ ,j ′ ,k ′ ,i ′ can be viewed as a matrix ˜ T ′ k ′ ,i ′ ifwe exchange the order of indices to k ′ , ˆ i ′ , j ′ , i ′ and groupthe indices ˆ i ′ , j ′ , i ′ as i ′ . As mentioned in previous sec-tion, the matrix ˜ T ′ k ′ ,i ′ is always non-singular, we replacethe Kronecker delta function in Eq. (19) by δ k ′ k ′′ = X i ′ ˜ T ′ k ′ ,i ′ ( ˜ T ′− ) i ′ k ′′ , (20)where ˜ T ′− is the inverse of the matrix ˜ T ′ . The partitionfunction is then expressed by Z = X i,j A i,i ′ M ˜ Ai,i ′ , (21)where ˜ A i,i ′ = P k ′ ˜ T k ′ ,i ˜ T ′ k ′ ,i ′ is defined in Eq. (9), and M ˜ A is the environment of ˜ A : M ˜ Ai,i ′ = X k ′′ (cid:16) ˜ T ′− (cid:17) i ′ ,k ′′ M k ′′ ,i (22)Since ˜ A is the lower rank approximation of A , where A i,i ′ = P ˆ k R ˆ k,i R ′ ˆ k,i ′ defined in Eq. (7), the environment M ˜ A is approximately the environment of M A M A ≈ M ˜ A (23)The physical explanation of the above approximation isthat, for a sub-system expressed by ˜ A and its environ-ment, if we replace this sub-system with another sub-system A , which interacts with the environment in a verysimilar way with more internal variable states, the envi-ronment will not change too much. From the relationshipof A and its environment, we get Z = X i,i ′ , ˆ k R ˆ k,i R ′ ˆ k,i ′ M Ai,i ′ . The environments of R and R ′ are obtained as M R ˆ k,i = X i ′ R ′ ˆ k,i ′ M Ai,i ′ , (24) M R ′ ˆ k,i ′ = X i R ˆ k,i M Ai,i ′ . (25)We expand the grouped indices i , i ′ of matrices M R and M R ′ and exchange the indices to get the environ-ments M R ˆ k,j , ˆ i ,i and M R ′ ˆ i ′ ,j ′ , ˆ k,i ′ of tensors R and R ′ . Thisis the backward iteration of the cut-off step.The backward iteration of the contraction step is morestraightforward. We unpack the indices ˆ k and ˆ i of R in the view before contraction, hence R ˆ k,j , ˆ i ,i → R ( i ,j ) ,j , ( i ,j ) ,i = P k T i ,k,i ,i T ′ j ,j ,j ,k as shown inEq. (6), in which the first index ˆ i is the index ˆ k here.From the relation of the tensor R and its environment Z = X i j j i j i k T i ,k,i ,i T ′ j ,j ,j ,k M R ( i ,j ) ,j , ( i ,j ) ,i , we can get the environments of T and T ′ as M Ti ,k,i ,i = X j j j T ′ j ,j ,j ,k M R ( i ,j ) ,j , ( i ,j ) ,i , (26) M T ′ j ,j ,j ,k = X i i i T i ,k,i ,i M R ( i ,j ) ,j , ( i ,j ) ,i . (27)After the above two steps, the environment matrix M ˆ T is calculated by knowing the environment matrix M ofhigher coarse-grained level tensor network. We repeatthis process until the environment tensors of the originaltensor network are obtained. Then the marginal proba-bility distributions can be calculated from Eq. (16). Inpractice, we reduce the computational complexity by uti-lizing the fact that the matrix A is at most rank D .The backward iteration is initially introduced to designa better way to do tensor coarse-graining by minimizingthe change of the whole system with the environments on the ferromagnetic Ising model. This improvement canalso be applied to EA model in the same way. We hereexploit the backward iteration to calculate local physicalquantities simultaneously. V. NUMERICAL RESULTS
We compared the partition function calculated by ourtopological invariant TRG with those obtained by theoriginal TRG and mean field approach, belief propa-gation and generalized belief propagation (GBP) ,on the pure spin glass model without external fields, i.e. p = 0 . h i = 0. The exact partition function iscalculated by the algorithm . The paramagnetic solu-tions of BP and GBP is included, which is the meanfield method under the Bethe-Peierls approximation and Kikuchi approximation respectively. We measurethe average error of the logarithm partition function as ǫ φ = 1 N h| log( Z exact ) − log( Z ) |i (28)over 64 instances with L=64 in Fig. 3 in the region β = 1 /T ∈ [0 , . D = 8, our -12-8-4 0 0 0.2 0.4 0.6 0.8 1 l og ( ε Φ ) β BPGBPTRG(D=8)pTRG(D=8)TRG(D=16)
FIG. 3. (Color online) Comparison of the error of N log Z ,calculated by our topological invariant TRG (pTRG) with D = 8, the original TRG method (TRG) with D = 8 , . The resultsare obtained by averaging, over 64 instances on a periodicsquare lattice with side length L = 64. topological invariant TRG are more accurate than theoriginal TRG. If one use a larger cut-off parameter D ,the results will be better, while the computation timewill increase dramatically.At low temperatures T , i.e. high inverse temperature β = 1 /T , we found that the TRG procedures may re-sult in a negative partition function. This phenomenonhappens both in the original TRG and our topologicalinvariant TRG. We tested 128 instances with the inversetemperature β ranging from 0 to 4 .
0. The probability ofnegative partition function is are showed in Fig. 4a. Abrief explanation is that the elements in the tensors donot constrained to be non-negative and the lower rankmatrix approximation makes the final result to fluctuatearound the exact partition function. At low tempera-tures, the error is so large that the scaled partition func-tion Z r of the finial 2 × D to reduce the probability of negative results. Ifone only cares about the asymptotic result for a largesystem, one could simply neglect the negative part, sincefor infinite system the log partition function is dominatedby the scaling factors Φ ( i ) l in each forward iteration steprather than the remaining contribution Z r . To clarifythis point, we define the ratio of remaining log partitionfunction and the leading part of the scaling factors as r , r = * log | Z r | P l P i Φ ( i ) l + , (29)where h·i means averaging over disorders. Numerically,we averaged 128 instances. As showed in Fig. 4, thecontribution of remaining free entropy decreases as thesystem size increase almost linearly in the log-log scale.For a large system, it will be even lower then the error, sothat we can safely discard this term. This phenomenon P r ob . o f neg . Z r β TRG(D=8 )pTRG(D=8)TRG(D=16) (a) -5-4-3-2-1 3 4 5 6 7 8 9 10 l og (r) log (L) (b) FIG. 4. (Color online) (a) The probability to obtain a nega-tive Z r and (b) The ratio r of log | Z r | to the leading part at β = 1 . -10-8-6-4-2 0 -1 -0.5 0 0.5 1 l og ( ε c ) < σ i σ j > exact D=8D=16
FIG. 5. (Color online) Comparing the nearest neighborcorrelations of a single instance L=64 with exact results at β = 1. The cut-off parameter is set to D = 8 and 16. also indicates that we can investigate the EA model inthe thermodynamic limit, similar to the work on ferro-magnetic Ising model . Because of the heterogeneity,the properties of the system are captured by infinite iter-ations of population of tensors rather than single tensoriteration. We leave the analysis of infinite systems in ourfuture work.We plot all the nearest-pair-spin correlations of a typi- U p L=16L=32L=64L=128 FIG. 6. (Color online) Estimation MNP point by finite sizescaling. Lines are got by fitting Eq. (33). cal single instance compared by the numerical exact val-ues which are calculated by numerical differential of thefree energy at β = 1 .
0. The error is defined by ǫ c = ( |h σ i σ j i pTRG − h σ i σ j i exact | ) (30)Larger cut-off parameter D will lead to better results, asshown in Fig. 5. We do not show the local magnetizationssince they are always zero because of the spin symmetryin the absence of no external fields.The p - T phase diagram of 2D EA model has been ex-tensively investigated in the papers and the ref-erences therein. There is no spin-glass phase occurred atfinite temperature , while it undergo a para-ferro mag-netic phase transition at low temperature T and large p .The system is in the paramagnetic phase when 0 . ≤ p
1. A special line p NL ( T ) = (tanh(1 /T ) + 1) / , on which some physical quantitiescan be calculated exactly. The multi-critical Nishimoripoint (MNP) is the crossing point of the Nishimori lineand the critical line p c ( T ). We compute the MNP bylocating the crossing point.We use the topological invariant TRG as a tool to cal-culate magnetizations, and compute susceptibility χ bynumerical differential χ = d P i h σ i i dh , (31)where h is the external field and h·i means averaging overthe Boltzmann distribution, which can be quickly calcu-lated by the marginal distribution Eq. (16) after the back-ward iteration. The MNP point is estimated by finitesize scaling stated in the work . We measure the RGinvariant quantity U , along Nishimori line near MNP,where U = [ χ ][ χ ] − , (32)where the square brackets are referred to as the averageover the disorder, i.e. the couplings { J ij } . We use 2 × TABLE I. Location of the multi-critical Nishimori pointMethods p ∗ BP instances for each point. Then, the MNP point is gotby fitting U = U ∗ + a ( p − p ∗ ) L y + a ( p − p ∗ ) L y , (33)where U ∗ , a n , p ∗ , y are fitting parameters. We fit thedata with the lattice size 16 ≤ L ≤
128 as showed inFig. 6. and estimiate the MNP point at p ∗ = 0 . ± . y = 0 . ± .
022 and other pa-rameters U = 0 . ± . a = − . ± . a = 6 . ± .
6. The chi-square test reports a small ratio ofchi-square to the degree of freedom χ /DOF = 7 . / L ≥
32, and L ≤
64. All test areconsistent with each other, excpet the data of L = 8,which has strong finite size effect so that we discard itin all fit. The susceptibilities χ are checked by usingdifferent differential steps δh ranging from 10 − to 10 − .For most of instances, they are insensitive to δh , andwe set δh = 10 − . While, a tiny fraction (about 10 − )depends on δh , and for these cases a larger δh is used.The location of MNP is not depend on the choices of δh .Small portions of instances are also verified by averag-ing the two-point correlations. The comparison of theestimation MNP is showed in Table I. The results agreewell with the recent Monte Carlo method with finite sizescaling , and the recent duality analysis inspired by hier-archical lattice . We leave the discussion of re-entrancephenomena and strong disorder universality as the futurework. We emphasis that the role of TRG here is a newtool to calculate physical quantities. Compared othermethods, the mean field estimation by BP and GBP on2D EA model , it improve quite lot. VI. DISCUSSIONS AND CONCLUSION
In this paper, we applied the TRG on the 2D EAmodel, and proposed a novel topological invariant ten-sor coarse-graining procedure, as well as an approach tocalculate local physical quantities simultaneously. Twoproblems hidden in the translation symmetric cases aresolved. We avoid to over-cut the freedom of indices inthe coarse-graining procedure and avoid to inverse a sin-gular matrix in backward iterations. The backward iter-ation process was used to compute single spin marginal probability distributions and nearest neighbor spin paircorrelations.We found that the TRG scheme is able to compute thefree energy and local correlations accurately if the tem-perature is not very low. At low temperatures the TRGscheme might lead to a negative value of the partitionfunction. We show that, for large systems, the main con-tribution of the partition function is the scaling factorsduring the coarse-graining iteration, and the negative re-maining scaled partition function of the final 2 × , thoughoriginally TRG is considered only be applied to gappedphase . The present TRG scheme can’t be applied tothe case at zero temperature, because the SVD only pre-serves local optimal coarse-graining mode, and they areorthogonal in the further coarse-graining iteration andfinally get zero partition function. It’s an open ques-tion that whether TRG can be used at zero temperatureproblem. A further improvement can be made by con-sidering the effect of environments, which is illustratedin the paper on the ferromagnetic Ising model. Inprinciple, one can investigate the fixed point of TRG.However, it may not get a good precision as we did inour paper, because the advantage of TRG is its excel-lence performance on compute physical quantities ratherthan analyzing the fix point of the renormalization .The topic on the nature of spin glass phase on 3Dlattice is still rather active . The main method inmost of the current studies is the Monte-Carlo Sampling.The topological invariant coarse-graining iteration canbe done in 3D cases by contracting tensors along one di-rection, and cut off the indices associated to the edgesalong the other two directions. Local physical quanti-ties, for example the Edward-Anderson parameter, canbe directly got as showed in this paper. The sample-to-sample overlap distribution or other non-local quantitieswould be estimated by TRG guided sampling, in the waythat we fix the spins one by one according to its marginalprobability. So, it presents an alternative way to investi-gate 3D spin glass models.Another application is on investigating combinatorialoptimization problems on finite-dimensional lattices orloopy random graphs. TRG can be immediately appliedon image segmentation and denoising . They share thesame mathematical structure as 2D spin glasses model.For random graph model, mean field method providedexcellent solutions on mean-field like systems, such aslocal-tree like structured graph and fully connectedgraph . While, for the system rich in local loops, themean field approximate may not quite accurate, for ex-ample small world networks, and many real networks.The extension of TRG on general graph provides a newinsight and maybe another physics-contributed solutionto such problems. Similar to the belief propagation,the decimation and reinforcement approaches can becombined with TRG to get optimization solutions. i j i i j j R i,k k R ′ k,j (a) U (2) k,l ′ V (1) k,l U (1) i,l j i j i j kl l ′ V (2) j,l ′ i U (1) i,l ˜ A l,l ′ j j i i j i V (2) j,l ′ (c) U (1) i,l i j i j i j k ′ V ˜ Al ′ ,k ′ V (2) j,l ′ U ˜ Al,k ′ (d) FIG. 7. Demonstration of simplifying cut-off step
ACKNOWLEDGMENTS
We thank Qiao-Ni Chen, Zhi-Yuan Xie, Jin-Fang Fan,Jack Raymond, Victor Martin-Mayor for helpful dis-cussions, and thank Jack Raymond for comments onan earlier version of the manuscript. H.-J. Zhou weresupported by the National Basic Research Program ofChina (No. 2013CB932804), the Knowledge InnovationProgram of Chinese Academy of Sciences (No. KJCX2-EW-J02), and the National Science Foundation of China(grant Nos. 11121403, 11225526). The numerical simu-lations were performed in the HPC computer cluster ofITP-CAS and we thank Dr. Hongbo Jin for technicalhelps.
Appendix A: Simplifying Singular ValueDecomposition of Matrix A We started from the definition of the matrix A inEq. (7), where R and R ′ are tensors with four indices.We exchange and combine the indices so that R ˆ k,j , ˆ i ,i , R ′ ˆ i ′ ,j ′ , ˆ k,i ′ change to be matrices ˆ R ( j , ˆ i ,i );ˆ k , ˆ R ′ (ˆ i ′ ,j ′ ,i ′ );ˆ k .For simplicity we write i = ( j , ˆ i , i ), and i ′ = (ˆ i ′ , j ′ , i ′ ).Instead of multiplying R and R ′ , here we firstly decom-pose them by the singular value decomposition R i,k = X l U i,l d l V l,k , (A1) R ′ k,i ′ = X l ′ U ′ k,l ′ d l V ′ l ′ ,i ′ . (A2)Let ˜ A l,l ′ = X k d l V l,k U ′ k,l ′ d l . (A3)We decompose ˜ A by the singular value decomposition˜ A l,l ′ = X k ′ U Al,k ′ d Ak ′ V Al ′ ,k ′ . (A4)Then tensors ˜ T i,k ′ , ˜ T ′ j,k ′ in Eq (9), could be calculatedby ˜ T i,k ′ = X l U i,l U Al,k ′ d A l , (A5)˜ T ′ i ′ ,k ′ = X l ′ d A l V Al ′ ,k ′ V ′ l ′ ,i ′ . (A6)The numerical SVD routines takes O ( mn ) flops todecompose a m × n matrix ( m ≥ n ) by Golub-Reinschalgorithm . The SVD routine in GNU Scientific Libraryis used in our numerical calculation. The maximum sizeof matrix ˆ R , and ˆ R ′ are D × D . The SVD of thesetwo matrices takes O ( D ) flops, which take most com-putational complexity in the coarse-graining step, whiledirectly decomposing the D × D matrix A in Eq. (8)takes O ( D ) flops. ∗ [email protected] S. F. Edwards and P. W. Anderson, J. Phys. F Met. Phys. , 965 (1975). G. Parisi, Phys. Rev. Lett. , 1946 (1983). D. S. Fisher and D. A. Huse, Phys. Rev. B , 386 (1988). B. Yucesoy, H. G. Katzgraber, and J. Machta, Phys. Rev.Lett. , 177204 (2012). A. Billoire, L. A. Fernandez, A. Maiorano, E. Marinari,V.Martin-Mayor, G. Parisi, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, and D. Yllanes, Phys. Rev. Lett. , 219701(2013). C. K. Thomas, D. A. Huse, and A. A. Middleton, Phys. Rev. Lett. , 047203 (2011). F. Parisen Toldin, A. Pelissetto, and E. Vicari, Phys. Rev.E , 051116 (2011). K. Tanaka, J.-i. Inoue, and D. Titterington, J. Phys. A , 11023 (2003). H. Derin and H. Elliott, IEEE Trans. Pattern Anal. Ma-chine Intell. , 39 (1987). J. Sun, N.-N. Zheng, and H.-Y. Shum, IEEE Trans. Pat-tern Anal. Mach. Intell. , 787 (2003). S. Z. Li, in
Computer Vision ECCV’94 (Springer, 1994)pp. 361–370. M. Levin and C. P. Nave, Phys. Rev. Lett. , 120601 (2007). U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005). Z.-Y. Xie, H.-C. Jiang, Q.-N. Chen, Z.-Y. Weng, andT. Xiang, Phys. Rev. Lett. , 160601 (2009). M. Hinczewski and A. N. Berker, Phys. Rev. E , 011104(2008). C. G¨uven, M. Hinczewski, and A. N. Berker, Phys. Rev.E , 051110 (2010). C. G¨uven and M. Hinczewski, Phys. A , 2915 (2010). H.-C. Jiang, Z.-Y. Weng, and T. Xiang, Phys. Rev. Lett. , 090603 (2008). Z.-C. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B ,205116 (2008). W. Li, S.-J. Ran, S.-S. Gong, Y. Zhao, B. Xi, F. Ye, andG. Su, Phys. Rev. Lett. , 127202 (2011). Z.-Y. Xie, J. Chen, M.-P. Qin, J.-W. Zhu, L.-P. Yang, andT. Xiang, Phys. Rev. B , 045139 (2012). A. Garc´ıa-S´aez and J. I. Latorre, Phys. Rev. B , 085130(2013). F. Barahona, J. Phys. A , 3241 (1982). J. S. Yedidia, W. T. Freeman, and Y. Weiss, IEEE Trans.Inf. Theory , 2282 (2005). H.-J. Zhou, C. Wang, J.-Q. Xiao, and Z.-D. Bi, J. Stat.Mech. Theor. Exp. , L12001 (2011). H.-J. Zhou and C. Wang, J. Stat. Phys. , 513 (2012). C. Wang and H.-J. Zhou, JPCS , 012004 (2013). A. Lage-Castellanos, R. Mulet, F. Ricci-Tersenghi, andT. Rizzo, J. Phys. A , 135001 (2013). M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, andE. Vicari, Phys. Rev. E , 051115 (2008). Y. Ozeki and H. Nishimori, J. Phys. Soc. Jpn. , 3265(1987). C. K. Thomas and A. A. Middleton, Phys. Rev. E ,043303 (2013). H. Nishimori and K. Nemoto, J. Phys. Soc. Japan , 1198(2002). M. Ohzeki, Phys. Rev. E , 021129 (2009). M. Ohzeki, C. K. Thomas, H. G. Katzgraber, H. Bombin,and M. Martin-Delgado, J. Stat. Mech. Theor. Exp. ,P02004 (2011). T. J¨org and F. Krzakala, J. Stat. Mech. Theor. Exp. ,L01001 (2012). H. Nishimori, Prog. Theor. Phys. , 1169 (1981). K. Takeda, T. Sasamoto, and H. Nishimori, J. Phys. A , 3751 (2005). P. W. Kasteleyn, Physica , 1209 (1961). H. Bethe, Proc. R. Soc. London. Ser. A , 552 (1935). R. Kikuchi, Phys. Rev. , 988 (1951). F. D. Nobre, Phys. Rev. E , 046108 (2001). R. N. Bhatt and A. P. Young, Phys. Rev. B , 5606(1988). F. Parisen Toldin, A. Pelissetto, and E. Vicari, J. Stat.Phys. , 1039 (2009). E. Efrati, Z. Wang, A. Kolan, and L. P. Kadanoff, Rev.Mod. Phys. , 647 (2014). M. M´ezard and G. Parisi, Eur. Phys. J. B , 217 (2001). D. Sherrington and S. Kirkpatrick, Phys. Rev. Lett. ,1792 (1975). M. M´ezard and R. Zecchina, Phys. Rev. E , 056126(2002). J. Chavas, C. Furtlehner, M. M´ezard, and R. Zecchina, J.Stat. Mech. Theor. Exp. , P11016 (2005). G. H. Golub and C. Reinsch, Numerische Mathematik14