Loop Series Expansions for Tensor Networks
Glen Evenbly Nicola Pancotti Ashley Milsted Johnnie Gray Garnet Kin-Lic Chan




evenbly@amazon.com


AWS Center for Quantum Computing, Pasadena, CA 91125, USA
California Institute of Technology, Pasadena, CA 91125, USA


California Institute of Technology, Pasadena, CA 91125, USA


§ ABSTRACT 
Belief propagation (BP) can be a useful tool to approximately contract a tensor network, provided that the contributions from any closed loops in the network are sufficiently weak. In this manuscript we describe how a loop series expansion can be applied to systematically improve the accuracy of a BP approximation to a tensor network contraction, in principle converging arbitrarily close to the exact result. More generally, our result provides a framework for expanding a tensor network as a sum of component networks in a hierarchy of increasing complexity.  We benchmark this proposal for the contraction of iPEPS, either representing the ground state of an AKLT model or with randomly defined tensors, where it is shown to improve in accuracy over standard BP by several orders of magnitude whilst incurring only a minor increase in computational cost. These results indicate that the proposed series expansions could be a useful tool to accurately evaluate tensor networks in cases that otherwise exceed the limits of established contraction routines.


Loop Series Expansions for Tensor Networks
    Garnet Kin-Lic Chan
    September 9, 2024
==========================================






§ INTRODUCTION
 
Tensor networks (TNs) are powerful tools for the classical simulation of quantum systems <cit.>. Although they were traditionally developed for model systems in the context of condensed matter physics, they have more recently been applied towards the benchmarking of nascent quantum computers <cit.>. Tensor network states were pioneered by the introduction of matrix product states <cit.> (MPS), which possess tensors connected in a 1D geometry. For several decades numeric algorithms based on the MPS, such as the density matrix renormalization group<cit.> (DMRG), have been widely regarded as the gold-standard of classical simulation algorithms for 1D quantum systems; producing high-accuracy and unbiased results for a wide range of problems. A large part of this success can be attributed to the computational efficiency with which MPS can be contracted. Subsequent developments have seen the proposal of sophisticated tensor networks, such as projected entangled pair states<cit.> (PEPS) and multi-scale entanglement renormalization ansatz<cit.> (MERA), designed to efficiently represent quantum states on large 2D systems. However, a significant challenge of these 2D tensor networks is the high computational expense required for their contraction <cit.> even when approximate contraction schemes are used <cit.>. In practice, this expense limits the amount of entanglement that can be encoded in these networks and therefore also limits the range of problems for which they produce accurate results. Thus the development of strategies for more efficient contraction of TNs remains a vital task.

A promising approach for the contraction of TNs is through the adoption of belief propagation<cit.> (BP), a message passing algorithm with a long history in both computer science<cit.> and in statistical physics<cit.>. In recent times the application of BP towards TNs has generated much interest <cit.>, and produced a notable success <cit.> of simulating a quantum experiment performed on IBMs 127-qubit Eagle processor <cit.>, a problem that was seen to be intractable for other 2D TN algorithms <cit.>. However, a major drawback of BP is that it is an uncontrolled approximation; one cannot systematically improve the accuracy of a BP result as is otherwise done in the context of tensor networks by simply growing the bond dimension, although there do exist schemes to improve the accuracy of BP through a preliminary clustering of a network <cit.>. 

In this manuscript we present a method to improve upon the accuracy of BP for the contraction of a tensor network by incorporating and building upon the loop series expansion proposed by Chertkov and Chernyak<cit.>. Assuming that a BP fixed point of the tensor network has been found, we describe how the dominant loop excitations, otherwise neglected by BP, can be re-incorporated into the evaluation of this network. This approach can also be understood as providing a general notion of a series expansion for tensor networks; that a complicated tensor network may be expanded as a sum of component networks in a hierarchy of increasing complexity, with the BP approximation as the zeroth-order term in the series. We demonstrate numerically that the proposed series expansion can provide both accurate and computationally cheap evaluations of PEPS, such that it could be a useful tool to contract tensor networks in cases that otherwise exceed the limits of established contraction routines.




§ BELIEF PROPAGATION
 
Here we provide a brief introduction to BP as a necessary precursor to our main results. Let 𝒯 be a closed tensor network defined of N tensors labelled {T_1, T_2, T_3, …, T_N } and with M indices labelled {i_1, i_2, i_3, … i_M }. Assume that our goal is to contract 𝒯 to a scalar 𝒵(𝒯),

    𝒵(𝒯) = ∑_i_1,i_2,… T_1 (ξ_1)  T_2 (ξ_2)  T_3 (ξ_3) …,

where ξ_r describes the subset of indices ξ_r ⊆{i_1, i_2, … i_M } connected to tensor T_r. The starting point of BP is to define a pair of messages for each edge of the network; if (r,s) denotes an edge shared between tensors T_r and T_s then we define messages {μ_r→ s, μ_s→ r} on this edge. We refer to a message μ_s→ r as outgoing from T_s and incoming to T_r. BP prescribes the following algorithm to update each of the messages: the new outgoing message for T_r towards tensor T_s is given by contracting the incoming messages to T_r on all other indices,

    μ_r→ s = ( ∏_q ∈ξ_r / sμ_q→ r) T_r.

The set of 2M messages realize a BP fixed point if all of the equations are simultaneously satisfied up to multiplicative constants. Let us assume that we have found a BP fixed point for our network 𝒯 and that each matching pair of messages have been normalized such that their dot product is unity, ⟨μ_r→ s·μ_s→ r⟩ = 1. Let us define scalar T̃_r as resulting from the contraction of a tensor T_r with all incoming messages,

    T̃_r = ( ∏_q ∈ξ_rμ_q→ r) T_r,

which we refer to as the (BP) vacuum contribution from tensor T_r. The BP approximation to the network contraction 𝒵 from Eq. <ref> is given by multiplication of these scalars

    𝒵(𝒯) ≈∏_rT̃_r.

It follows that the free energy ℱ = -log𝒵 can be approximated as

    ℱ(𝒯) ≈ - ∑_r log( T̃_r ),

