aa r X i v : . [ c ond - m a t . s t r- e l ] F e b Quantum Belief Propagation
M. B. Hastings Center for Nonlinear Studies and Theoretical Division,Los Alamos National Laboratory, Los Alamos, NM 87545, [email protected]
We present an accurate numerical algorithm, called quantum belief propagation (QBP), for sim-ulation of one-dimensional quantum systems at non-zero temperature. The algorithm exploits thefact that quantum effects are short-range in these systems at non-zero temperature, decaying on alength scale inversely proportional to the temperature. We compare to exact results on a spin-1 / − , while more elaboratecalculations further reduce the error. The fact that interactions are short-ranged in manyphysical systems is a major simplification in finding quan-tum ground states. The classic example of this is thedensity matrix renormalization group(DMRG)[1], whichrelies on the ability to approximate the ground stateby a matrix product state. While it has long been be-lieved that such an approximation is possible wheneverthere is a spectral gap, due to conformal field theorycalculations[2] and decay of correlations[3, 4], only veryrecently has a general proof been given that such an ap-proximation is possible whenever there is a local Hamil-tonian and a gap[5].At non-zero temperature, the system is in a mixedstate instead. Here, matrix product density operators[6]play the same role that matrix product states do in study-ing pure states, and a non-zero temperature plays a sim-ilar role when it comes to representing mixed states asmatrix product density operators as does a gap for rep-resenting pure states as matrix product states. Indeed, ithas been shown that good matrix product operator rep-resentations exist for quantum systems at any non-zerotemperature[7]. In this paper we present quantum be-lief propagation (QBP), another method for finding ma-trix product density operators for thermal states, whichavoids the problem of Trotter error in other methods.Belief propagation[8] for classical systems (CBP) is es-sentially a Bethe-Peierls solution of a classical statisticalmechanics model. CBP is exact on trees, and is often avery good approximation on lattices with few loops. Atthe end of this paper we will discuss application of theQBP equations to higher dimensional systems and trees,but for now we will focus discussion on belief propagationfor one dimensional lattices. In this case, CBP becomeequivalent to a transfer matrix technique: one solves theproblem on a chain of N -sites to get a partition func-tion which depends on the value of the spin on the N -thsite. Then, one adds a coupling of the N -th spin to oneadditional spin, traces out the N -th spin, arriving at apartition function which depends on the value of the spinon the N + 1-st site. One proceeds in this way, iterativelysolving longer and longer chains.In a quantum system, we will proceed in a similar way.However, now the operator coupling the N -th spin to the N + 1-st spin need not commute with the rest of theHamiltonian. This means that to study properties of thechain of N + 1-spins, it is not sufficient simply to knowthe density matrix of the N -th spin in an N -site chain.However, our physical intuition tells us that at a non-zerotemperature, quantum effects should be short range. Ourmain result in the next section realizes this intuition inthe QBP equations, which involve two terms. One is a“classical” term which is local, coupling the N -th spinto the N + 1-st spin. The other is a “quantum” term,which is non-local, coupling the N + 1-st spin to severalother spins; however, this quantum term is exponentiallydecaying on a length scale set by the inverse tempera-ture. This will allow us to accurately describe statisticalproperties of the N + 1-st spin knowing only the reduceddensity matrix of spins N, N − , ..., N − l + 1, for some l , so that we keep track of a reduced density matrix for l spins. We then iterate these equations by tracing outspin N − l + 1 and then computing statistical propertiesof spin N + 2, keeping always a reduced density matrixon the last l spins on the chain. Quantum Belief Propagation Equations—
TheQBP equations describe how the partition function,exp( − βH ), changes when H is changed by some per-turbation A . We take A, H
Hermitian throughout.We will apply this result to the following case: wehave a nearest-neighbor Hamiltonian, with h i,i +1 theterm coupling spin i to spin i + 1. We will take H = h ( N ) ≡ P N − i =1 h i,i +1 , so that H is the Hamiltonianfor the first N spins, and A = h N,N +1 . We define H s = H + sA. (1)Then, H = h ( N ) + A = h ( N +1) .Although ultimately we want to compute exp( − βH s )at s = 1, we begin by computing the derivative ∂ s exp( − βH s ). Let A ab ( s ) denote matrix elementsof A in a basis of eigenstates of H s , with energies E a ( s ) , E b ( s ). Define the operator A ω,s by its matrix ele-ments: ( A ω,s ) ab = A ab ( s ) δ ( E a ( s ) − E b ( s ) − ω ). We willdo most of the calculations in terms of A ω,s rather than A for simplicity, although we will convert certain resultsback into results in terms of A itself using the integral (cid:16) ∂ s exp( − βH s ) (cid:17) = Z d ω (cid:16) ∂ ǫ exp[ − β ( H s + ǫA ω,s )] (cid:17) , (2)where all derivatives with respect to ǫ are taken at ǫ = 0throughout.One result for the derivative on the right-hand side ofEq. (2) (not the result we will use!) is (cid:16) ∂ ǫ exp[ − β ( H s + ǫA ω,s )] (cid:17) = − exp( − βH s ) (cid:16) exp( βω ) − ω A ω,s (cid:17) . One prob-lem with this is that the operator norm k exp( βω ) − ω A ω,s k may be exponentially large. The QBP equations will im-prove on this, writing ∂ s exp( − βH s ) = η exp( − βH s ) +exp( − βH s ) η † , where k η s k ≤ ( β/ k A k . Finally, η willbe a local operator as discussed below[11].To find the QBP equations, we begin by studying acertain correlation function. Let B be an arbitrary oper-ator. Then, ∂ ǫ tr(exp[ − β ( H s + ǫA ω,s )] B ) (3)= − Z β d τ tr(exp( − βH s ) A ω,s ( − iτ, H s ) B )= − exp( βω ) − ω tr(exp( − βH s ) A ω,s B ) . Adding and subtracting − ( β/ { exp( − βH s ) , A ω,s } B )to Eq. (3) gives ∂ ǫ tr(exp[ − β ( H s + ǫA ω,s )] B ) (4)= − β { exp( − βH s ) , A ω,s } B )+ (cid:16) β e βω ) − e βω − ω (cid:17) tr(exp( − βH s ) A ω,s B )where we used tr( A ω,s exp( − βH s ) B ) =exp( βω )tr(exp( − βH s ) A ω,s B ).We now focus on the correlation functiontr(exp( − βH s ) A ω,s B ), following a procedure very similarto that in [9]. In [9] the result was expressed in terms ofanti-commutators, as A, B were fermionic operators andhence anti-commutated if they were separated in space,while here we will express the result in terms of commu-tators, since we intend to apply it to bosonic operatorswhere the Lieb-Robinson bound is expressed as a boundon the commutator. Using tr(exp( − βH s ) A ω,s B ) = − exp( βω ) tr(exp( − βH s )[ A ω,s , B ]), we have (cid:16) β e βω ) − e βω − ω (cid:17) tr(exp( − βH s ) A ω,s B )= βF ( ω )tr(exp( − βH s )[ A ω,s , B ]) , (5)where F ( ω ) ≡
12 1 + e βω − e βω + 1 βω = − coth( βω/ βω . (6) Thus for any operator B , we have ∂ ǫ tr(exp[ − β ( H s + ǫA ω,s )] B ) = − β tr( { exp( − βH s ) , A ω,s } B ) + βF ( ω )tr(exp( − βH s )[ A ω,s , B ]), and hence ∂ ǫ exp[ − β ( H s + ǫA ω,s )] (7)= − β { exp( − βH s ) , A ω,s } + βF ( ω )[exp( − βH s ) , A ω,s ] . We define η ωs = − ( β βF ( ω )) A ω,s , (8) η s = Z d ωη ωs and integrate Eq. (7) over ω to get ∂ s exp( − βH s ) = η s exp( − βH s ) + exp( − βH s ) η † s . (9)Eq. (9) is the QBP equation for ∂ s exp( − βH s ). Theanti-commutator in Eq. (7) is a “classical term”. It is theonly term present if [ A, H ] = 0, in which case these equa-tion reproduce the classical belief propagation equations.The commutator in Eq. (7) is a “quantum” term. It isodd in ω and vanishes at ω = 0. We now discuss local-ity properties of the quantum term, assuming that H islocal in the sense of having a Lieb-Robinson bound[10].We have[9] β R d ωF ( ω ) A ω,s = Z d ω (cid:16) − β βω/
2) + 1 ω (cid:17) A ω,s (10)= − Z d ω X n =0 ( ω − πni/β ) − A ω,s = i X n ≥ Z ∞−∞ d t sign( t ) exp( − πnt/β ) A ( t, H s ) , where the sum is over integer n and A ( t, H s ) =exp( iH s t ) A exp( − iH s t ). The integral over t is exponen-tially decaying for t & β and so, using a Lieb-Robinsonbound, η is local. Numerical Implementation and Results—
In this sec-tion we discuss the numerical implementation of the QBPequations, and the results of their application to the an-tiferromagnetic spin-1 / A and H = h ( N ) , and instead set H = P N − i = N − l +2 h i,i +1 , forsome constant l ≥
2. Since η is local, this approxima-tion is justified for small enough β ; as β increases, l must increase and the numerical effort is exponential in l , since one must diagonalize matrices of size 2 l .Although translation invariance is not necessary, itis useful as we can then apply translations to H = P N − i = N − l +2 H i,i +1 to make H equal to h ( l − , so forthe Heisenberg chain with coupling constant J = 1, H = P i, ≤ i ≤ l − ~S i · ~S i +1 . We set A = h l − ,l . Weset T to be the operator which translates one site to theright. We define ρ to be a reduced density matrix forsites 1 ...l , so that ρ is a 2 l dimensional matrix. Define O = S ′ exp( Z η s ′ d s ′ ) , (11)where the S ′ denotes that the integral is s ′ -ordered.The algorithm to compute the free energy per site of aninfinite chain proceeds through the following four steps.(1) Initialize ρ to exp( − βH ). (2) Approximately com-pute the operator O as described below. (3) Go througha series of n it iterations of the following three steps:(A) replace ρ with OρO † . (B) trace out the first site, sothat ρ is replaced with tr ( ρ ) ⊗ l +1 . Here, tr denotesa partial trace over the first site, and l +1 is the unitoperator on the l + 1-th site. At this stage ρ is now adensity operator on sites 2 ...l + 1. Translate by one site,so that ρ becomes a density operator on sites 1 ...l again.(C) Define Z = tr( ρ ) and then replace ρ with Z − ρ .After these iterations, we (4) Output ln( Z ) from thelast step as the free energy. This procedure relies on aseries of iterations of steps ( A ) , ( B ) , ( C ) to find the ρ which is the fixed point of the map ρ → Z − T † (cid:16) tr ( OρO † ) ⊗ l +1 (cid:17) T. (12)The number of iterations required for convergence ap-pears to increase roughly linearly with l and β , but thecomputational effort grows only linearly in the numberof iterations. The locality of the operators η, O justifiestracing out the first site in step (B). After n it iterations, ρ is approximately proportional to the reduced density ma-trix tr ...n it (exp( − βh ( l − n it ) )), where the partial traceis over sites 1 ...n it .To compute correlation functions of operators such as S zi S zi + j we follow the following procedure. We go throughsteps (1) , (2) , (3) as above, with n it iterations in step (3).On the last of the n it iterations, after step (A) we copy ρ to a new matrix ρ corr . We then replace ρ corr with S z ρ corr , thus inserting the first of the two operators. Onsteps (B,C) of the last iteration we replace of the map ρ → Z − tr ( ρ ) ⊗ ρ corr → Z − tr ( ρ corr ) ⊗
1, wherethe same Z is used for both ρ and ρ corr . When then pro-ceed through j more iterations of steps (A),(B),(C), andon the last iteration after step (A) we replace ρ corr with S z ρ corr , inserting the second of the two operators, beforeproceeding with steps (B),(C). We then proceed throughseveral more iterations in which at each step we map ρ → Z − tr ( OρO † ) ⊗ ρ corr → Z − tr ( Oρ corr O † ) ⊗ ρ corr ) / tr( ρ ).To compute the matrix O , we approximate by dividingthe integral over s ′ into n slice different slices: O ≈ exp( η s ( n slice ) n slice ) ... exp( η s (2) n slice ) exp( η s (1) n slice ) , (13)where s ( m ) = ( m − / /n slice . To computethe matrix exponential for each slice, exp( η s ′ ) for s ′ = (1 / /n slice ) , (3 / /n slice ) , (5 / /n slice ) , ... we used a Taylor series method, while to compute η s ′ itself we diagonalize H s ′ and transform A into a basis ofeigenvector of H s ′ . We use the fact that H and A con-serve total S z to speed up both this diagonalization andthe multiplication ρ → OρO † .Since η s is non-Hermitian, but || η s || is not too large,the Taylor series method is a good choice for computingthe matrix exponential. Eq. (13) approximates the s ′ -dependence of η s ′ by taking its value halfway through theslice, which gives us first order accuracy in ∂ s ′ η s ′ for free.Further, η s ′ is in fact only very weakly dependent on s ′ and so a small n slice suffices, as seen by the following test:set l = 3. Then, after one iteration of steps (A),(B) andbefore normalizing in step (C), the matrix ρ should beequal to, the thermal density matrix for a three site chain.At β = 1, the largest eigenvalue should be equal to e . Fora calculation with n slice = 1, the largest eigenvalue wasfound to be 2 . n slice = 2, we find 2 . n slice = 3 we find 2 . β P j h S zi S zi + j i , at the susceptibility peak, using theknown location[12] of the peak at T max = 1 /β =0 . χ exact ( T max ) = 0 . ... (14)while a calculation using l = 5, n it = 20, n slice = 20,and correlations up to j = 20 gives χ QBP ( T max ) = 0 . ... (15)for a relative error of ≈ ∗ − . The calculation took . .
08 seconds on a 1.5 GHz PowerPC G4 processor,and using conservation of S z the largest matrix diago-nalized was 10-by-10. A larger calculation, with l = 9, n slice = 50, n it = 30 improves this to χ QBP ( T max ) =0 . ... for a relative error of ≈ ∗ − .We calculated the specific heat C as a function of tem-perature in two ways: first, by calculating β ∂ β ln( Z ) us-ing ln( Z ) from the algorithm, and, second, by calculating − β ∂ β h S zi S zi +1 i . The results are shown in Fig. 1, whereto take derivatives we calculated ln( Z ) and h S zi S zi +1 i for β = 0 . , . , . , ..., .
0. As l gets larger, the curves re-main accurate to lower temperature. The peak specificheat for l = 7 was 0 . ... from the second derivativecalculation and 0 . ... from the first derivative, bothof which compare very well with the Bethe ansatz resultof 0 . ... .The accuracy can be improved by going to larger l . Another improvement is to take H = h ( l − +(1 / h l − ,l − and A = (1 / h l − ,l − + h l − ,l ) in-stead of H = h ( l − and A = h l − ,l . We still have H + A = T † HT + h , in this case, but the slightly differ-ent form of the perturbation seems to work better. Thefigure inset shows a comparison of Bethe ansatz data to Temperature S p ec i f i c h ea t FIG. 1: Specific heat against temperature for l = 3 (dashedline),5 (dotted line), 7 (solid line). Curves that go negativeare from − β ∂ β h S zi S zi +1 i , while those that diverge positivelyare from β ∂ β ln( Z ). Inset: l = 9 and Bethe ansatz. l = 9 (where the largest matrix diagonalized is 126 di-mensional). Discussion—
The implementation of the QBP equa-tions here must be considered as preliminary. More workis needed to optimize the algorithm and, most impor-tantly, to quantify the sources of error. Despite this,the method yields accurate results, giving qualitativelycorrect behavior even at l = 3 where QBP can be im-plemented as a “pen-and-paper” technique.In contrast to QBP, thermodynamic DMRG[13] com-putes low-lying eigenvalues for finite size chains up to size N = 30; the results were accurate to quite low temper-atures, but required an extrapolation to avoid finite sizeeffects. In that work, chains of size up to N = 14 wereexactly diagonalized, so if we improve the linear algebraroutines used in our implementation it should certainlybe possible to do QBP with an l ∼
14. It may be pos-sible to do QBP with l ∼
30 using DMRG techniquesto avoid keeping all of the eigenstates of H and insteadtruncate to some smaller subset of eigenstates of H . Ex-trapolation of our results suggests that this should permitaccess to temperatures T /J ∼ / l once l becomes of or-der Jβ . We now describe a procedure that combinessome of the ideas of QBP and transfer matrix DMRG. Introduce M copies of the system, each with density ma-trix exp( − βh ( N ) /M ), so that the joint density matrixis exp( − βh ( N ) /M ) ⊗ ... ⊗ exp( − βh ( N ) /M ). Let P i bethe operator that cyclically permutes the value of thespin on site i between the M different copies. Then,tr( P P ...P N exp( − βh ( N ) /M ) ⊗ ... ⊗ exp( − βh ( N ) /M )) =tr exp( − βh ( N ) ). Let ρ be a 2 l M dimensional matrix.We will define a QBP procedure such that after n it it-erations, ρ is approximately proportional to the reduceddensity matrix tr ...n it ( P ...P n it exp( − βh ( l − n it ) /M ) ⊗ ... ⊗ exp( − βh ( l − n it ) /M )). Define O by Eq. (11), for H = h ( l − at inverse temperature β/M , and map ρ → Z − T † (cid:16) tr ( P ( O ⊗ ... ⊗ O ) ρ ( O † ⊗ ... ⊗ O † )) ⊗ l +1 ⊗ ... ⊗ l +1 (cid:17) T. (16)For M = 1, this reduces to the QBP implementationdescribed here, while for l = 2, this becomes very similarto the transfer matrix used in transfer matrix DMRG.The question is whether for M > , l > ρ of this transfer matrix.QBP can be directly applied to finite size chains andto infinite or finite trees. The ability to compute real-space correlation functions and handle translationallynon-invariant systems are advantages of this method, andin future this method will be applied to disordered sys-tems where TDMRG will have problems. Probably themost interesting question is the possible application ofQBP to higher dimensional systems, by replacing thehigher dimensional lattice with a Cayley tree or Husimicactus with the correct local structure. Acknowledgments—
I thank M. Chertkov for introduc-ing me to classical belief propagation and A. Kl¨umper forsupplying the Bethe ansatz data. This work supportedby U. S. DOE Contract No. DE-AC52-06NA25396. [1] S. R. White, Phys. Rev. Lett. , 2863 (1992).[2] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys.Rev. Lett. , 227902 (2003).[3] M. B. Hastings, Phys. Rev. Lett. , 140402 (2004).[4] M. B. Hastings, Phys. Rev. B , 104431 (2004).[5] M. B. Hastings, arXiv:0705.2024.[6] F. Verstraete, J. J. Garc´ıa-Ripoll, and J. I. Cirac, Phys.Rev. Lett. , 207204 (2004); M. Zwolak and G. Vidal,Phys. Rev. Lett. , 207205 (2004).[7] M. B. Hastings, Phys. Rev. B , 085115 (2006).[8] R. G. Gallager, Low density parity check codes (MITPress Cambridge, MA, 1963).[9] M. B. Hastings, Phys. Rev. Lett. , 126402 (2004).[10] E. H. Lieb and D. W. Robinson, Commun. Math. Phys. , 251 (1972); M. B. Hastings and T. Koma, Commun.Math. Phys. , 781 (2006); B. Nachtergaele and R.Sims, Commun. Math. Phys. , 119 (2006). [11] For small enough β , the operator R d ω exp( βω ) − ω A ω,s islocal as may be shown by a power series expansion, butat larger β this fails.[12] A. Kl¨umper and D. C. Johnston, Phys. Rev. Lett. ,4701 (2000); M. Shiroishi and M. Takahashi, Phys. Rev.Lett. , 117201 (2002).[13] S. Moukouri and L. G. Caron, Phys. Rev. Lett. , 4640(1996).[14] X. Wang and T. Xiang, Phys. Rev B56