which is known as the Bethe free energy <cit.>.








§ CONFIGURATIONS IN THE BP BASIS
 
We now outline a methodology to improve upon the accuracy of a known BP fixed point by introducing additional terms into Eq. <ref>. For any edge (r,s) of the network, assumed to be of dimension d, we define a rank 1 projector 𝒫_rs onto message subspace and its rank (d-1) compliment 𝒫_rs^C,

    𝒫_rs = μ_s→ r⊗μ_r→ s,     𝒫_rs^C = 𝕀 - 𝒫_rs.

with 𝕀 as the d-dim identity, see also Fig. <ref>(a).
We henceforth refer to 𝒫_rs as the projector onto the (BP) ground state of edge (r,s) and 𝒫_rs^C as the projector onto the excited subspace. It follows that the tensor network 𝒯 can be resolved as a sum over 2^M (BP basis) configurations, formed from all combinations of projecting each edge into either the ground or excited sub-space. We define the degree of a configuration as the number of excited edges that it contains, using δ_x to denote a degree-x configuration and W(δ_x) to denote its weight (i.e. the scalar resulting from its contraction). Thus we can recognise δ_0 as the BP vacuum; the unique configuration with all edges projected into their BP ground state. It follows that W(δ_0) evaluates to the standard BP result of Eq. <ref>. Our goal is now to characterise the weights of the remaining configurations, which may be thought of as excitations on top of the BP vacuum, such that the high-weight contributions can then be re-introduced into the approximation. 

A significant step towards this goal comes from the understanding derived in Ref. LOOP1 that any configuration with a dangling excitation, i.e. that contains a tensor with a single excited index, has weight zero. This result follows directly from the definition of the BP fixed point; let us assume that edge (r,s) of tensor T_r is projected via 𝒫^C_rs into the excited subspace and all other edges are in the BP ground. It follows that

    ( ∏_q ∈ξ_r / sμ_q→ r) T_r  𝒫^C_rs = μ_r→ s𝒫^C_rs = 0,

where Eq. <ref> was used to replace T_r with the outgoing message μ_r→ s which is, by construction in Eq. <ref>, orthogonal to the projector 𝒫^C_rs. This result, which implies that non-trivial excitations can only occur in closed loops, is useful in reducing the number of basis configurations that need be considered. For instance, when applied to the network of M=7 edges in Fig. <ref>(b), the number of configurations is reduced from 2^7=128 to only 5 configurations of non-zero weight labelled {δ_0, δ_4, δ_4', δ_6, δ_7 }. Another useful observation, as examined in Ref. LOOP3, is that the weight of a configuration containing separate excitations, i.e. excitations that do not share any common tensors, is given by the product of the individual weights

    W(δ_x ⊗δ_y) = W(δ_x) × W(δ_y),

see also the example shown in Fig. <ref>(c). This result is further simplified by re-scaling each tensor T_r so that its vacuum state contribution, i.e. T̃_r from Eq. <ref>, is unity which allows the vacuum component of any configuration to be ignored. This observation implies that one only needs to understand the connected excitations in order to fully characterize a tensor network contraction.

For any reasonably large tensor network the full number of non-trivial excitations will still exceed practical consideration. However, as detailed further in Appendix <ref>, we posit that in most spatially homogeneous tensor networks the weight of a configuration δ_x will be exponentially suppressed in its degree x, 

    W(δ_x) ≈exp (-k x)

with k a positive scalar. While it is easy to construct counterexamples where this scaling will not hold, it is argued in Appendix <ref> that this behavior will hold for networks where the BP fixed point already provides a reasonably accurate approximation to the network contraction. Thus one can understand that the (BP basis) configurations, when ordered in increasing degree, form a series expansion for a network contraction. While the number of δ_x configurations may still grow (exponentially) quickly with x, the exponential suppression of the contribution of each configuration provides the possibility that a partial summation will rapidly converge to the exact result as higher-degree excitations are added. 



§ CONTRACTION OF IPEPS
 
In this section we discuss how the proposed series expansion can be applied to contract infinite projected entangled pair states (iPEPS). For concreteness we focus on a hexagonal lattice PEPS |ψ⟩, local dimension d and virtual dimension m, formed from a single pair of unique tensors A and B tiled on alternate sites (such that the state has a 2-site unit cell). A closed hexagonal lattice tensor network 𝒯 = ⟨ψ|ψ⟩ of bond dimension m^2 is formed by taking the product of the PEPS with its conjugate and contracting over the site indices, see also Fig. <ref>(a). The network 𝒯 is composed of two unique tensors a and b, which result from contraction of the physical index in (A⊗ A^*) and (B⊗ B^*) respectively. In what follows we assume that we have found a BP fixed point for 𝒯 and that tensors a and b have been re-scaled such that their vacuum contributions, as defined in Eq. <ref>, are unity.







 §.§ Free-energy density
 
A common task in tensor network algorithms is to evaluate the scalar 𝒵 associated to contracting a closed tensor network, as per Eq. <ref>. However in the context of an infinite tensor network this scalar would diverge so we instead evaluate an intensive equivalent: the free-energy density f,

    f = lim_N →∞( -log(𝒵_N)/N).

where 𝒵_N represents the evaluation of a patch of N tensors from the full network. Note that, due to our normalization of the BP vacuum, the BP contribution to 𝒵_N is unity for all N (or equivalently the Bethe free energy is zero <cit.>).

Our first step is to evaluate the weights W(δ_x) of the lowest degree excitations δ_x; in the calculations performed in this manuscript we evaluate the excitations up to degree x=12 as shown in Fig. <ref>(b-e). Each of these weights is given by contracting a small tensor network formed by capping the external indices (i.e. of tensors within the excitation) with their fixed point messages and projecting the internal indices onto their excited subspace with projectors 𝒫_rs^C as defined in Eq. <ref>. However, even once these weights W(δ_x) are known, it still remains a non-trivial combinatorial task to compute their impact on the free energy density f. An approximate counting of the possible combinations of loop excitations, as derived in Appendix <ref>, gives the following result for the free energy correction

    f ≈ -∑_{δ_x} L(δ_x) W(δ_x) e^S(δ_x) f.

with L(δ_x) as the number of locations per lattice site that excitation can be placed and S(δ_x) as the number of lattice sites that the excitation occupies. For instance, the δ_6 excitation shown in Fig. <ref>(b) has L(δ_6) = 1/2 (since there is one hexagon per two lattice sites) and S(δ_6) = 6. Notice that the terms exp(S(δ_x) f) in Eq. <ref> act to suppress excitations based on the amount of space they occupy. Given that Eq. <ref> contains the unknown f on both sides we cannot solve for f directly. Instead we solve Eq. <ref> iteratively by first setting f=0 on the r.h.s, computing a new value for f, then feeding the new f back into r.h.s and repeating until f is converged. Usually this requires only a few iterations to converge to high precision. 



 §.§ Transfer matrix
 
We now examine how one can compute the transfer matrix T_AB associated to cutting open a single edge from the network 𝒯 = ⟨ψ|ψ⟩ as shown in Fig. <ref>(a). Such transfer matrices are utilized in determining the optimal truncation of an internal index from the PEPS and thus play an important role in many PEPS optimization algorithms <cit.>.

The BP vacuum contribution to T_AB is shown in Fig. <ref>(b) together in Fig. <ref>(c) with examples of some of the contributions from loop excitations, all of which individually evaluate to a matrix with the same dimensions as T_AB. The series approximation to T_AB is constructed by evaluating each of the loop-excitation networks (up to some desired degree), then adding the contributions together after weighting each by a suppression factor exp(S(δ_x)f) similar to Eq. <ref>. Note this requires that the free energy density f be evaluated prior to computing T_AB. Finally, we remark that excitation terms that do not include either of the tensors adjacent to an open index do not need to be considered explicitly; their contributions, which are proportionate to the vacuum term, are already accounted for through the suppression factors. 







 §.§ Density matrix
 
We now describe how the density matrix ρ_AB associated to a pair of adjacent lattice sites can be computed. Let us first define tensor a' from taking (A⊗ A^*) while leaving the physical indices uncontracted and similarly for b', see Fig. <ref>(e). The network for ρ_AB can be realized from a substitution of a' and b' as impurities in the original network 𝒯 = ⟨ψ|ψ⟩ as shown in Fig. <ref>(d).

Similarly to the transfer matrix, only the loop excitations that include one (or both) of the tensors a' and b' will give non-trivial contribution to ρ_AB. The BP vacuum contribution to ρ_AB is shown in Fig. <ref>(f) together in Fig. <ref>(g) with an example of a loop excitation. The series approximation to ρ_AB is thus constructed in the same way as the transfer matrix; the contributing loop-excitation networks are evaluated up to some desired degree and then added together after weighting each by a suppression factor exp(S(δ_x)f). 



§ BENCHMARK RESULTS
 
We now benchmark our approach, beginning with contraction of a PEPS representing the ground state |ψ⟩ of the hexagonal lattice AKLT model<cit.> in the thermodynamic limit. The AKLT model is particular useful as a test bed, since |ψ⟩ can be represented exactly<cit.> with an iPEPS of bond dimension m=2. Following the methodology prescribed in Sect. <ref> we compute the free energy density f, transfer matrix T_AB, and the density matrix ρ_AB for loop expansions containing up to 12^th degree excitations δ_12 shown in Fig. <ref>. In all cases we compare against evaluations using a boundary MPS contraction<cit.> whose dimension χ has been chosen large enough, χ≈ 30, such that any truncation errors are negligible and the results can be considered as numerically exact. Fig. <ref> shows that the addition of the loop corrections improve the accuracy of the BP results dramatically, such that they appear to converge towards the exact results as the degree of the expansion is increased. The expansion to 12^th degree excitations δ_12 gave the free energy density f accurate to ε = 2×10^-7, which represents an improvement of four orders of magnitude over the BP result. That the infinite hexagonal tensor network can be characterised so accurately from only the BP vacuum plus a few of the low-degree loop excitations is remarkable. Note that we also compare against results obtained from the exact contraction of the PEPS on a finite N× N geometry (with periodic boundaries) for N={6,8,10,12 }. In all cases we see that an N^th degree loop expansion exceeds the accuracy of an exact evaluation on a N× N geometry. This is expected since the N× N systems would contain degree-N corrections (i.e. that `wrap' the whole system) whereas a degree-N loop expansion is only neglecting excitations of degree larger than N. None-the-less, this provides additional confirmation that the loop expansions are properly accounting for all short-range contributions.

In Fig. <ref>(a) we present results for the accuracy of the free-energy density evaluated from randomly defined iPEPS (with 2-site unit cell) of local dimension d and bond dimension m. These were prepared by placing on each edge of the lattice a maximally entangled pair of dim-m particles and then projecting each vertex onto a dim-d space via a (m^3× d) isometry. The isometries were formed by first populating a (m^3× d) matrix with random elements chosen uniformly from the interval [0,1], then enforcing rotational symmetry (useful to reduce the number of unique terms in the loop expansion), and finally orthogonalizing the columns of the matrix. The results from random iPEPS, which are broadly similar to those from the AKLT model, demonstrate that the series expansion seems to be valid for typical instances of iPEPS. Notice that increasing the local dimension d (while keeping a fixed bond dimension m) was generally found to reduce the accuracy of our approach; this could be expected as increasing d also increases the per-site entanglement entropy of the random PEPS.

Fig. <ref>(b-c) presents equivalent results for kagome and square lattices. For iPEPS defined on kagome lattices we perform local manipulations to exactly transform the network into a hexagonal geometry (similar to the form of PEPS described in Ref. PESS1) prior to applying the series expansion, see Appendix <ref> for details. The square lattice results were obtained by applying the expansion directly to the square lattice, evaluated up to the δ_8 excitations as shown in Appendix <ref>. Similar to the hexagonal lattice, the results from the kagome and square lattices also show robust improvements over the BP approximation. However we see that the convergence on the square lattice is slower with respect to increasing the degree of expansion; we attribute this to the fact that the quantity of distinct loop configurations at any given degree grows much faster than on the hexagonal lattice. 











§ DISCUSSION
 
A key result presented in this manuscript is a systematic way to improve on BP results for the contraction of tensor networks, turning BP from an uncontrolled approximation into a well-controlled approximation. This result would be directly useful in contexts where BP has been employed for the contraction of tensor networks, such as the simulation of quantum experiments<cit.>.

Another direct application of these results could be for PEPS optimization algorithms. Given that a recent work<cit.> established the equivalence of BP to the `simple update' strategy for PEPS, the series expansion could significantly enhance the accuracy of the `simple update' strategy while still remaining relatively computationally efficient. For instance, the genus-1 loop excitations on the hexagonal lattice, i.e. the δ_6 and δ_10 terms in Fig. <ref>, can be evaluated with O(m^6) cost for PEPS of bond dimension m, and this cost could also potentially be reduced via approximation strategies like sparse diagonalization. Thus the dominant terms in the expansion are still feasible to compute even for PEPS with bond dims in excess of m=100, well beyond the computational limits of boundary MPS or corner transfer matrix methods for PEPS contractions<cit.>. 

We also expect the series expansions to be useful in the context of coarse-graining tensor algorithms for tensor networks, especially in relation to proposals that extend the tensor renormalization group<cit.> (TRG) by including mechanisms to identify and remove short-range loop correlations<cit.>. Given that our results provide a general framework precisely for efficiently identifying loop correlations it is likely that they would find significant utility in this setting. 

Although our proposal worked well for the test cases considered, it is clear that there are also several potential ways that it could fail. Most obviously, application of the series expansion first requires finding a BP fixed point. It is unknown how difficult this task is for a PEPS representing the ground state of a local Hamiltonian or other physically relevant tensor networks. Even if a BP fixed point is found, the series expansion would not be expected to converge if the network under consideration possessed multiple degenerate, or almost degenerate, fixed points. In this scenario Eq. <ref> would fail; relative to any single BP fixed point there would exist high-degree configurations that still carry large weight (i.e related to a different BP fixed point). An example of this would be a PEPS encoding of a GHZ state, where a BP fixed point could realize either the |↑↑↑…⟩ or |↓↓↓…⟩ component of the network; relative to one of these fixed points the other component would represent a maximal-degree excitation with equal weight.  

Most existing tensor network algorithms are built on the twin pillars of tensor contractions and tensor decompositions/truncations. The results presented in this manuscript, which prescribe a general procedure for resolving a complicated network into its fundamental components, could represent a third pillar for tensor network algorithms; one that could have interesting synergies with established routines. For instance, in evaluating a network 𝒯, one could begin by expanding in terms of the BP vacuum, the short-loop excitations, and the remaining `residual'. Potentially, one could then apply a traditional (approximate) contraction strategy on the residual which, presumably, would be more amenable to truncations due to the separation from the short-loop excitations. Exploring the ways that series expansions could compliment or enhance existing tensor network algorithms remains an interesting avenue for future research. 


We thank AWS for supporting the quantum computing program. GKC acknowledges support from the US DOE, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator (QSA), and from a generous gift from AWS. 


99





TN1
J. I. Cirac and F. Verstraete, Renormalization and tensor product states in spin chains and lattices, J. Phys. A: Math. Theor. 42, 504004 (2009).

TN2
T.G. Kolda and B.W. Bader, Tensor decompositions and applications, SIAM review, 51(3), pp.455-500 (2009).

TN3
G. Evenbly and G. Vidal, Tensor network states and geometry, J. Stat. Phys. 145, 891-918 (2011).

TN4
R. Orus, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. 349, 117 (2014).

TN5
J. C. Bridgeman and C. T. Chubb, Hand-waving and Interpretive Dance: An Introductory Course on Tensor Networks, J. Phys. A: Math. Theor. 50, 223001 (2017).

TN6
S. Montangero, Introduction to Tensor Network Methods - Numerical Simulations of Low-dimensional Many-body Quantum Systems, (Springer, Berlin, 2018).

TN7
R. Orus, Tensor networks for complex quantum systems, Nat. Rev. Phys. 1, 538-550 (2019).

TN8
P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Phys. Lect. Notes 8 (2019).

TN9 
S.-J. Ran, E. Tirrito, C. Peng, X. Chen, L. Tagliacozzo, G. Su, and M. Lewenstein, Tensor Network Contractions Methods and Applications to Quantum Many-Body Systems, (Springer, 2020).

TN10
J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).





QC1
E. S. Fried, N. P. D. Sawaya, Y. Cao, I. D. Kivlichan, J. Romero, and A. Aspuru-Guzik, qTorch: The quantum tensor contraction handler, PLoS ONE 13(12): e0208510. (2018).

QC2
B. Villalonga, S. Boixo, B. Nelson, C. Henze, E. Rieffel, R. Biswas and S. Mandrà, A flexible high-performance simulator for verifying and benchmarking quantum circuits implemented on real hardware, npj Quantum Inf 5, 86 (2019).

QC3
R. Schutski, D. Lykov, and I. Oseledets, Adaptive algorithm for quantum circuit simulation, Phys. Rev. A 101, 042335 (2020).

QC4
Y. Zhou, E. M. Stoudenmire, and X. Waintal, What Limits the Simulation of Quantum Computers?, Phys. Rev. X 10, 041038 (2020). https://doi.org/10.1103/PhysRevX.10.041038DOI: 10.1103/PhysRevX.10.041038. 

QC5
M. Levental, Tensor Networks for Simulating Quantum Circuits on FPGAs, arXiv:2108.06831 (2021).

QC6
F. Pan and P. Zhang, Simulation of Quantum Circuits Using the Big-Batch Tensor Network Method, Phys. Rev. Lett. 128, 030501 (2022)

QC7
T. Vincent, L. J. O'Riordan, M. Andrenkov, J. Brown, N. Killoran, H. Qi, and I. Dhand, Jet: Fast quantum circuit simulations with parallel task-based tensor-network contraction, arXiv:2107.09793 (2021).

QC8
F. Pan and P. Zhang, Simulation of Quantum Circuits Using the Big-Batch Tensor Network Method, Phys. Rev. Lett. 128, 030501 (2022). https://doi.org/10.1103/PhysRevLett.128.030501DOI: 10.1103/PhysRevLett.128.030501. 

QC9
F. Pan, K. Chen, and P. Zhang, Solving the Sampling Problem of the Sycamore Quantum Circuits, Phys. Rev. Lett. 129, 090502 (2022). https://doi.org/10.1103/PhysRevLett.129.090502DOI: 10.1103/PhysRevLett.129.090502. 

QC10
T. Ayral, T. Louvet, Y. Zhou, C. Lambert, E. M. Stoudenmire, and X. Waintal. Density-Matrix Renormalization Group Algorithm for Simulating Quantum Circuits with a Finite Fidelity, PRX Quantum 4, 020304 (2023). https://doi.org/10.1103/PRXQuantum.4.020304DOI: 10.1103/PRXQuantum.4.020304. 

QC11
J. Tindall, M. Fishman, M. Stoudenmire, and D. Sels. Efficient tensor network simulation of IBM’s kicked Ising experiment, arXiv: 2306.14887 (2023).

QC12
T. Begušić, J. Gray, and G. K.-L. Chan. Fast and converged classical simulations of evidence for the utility of quantum computing before fault tolerance, Sci. Adv. 10, eadk4321 (2024).





MPS1 
M. Fannes, B. Nachtergaele, and R. F. Werner, Finitely correlated states on quantum spin chains, Commun. Math. Phys. 144, 443 (1992).

MPS2
S. Ostlund and S. Rommer, Thermodynamic limit of density matrix renormalization, Phys. Rev. Lett. 75, 3537 (1995).





DMRG1
S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).

DMRG2
S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48, 10345 (1993). 

DMRG3
U. Schollwoeck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).





PEPS1 
F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many-body systems in two and higher dimensions, arXiv:cond-mat/0407066.

PEPS2 
F. Verstraete, J.I. Cirac, and V. Murg, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57, 143 (2008).





MERA1
L. Cincio, J. Dziarmaga, and M. M. Rams, Multiscale entanglement renormalization ansatz in two dimensions: quantum Ising model, Phys. Rev. Lett. 100, 240603 (2008).

MERA2 
G. Evenbly and G. Vidal, Entanglement renormalization in two spatial dimensions, Phys. Rev. Lett. 102, 180406 (2009).





PEPS3 
N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Computational Complexity of Projected Entangled Pair States, Phys. Rev. Lett. 98, 140506 (2007). https://doi.org/10.1103/PhysRevLett.98.140506DOI: 10.1103/PhysRevLett.98.140506. 

PEPS4 
T. Nishino and K. Okunishi, Corner Transfer Matrix Renormalization Group Method, JPSJ 65, 891 (1996). https://doi.org/10.1143/JPSJ.65.891DOI: 10.1143/JPSJ.65.891. 

PEPS5 
J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Classical simulation of infinite-size quantum lattice systems in two spatial dimensions, Phys. Rev. Lett. 101, 250602 (2008).

PEPS6 
R. Orús and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009).

PEPS7 
P. Corboz, T. M. Rice, and M. Troyer, Competing States in the t-J Model: Uniform d-Wave State versus Stripe State, Phys. Rev. Lett. 113, 046402 (2014).

PEPS8 
M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Faster methods for contracting infinite two-dimensional tensor networks, Phys. Rev. B 98, 235148 (2018).





BP0 
H. A. Bethe, Statistical theory of superlattices, Proc. R. Soc. Lond. A 150:552–575 (1935). http://doi.org/10.1098/rspa.1935.0122DOI: 10.1098/rspa.1935.0122. 

BP1 
J. Pearl. Reverend Bayes on Inference Engines: A Distributed Hierarchical Approach, Proceedings of the Second AAAI Conference on Artificial Intelligence. AAAI’82. Pittsburgh, Pennsylvania: AAAI Press, 1982, pp. 133–136. 

BP2 
J. Pearl. Fusion, propagation, and structuring in belief networks, Artificial Intelligence 29, 3 (1986), pp. 241–288. https://doi.org/10.1016/0004-3702(86)90072-XDOI: 10.1016/0004-3702(86)90072-X.   





BP3 
M. Mézard, G. Parisi, and R. Zecchina, Analytic and Algorithmic Solution of Random Satisfiability Problems, Science 297, 812 (2002). https://www.science.org/doi/pdf/10.1126/science.1073287DOI: 10.1126/science.1073287.

BP4 
J. S. Yedidia, W. T. Freeman, and Y. Weiss, Understanding Belief Propagation and Its Generalizations, Exploring Artificial Intelligence in the New Millennium. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 2003, pp. 239–269. ISBN: 1558608117.

BP5 
J. S. Yedidia, W. T. Freeman, and Y. Weiss, Constructing free-energy approximations and generalized belief propagation algorithms, IEEE Transactions on Information Theory 51.7 (2005), pp. 2282–2312. https://doi.org/10.1109/TIT.2005.850085DOI: 10.1109/TIT.2005.850085.





BP6 
M. Mézard, G. Parisi, and M. A. Virasoro, Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, Vol. 9 (World Scientific Publishing Company, 1987). https://doi.org/10.1142/0271DOI: 10.1142/0271.

BP7 
S. Yoon, A. V. Goltsev, S. N. Dorogovtsev, and J. F. F. Mendes, Belief-propagation algorithm and the Ising model on networks with arbitrary distributions of motifs, Phys. Rev. E 84, 041144 (2011).

BP8 
B. Karrer, M. E. J. Newman, and L. Zdeborová, Percolation on Sparse Networks, Phys. Rev. Lett. 113, 208702 (2014).





BP9 
M. S. Leifer and D. Poulin, Quantum Graphical Models and Belief Propagation, Annals of Physics 323.8 (2008), pp. 1899–1946. https://doi.org/10.1016/j.aop.2007.10.001DOI: 10.1016/j.aop.2007.10.001.
 
BP10 
D. Poulin and E. Bilgin, Belief propagation algorithm for computing correlation functions in finite-temperature quantum many-body systems on loopy graphs, Phys. Rev. A 77, 052318 (2008). https://doi.org/10.1103/PhysRevA.77.052318DOI: 10.1103/PhysRevA.77.052318.

BP11 
A. Wrigley, W. S. Lee, and N. Ye, Tensor Belief Propagation, Proceedings of the 34th International Conference on Machine Learning. Ed. by Doina Precup and Yee Whye Teh. Vol. 70. Proceedings of Machine Learning Research. PMLR, June 2017, pp. 3771–3779. 

BP12 
E. Robeva and A. Seigal, Duality of graphical models and tensor networks, Information and Inference: A Journal of the IMA 8.2 (June 2018), pp. 273–288. ISSN: 2049-8772. https://doi.org/10.1093/imaiai/iay009DOI: 10.1093/imaiai/iay009.

BP13 
R. Alkabetz and I. Arad, Tensor networks contraction and the belief propagation algorithm, Phys. Rev. Research 3, 023073 (2021). https://doi.org/10.1103/ PhysRevResearch.3.023073DOI: 10.1103/ PhysRevResearch.3.023073.

BP14 
S. Sahu and B. Swingle, Efficient tensor network simulation of quantum many-body physics on sparse graphs, arXiv: 2206.04701 (2022). 

BP15 
C. Guo, D. Poletti, and I. Arad, Block Belief Propagation Algorithm for 2D Tensor Networks, arXiv: 2301.05844 (2023). 

BP16 
N. Pancotti and J. Gray, One-step replica symmetry breaking in the language of tensor networks, arXiv: 2306.15004 (2023).

BP17 
Y. Wang, Y. E. Zhang, F. Pan, and P. Zhang, Tensor Network Message Passing, arXiv: 2305.01874 (2023). 

BP18 
J, Tindall, M. Fishman, Gauging tensor networks with belief propagation, SciPost Phys. 15, 222 (2023).





IBM1 
Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500–505 (2023).





ISO1 
M. P. Zaletal and F. Pollmann, Isometric Tensor Network States in Two Dimensions, Phys. Rev. Lett. 124, 037201 (2020).
https://doi.org/10.1103/PhysRevLett.124.037201DOI:10.1103/PhysRevLett.124.037201.





LOOP1 
M. Chertkov and V. Chernyak, Loop calculus in statistical physics and information science, Phys. Rev. E 73, 65102(R) (2006). 

LOOP2 
M. Chertkov and V. Chernyak, Loop series for discrete statistical models on graphs, J. Stat. Mech: Theory Exp. 6, P06009 (2006).

LOOP3
V. Gómez, J. M. Mooij, and H. J. Kappen, Truncating the Loop Series Expansion for Belief Propagation, Journal of Machine Learning Research 8 (2007) 1987-2016.





AKLT1
I Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Rigorous results on valence-bond ground states in antiferromagnets, Phys. Rev. Lett. 59, 799 (1987).

AKLT2
I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Valence Bond Ground States
in Isotropic Quantum Antiferromagnets, Comm. Math. Phys. 115, 477 (1988).





PESS1
Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, Tensor Renormalization of Quantum Many-Body Systems Using Projected Entangled Simplex States, Phys. Rev. X 4, 011025 (2014).





TRG1 
M. Levin and C. P. Nave, Tensor Renormalization Group Approach to Two-Dimensional Classical Lattice Models, Phys. Rev. Lett. 99, 120601 (2007).

TRG2 
G. Evenbly and G. Vidal, Tensor Network Renormalization, Phys. Rev. Lett. 115, 180405 (2015).

TRG3 
S. Yang, Z.-C. Gu, and X.-G. Wen, Loop Optimization for Tensor Network Renormalization, Phys. Rev. Lett. 118, 110504 (2017). https://doi.org/10.1103/PhysRevLett.118.110504DOI: 10.1103/PhysRevLett.118.110504.

TRG4 
G. Evenbly, Gauge fixing, canonical forms, and optimal truncations in tensor networks with closed loops, Phys. Rev. B 98, 085155 (2018). https://doi.org/10.1103/PhysRevB.98.085155DOI: 10.1103/PhysRevB.98.085155.

TRG5 
M. Hauru, C. Delcamp, and S. Mizera, Renormalization of tensor networks using graph-independent local truncations, Phys. Rev. B 97, 045111 (2018).





LBP1 
A. T. Ihler, J. W. Fischer, and A. S. Willsky, Loopy Belief Propagation: Convergence and Effects of Message Errors, J. Mach. Learn. Res. 6, 905 (2005).








§ EXPONENTIAL SUPPRESSION OF HIGH-DEGREE EXCITATIONS
 
The utility of the proposed loop expansions is built on the premise that the weight of excitations (i.e. on top of the BP vacuum) are exponentially diminishing in their degree; this property allows us to accurately approximate a tensor network contraction using only a few low-degree excitations. In this appendix we justify this premise.

Let 𝒯 be a finite tensor network containing at least a single closed loop. We assume that a BP fixed point has been found and that the tensors in network 𝒯 have been normalized such that their vacuum contributions are unity (or equivalently that the Bethe free energy is zero). Finally, we also assume that the BP fixed point accurately approximates the network contraction; that the difference between the full network contraction 𝒵(𝒯) and the BP approximation 𝒵_BP is small. We now consider isolating a single length-x loop of the network by absorbing the fixed point messages on all edges outside of the loop into their adjoining tensors, such that the reduced network is realized as the trace of a product of matrices as shown in the example of Fig. <ref>(a). For simplicity we assume that the matrices within the loop are copies of a single matrix A, which would only occur in practice if the original network possessed sufficient spatial symmetries, although the essence of our argument does not rely on this assumption. It follows that the network contraction can be written as 𝒵(𝒯) = {A^x}. The reduced network can be decomposed as the BP vs the non-BP component

    𝒵(𝒯) = {(A𝒫)^x} + {(A𝒫^C)^x},

with 𝒫 as the rank-1 projector onto message subspace and 𝒫^C its complement, see also Fig. <ref>(b). Note that 𝒵_BP = {(A𝒫)^x} = 1 due to the normalization condition imposed. Given that BP fixed points are known to correspond to stationary points of the Bethe free energy, and that we have assumed our known fixed point accurately approximates the network contraction, it follows that the projector 𝒫 onto message subspace should maximize 𝒵_BP. Thus 𝒫 must equal the outer product of the dominant left/right eigenvectors of matrix A (or, equivalently, the fixed point messages within the loop must by the left/right eigenvectors of A). Let {λ_0, λ_1, λ_2,…} be the eigenvalues of A ordered with descending magnitude, where the normalization condition implies that |λ_0|=1. It follows that the weight W(δ_x) of the loop correction term can be written as

    W(δ_x) = {(A𝒫^C)^x} = ∑_r=1^d( (λ_r)^x).

Let us now consider how Eq. <ref> would scale in the limit of large length-x (still assuming that the matrices within the loop are identical). Firstly, we should realize that the sub-leading eigenvalues should have magnitude strictly less than unity, i.e. |λ_1| < 1, otherwise the loop correction term would out-weight the BP term, violating the prior assumption that the BP fixed point accurately approximates the full network contraction 𝒵(𝒯). In the large-x limit Eq. <ref> will be dominated be the largest magnitude eigenvalues, so if λ_1 of A is n-fold degenerate it follows

    W(δ_x) ≈ n (λ_1)^x

which may be rewritten as 

    W(δ_x) ≈ e^-kx

with k=-log(λ_1) - (1/x)log(n). Notice that k tends to a positive constant with large-x, given that λ_1 < 1, such that the weight W(δ_x) tends to purely exponential decay in degree-x as was posited Eq. <ref> of the main text.

Although our argument assumed, for simplicity, that the transfer matrices with any loop were identical this condition is not strictly necessary; it is only required that the transfer matrices have a gap between their dominant and sub-leading eigenvalues in order to realize (approximate) exponential decay. Finding a general characterisation of the networks that are expected to have a gap in the transfer matrix spectrum associated to any closed loop is likely a difficult task (and certainly beyond the scope of this manuscript). However, given that the loop series expansion is only expected to be viable for networks where the loop corrections are small in comparison to the BP approximation, it follows that the gap condition must hold in this circumstance. Finally we remark that our argument is only valid for excitations that realize genus-1 loops not, for instance, higher genus loops like the δ_11 excitations shown in Fig. <ref>(d). None-the-less, it seems to be the case that the scaling of Eq. <ref> does hold for higher genus loop excitations.







§ LOOP EXCITATIONS IN THE THERMODYNAMIC LIMIT
 
In any finite tensor network it is possible, at least in principle, to enumerate over the different possible loop configurations, organize them by degree, and then evaluate the terms out to some desired degree. In contrast, in the thermodynamic limit, even counting the density of terms of some fixed degree may be an intractable combinatorial problem. To this end we now discuss two different approximations which can be applied to use the loop expansion to calculate the free energy density f, as defined in Eq. <ref>, of a homogeneous tensor network on a lattice in the thermodynamic limit (e.g. an iPEPS).
 


 §.§ Single-excitation approximation
 
The first method that we propose to estimate the contribution of loop excitations to the free energy density involves counting only the configurations that involve only a single excitation on top of the BP vacuum. Let 𝒯_N be a homogeneous tensor network on a square lattice of N sites which is assumed to evaluate to some scalar 𝒵_N. We focus on the task of evaluating the free energy density f_N defined

    f_N ≡ -log(𝒵_N)/N.

Let us assume that we have already found the BP fixed point of 𝒯_N, and that the network has been normalized such that the BP vacuum is unity, 𝒵_BP = 1 (or equivalently that the Bethe free energy is zero). 

We shall compute the contributions from configurations that have a single excitation, with all other tensor edges set in their BP ground state. For simplicity, here we only consider the corrections arising from the δ_4 excitations (occupying a single square on the lattice) although the method we describe is easily extendable to include arbitrary loop excitations. On the lattice of N sites there are N distinct locations where the δ_4 excitation could be placed, see Fig. <ref>(a), with each configuration having a weight of W(δ_4). It follows that the scalar 𝒵_N of the network can be approximated as 

    𝒵_N≈ 1 + N W(δ_4),

thus the free energy density can be further approximated as 

    f_N    ≈-log(1 + N W(δ_4))/N
       ≈ -W(δ_4).

To get to the second line of working in Eq. <ref> we have used log(1+ϵ) ≈ϵ which requires 1/N ≫ W(δ_4) to be valid. That the single-excitation approximation is only valid for small N makes sense; for larger N lattices the contributions from configurations with double-excitations (or higher) will dominate as their multiplicity scales at least as N^2. However, given that the result from Eq. <ref> is independent of N, it may still be used as an approximation to the free energy density f in the thermodynamic limit. Another way of interpreting this result is that we use the single-excitation approximation to calculate the free energy density on only a small patch of the otherwise infinite system. We can easily generalize Eq. <ref> to include arbitrary excitations δ_x,

    f ≈∑_{δ_x} -L(δ_x) W(δ_x),

where L(δ_x) is the number of locations per lattice site that the excitation δ_x can be placed.







 §.§ Multi-excitation approximation
 
The second method that we propose to estimate the contribution of loop excitations to the free energy density attempts to directly include contributions from configurations that have multiple excitations. For simplicity we again focus on the δ_4 excitations, noting that the derivation is easily extended to include arbitrary excitations. We now resolve the tensor network 𝒯_N as the BP contribution 𝒵_BP plus N distinct contributions arising from the placement of a δ_4 excitation on the lattice as shown Fig. <ref>(b). However, here in each contribution we freeze only edges connecting the excitation to their BP ground while leaving the other edges in the network free. It follows that each of these contributions has total weight W(δ_4) 𝒵̃_N-4, where 𝒵̃_N-4 represents the scalar from contracting the remaining tensor 𝒯̃_N-4 network on (N-4) sites. Thus the scalar 𝒵_N from the network contraction is approximated as

    𝒵_N≈ 1 + N W(δ_4) 𝒵̃_N-4.

Before continuing our derivation it is informative to reflect on the advantages and disadvantages of the proposed approximation. A key advantage of this proposal is that it implicitly includes multi-excitations; when resolving a contribution containing a single W(δ_4) excitation the remaining network 𝒯̃_N-4 could subsequently be resolved into a further sum of excitations. However, a significant disadvantage is that this approximation counts certain configurations multiple times, see Fig. <ref>(c-d) for examples. While it could be possible to use a more sophisticated strategy that subtracts out contributions from multiply-counted configurations, here we simply ignore them as in practice they seem to have minimal impact on accuracy.

In order to make headway with Eq. <ref> another approximation is needed; we assume that the remaining network 𝒯̃_N-4 has the same free energy density f_N as the original network 𝒯_N. In terms of the respective contraction scalars, 𝒵̃_N-4 and 𝒵_N, this implies 

    𝒵̃_N-4≈ e^4f_N𝒵_N.

Note that this approximation for 𝒵̃_N-4 ignores the boundary effects (i.e. resulting from fixing some of the network edges in their BP ground). Substitution of Eq. <ref> into Eq. <ref> gives

    𝒵_N≈ 1 + 𝒵_N(e^4f_N N W(δ_4)),

which can be rearranged to give

    𝒵_N≈(1 - e^4f_N N W(δ_4))^-1.

Finally, the free energy density f_N can be expressed as

    f_N    ≈log(1-e^4f_NN W(δ_4))/N
       ≈ -e^4f_N W(δ_4),

where the approximation log(1+ϵ) ≈ϵ is used to obtain the second line of working. This formula is easily generalized to include arbitrary excitations δ_x

    f ≈∑_{δ_x} -L(δ_x) W(δ_x) e^S(δ_x) f,

with L(δ_x) as the number of locations per lattice site that the excitation can be placed and S(δ_x) as the number of lattice sites that the excitation occupies. Notice that Eq. <ref> is identical to Eq. <ref> apart from the addition of the extra terms exp(S(δ_x) f), which we refer to as bulk suppression factors as they act to suppress excitations that occupy larger regions of the network. The appearance of these factors in the multi-excitation derivation makes intuitive sense: excitations that occupy larger regions leave less `room' for other excitations to occur and are penalised accordingly (whereas this consideration is irrelevant if only single excitations are allowed).

Finally, we remark that, in general, it will not be possible to solve Eq. <ref> directly for (the loop corrections to) the free energy density f as it appears on both sides of the equation. However, given that these corrections are assumed to be relatively weak (otherwise the series expansion itself would not be valid), we are able to use an iterative strategy to approximately solve for f. This strategy is as follows: first we set f=0 on the r.h.s of Eq. <ref> to compute a new value for f, before then feeding the new f back into the r.h.s and repeating until f is converged. In our benchmark calculations we observed that this strategy only required a few iterations in order to converge to high precision. 







 §.§ Comparison
 
We have presented two approaches, each based on different approximations, for using the loop expansion to estimate the free energy density on an infinite, homogeneous tensor network. A comparison of these approaches for the hexagonal lattice AKLT model is shown in Fig. <ref>. While both approximations show similar accuracy when the expansion includes only the δ_6 excitation, the accuracy of single-excitation approximation shows minimal improvement as higher-degree terms are included in the expansion. In contrast, the multi-excitation approximation appears to converge to the exact result as the expansion is taken to include higher-degree terms. Similar results were also observed for the AKLT models on square and kagome lattices, and for randomly initialized PEPS tensor networks. These results demonstrate that the bulk suppression factors exp(S(δ_x) f), which appear only in Eq. <ref> derived using the multi-excitation approximation, are vital to converge the loop expansions to high accuracy.



§ KAGOME AND SQUARE LATTICES
 
In this appendix we provide some of the details of the benchmark results provided for iPEPS on kagome and square lattices. An iPEPS for the kagome lattice is shown in Fig. <ref>(a). Although we could perform the series expansion on the (square norm) of this network directly, we instead found it preferable to perform some preliminary restructuring. The reasoning for this was to remove the triangular loops from the network which are otherwise problematic since (i) the weight of the excitations around the triangular loops are likely to be large (since the loops are short) and (ii) they generate a multitude of low-degree terms that would need to be included in the expansion. The restructuring is depicted in Fig. <ref>(b), where each (5-index) PEPS tensor is first decomposed as a product of three 3-index tensors, then the tensors surrounding each triangle of the kagome lattice are subsequently contracted into a single tensor (located at the center of each kagome lattice triangle). These manipulations result in a PEPS with a hexagonal lattice geometry but where the physical indices are located on tensors set on the edges of the hexagonal lattice as shown in Fig. <ref>(c). Upon forming the square norm ⟨ψ|ψ⟩ a purely hexagonal lattice tensor network is obtained, as shown in Fig. <ref>(d), which can be treated with the same expansion as described for the hexagonal lattice PEPS. Note that the manipulations performed in Fig. <ref>(b) may require an increase in PEPS bond dimension from m→ m^2, although in many cases the new bond dimension could be truncated without incurring significant error. This example of manipulating the kagome lattice into a hexagonal lattice highlights an interesting point: that the accuracy of the series expansion can depend on the particulars of the tensor network representation of a quantum state |ψ⟩, rather than solely on the properties of |ψ⟩ itself (e.g. entanglement and correlations). 

The benchmark results of Sect. <ref> also explore contraction of PEPS on the square lattice, both for the ground-state of the AKLT model as well as for randomly initialized networks. The square lattice results were computed similarly to the methods prescribed for the hexagonal lattice, both in terms of evaluating the series expansion and of evaluating the numerically exact results using boundary MPS. Several of the loop correction terms are shown in Fig. <ref>, where it can be seen that the number of distinct loop excitations of degree-x scales more quickly as a function of x than was seen with the hexagonal lattice. 









