Tensor renormalization group approach to 2D classical lattice models

Michael Levin1 and Cody P. Nave2

1
Department of Physics, Harvard University, Cambridge, Massachusetts 02138
2
Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
arXiv:cond-mat/0611687v2 [cond-mat.stat-mech] 29 Dec 2008

We describe a simple real space renormalization group technique for two dimensional classical
lattice models. The approach is similar in spirit to block spin methods, but at the same time it
is fundamentally based on the theory of quantum entanglement. In this sense, the technique can
be thought of as a classical analogue of DMRG. We demonstrate the method - which we call the
tensor renormalization group method - by computing the magnetization of the triangular lattice
Ising model.

Introduction: The density matrix renormalization n

group (DMRG) technique has proved extraordinarily p
j m
powerful in the analysis of one dimensional quantum sys- T T T
tems. [1, 2] Thus it is natural to try to develop an analo- i
q T k l
gous renormalization group method in higher dimensions.
Such a method could solve many currently intractable r
problems (such as the 2D Hubbard model).
Recent work has focused on generalizing DMRG to FIG. 1: A tensor network model on the honeycomb lattice.
higher dimensional quantum systems. [3] But it is also
natural to try to generalize to higher dimensional clas-
sical lattice models. While classical real space renor- by
malization group methods (such as block spin methods
[4]) have been around for many years, they have never e−S(i,j,k,...) = Tijk Tilm Tjnp Tkqr ... (1)
achieved the generality or precision of DMRG.
In this paper, we address this problem in the two di- where the product includes a tensor for each site of the
mensional case. We use ideas from quantum informa- lattice (Fig. 1). The partition function is the sum of all
tion theory to develop a numerical renormalization group the weights:
method that can effectively solve any two dimensional X X
classical lattice model. The technique - which we call Z= e−S(i,j,k,...) = Tijk Tilm Tjnp Tkqr ... (2)
the tensor renormalization group (TRG) method - has ijk... ijk...
no sign problem and works equally well for models with
complex weights. In other words, the partition function is obtained by tak-
Accurate numerical methods based on transfer matri- ing the product of all the tensors, contracting the pairs
ces [5, 6] have already been developed for 2D classical of indices on each bond.
systems. The advantage of the approach described here is The TRG method: The tensor renormalization group
that it is a fully isotropic coarse graining procedure, sim- method is a way to compute the partition function Z
ilar in spirit to block spin methods. It is thus naturally using a real space RG flow. We explain the method in
suited to investigating universal long distance physics. the case of the honeycomb lattice. Each coarse graining
Also, on a more theoretical level, the method reveals the iteration is made up of two separate steps. The first is
relationship between classical RG and quantum entangle- approximate, and the second is exact (Fig. 2). The first
ment. Finally, if only for its simplicity, we feel that the step is to find a tensor S such that
method is a useful numerical tool in two dimensions as X X
well as a natural candidate for higher dimensional gener- Slin Sjkn ≈ Tijm Tklm (3)
alizations. n m
Tensor network models: The tensor renormalization
group method applies to a set of classical lattice mod- (We explain how to find such a tensor S at the end of
els called ”tensor network models.” [7] Many well known this section). The relation (3) is a very useful property.
statistical mechanical models, such as the Ising model, Geometrically, it means that we can reconnect the lattice,
Potts model, and the six vertex model, can be written making the replacement
naturally as tensor network models. In fact, as we show i l i l
S
later, all classical lattice models with local interactions T T

can be written as tensor network models. j k j S k (4)

To describe a tensor network model on the honeycomb
lattice, one must specify a (cyclically symmetric) tensor wherever we like, without affecting the partition func-
Tijk with indices i, j, k running from 1 to D for some D. tion Z. Applying this to the bonds shown in Fig. 2a, we
The corresponding tensor network model has a degree change the honeycomb lattice to the new lattice shown
of freedom i = 1, ..., D on each bond of the honeycomb in Fig. 2b. The partition function is now given by con-
lattice. The weight for a configuration (i, j, k, ...) is given tracting the S tensors on the new lattice.
 2

T T
T T T T S
T’
S S S S S S
T T
T’ T’
{σa } {σb }
T T T T S S
S S
T T S S S S S T’ T’ n n
S
T T T T
T’
T T S
T T T T
S
T’ {σcn }
S S S S S S
T T
T T T T S S
T’ T’ Ψ({σ a },{σb },{σc })
n n n
(a) (b) (c) R
(a) (b)
FIG. 2: A TRG transformation on the honeycomb lattice.
FIG. 4: After (a) dividing R into triangles, (b) the parti-
S S S S T’ T’
tion function of each triangle can be written as a function
T T T T Ψ({σan }, {σbn }, {σcn }) of the boundary spins.
S S S S S S S S T’ T’
T T T T S S
S S S S S
S
T’ T’
T T T T S S S S
S S S S T’ T’
T T T T S S
S S
S S
S S T’ T’
However an approximate solution can be obtained as
T T T T
S S
S S
S S
S S T’ T’
follows. The idea is to choose the matrix S that mini-
T T T T S S S S T’
mizes the error, |M − S · S T |2 ; the optimal S can then
T’
be found using the singular value decomposition P of M .
∗
In more detail: first one writes Mli,jk = n sn Uli,n Vjk,n
FIG. 3: A TRG transformation on the square lattice.
(here, sn are the singular values and U, V are unitary
matrices). Second, one truncates the matrices Uli,n , Vjk,n
keeping only those columns corresponding to the largest
The second step is now clear. We group together
D singular values. The result are D2 × D matrices
triplets of neighboring points replacing them by a sin- √
Ũli,n , Ṽjk,n . Finally, one sets Slin = sn Ũli,n . This gives
gle lattice point with a coarse-grained tensor T ′ :
the required factorization - provided that we adjust the
i S S k i k phase ambiguity Uli,n → Uli,n eiφn , Vjk,n → Vjk,n e−iφn
S
T’ appropriately. In practice, it is often more convenient
j j (5) to ignore the phase adjustment issue and set Slin A
=
√ B √ ∗
sn Ũli,n , Sjkn = sn Ṽjk,n . The result is a factorization
Here the tensor T ′ is given by contracting over the three P A B P A B
bonds of the triangle: n Slin Sjkn ≈ m Tijm Tklm where S and S differ by
some phase factors. The TRG procedure can be applied
X as before - the only difference being that we have to keep
′
Tijk = Skpq Sjqr Sirp (6) track of two different tensors T A , T B for the A and B
pqr
sublattices.
Making this replacement everywhere gives a new We will show in the next section that the error for this
(coarser) honeycomb lattice (see Fig. 2c). This com- optimal decomposition is independent of the number of
pletes the coarse graining transformation. The end result iterations and can be made arbitrarily small by increas-
is that the number of points in the lattice has decreased ing D. Indeed, the error vanishes as ǫ ∼ exp(−const ·
by a factor of 3 and T has been replaced by T ′ . (log D)2 ) - the same scaling behavior as the truncation
Iterating this procedure, one can compute the parti- error in DMRG. [2]
tion function of an arbitrarily large finite lattice. Ther- Physical picture: In this section we explain the physics
modynamic observables and correlation functions can be behind the TRG method, and give a physical inter-
obtained by taking numerical derivatives of F = − log Z, pretation of the fixed point tensor T ∗ . We begin by
or by evaluating the free energy of more general models showing how arbitrary classical models in two dimen-
where the tensors Tijk vary from site to site. An addi- sions can be thought of as tensor models on the hon-
tional feature is that the tensors T converge to a fixed eycomb lattice. For concreteness we frame our discus-
point tensor T ∗ - whose physical significance we explain sion around
P the case
P of the square lattice Ising model,
in the next section. The method is not limited to the hon- Z = {σi } exp(K hiji σi σj ). Consider the partition
eycomb lattice and can easily be implemented on other function ZR for some finite region R in the plane. One
lattices (see Fig. 3). way to compute ZR is to triangulate R, dividing it into
To complete the discussion, we address the issue of triangles of size L much larger than the lattice spacing l
finding aP tensor S which satisfies (3). The first step is to (Fig. 4a).
2 2
think
P of m Tijm Tklm as a D × D matrix M : Mli,jk = Consider one of the triangles. Imagine summing over
m Tijm Tklm . It is also useful to think of the tensor Slin all the lattice degrees of freedom within the triangle. The
as a D2 × D matrix Sli,n . Then the problem of satisfying result will yield some number Ψ = Ψ({σin }) that depends
(3) is the problem of finding a D2 × D matrix S such on the values of {σin } at the boundary of the triangle
that M = S · S T . In general, this factorization cannot be (Fig. 4b). It is convenient to separate out these boundary
done exactly since the left hand side typically has rank degrees of freedom into three groups {σan }, {σbn }, {σcn }
D2 , while the right hand side has rank at most D. corresponding to the three sides a, b, c of the triangle.
 3

(a) (b) (c) M

D=4
D = 10
FIG. 5: The tensor renormalization group transformation D = 16
D = 22
can be viewed as a two step change in the triangulation of R. D = 28
D = 34
Exact

Denoting {σan } schematically by α, and similarly for 1e-05 0.0001 0.001 0.01 0.1 1
{σbn }, {σcn }, we can think of Ψ as a three index tensor (αc- α)/αc
Ψ = Ψαβγ .
To obtain the partition function for the region R, FIG. 6: The magnetization of the triangular lattice Ising
we simply need to glue together all the triangles and model as a function of α = e−2βJ , obtained using the TRG
then method.
P sum over the spins at their boundaries: ZR =
αβγδǫ... Ψαβγ Ψαδǫ .... This is nothing but the partition
function of a tensor model (2). (Actually, a more care- Ψ({σan }, {σbn }, {σcn }) of the spins on the three sides of
ful analysis shows that the resulting tensor model has the triangle. In fact, this function should be thought
different tensors ΨA , ΨB for the A, B sublattices). of as a wave function for a one dimensional quantum
It is useful to think about the tensor renormalization spin system with spins living on the boundary of a tri-
group transformation in this context. Recall P that the angle. This interpretation comes from thinking of the
first
P step is to find tensors S satisfying n S lin Sjkn ≈ original 2 dimensional classical model (e.g. the Ising
m Tijm Tklm . Thinking in terms of the triangles, model) as a (1 + 1) dimensional quantum model. We
the right hand side is simply the partition function of think of the direction parallel to the boundary of the
the rhombus obtained by gluing two triangles together. triangle as space, and the (radial) direction perpendic-
Thus, a solution S can be constructed by setting S equal ular to the boundary as time. We imagine constructing
to the partition function of one of the obtuse triangles a one dimensional transfer matrix/quantum Hamiltonian
obtained by dividing the rhombus in the other direction: H living on the boundary of the triangle, whose (radial)
time evolution generates the two dimensional classical
S model in question. Then in this picture, the function
TT
S Ψ({σan }, {σbn }, {σcn }) is the result of evolving H for a
(7)
long time (the triangle is large). Hence Ψ is simply the
Thus, the first step of the renormalization transforma- ground state of H - up to exponentially small corrections.
tion process simply changes the triangulation as shown If we assume that the original classical model (e.g. the
in Fig. 5b. The second step also has a simple interpreta- Ising model) is not critical, then Ψ is the ground state
tion. Examining the definition of T ′ , it is not hard to see of a gapped Hamiltonian. Gapped ground states in one
that T ′ is simply the partition function of a large equi- dimension have an important property: they are only
lateral triangle, obtained by gluing together three obtuse weakly entangled. More specifically, it is known that the
triangles: density matrix ρ(x, L) of a region of size x converges to
a fixed density matrix ρ∞ as x, L → ∞. Moreover, the
size of the mth eigenvalue λm of ρ∞ falls off rapidly with
S S
T’ increasing m: λm ∼ exp(−const · log(m)2 ). [8]
S (8)
It is this property which guarantees the accuracy of the
Thus, the second step simply glues together triplets of TRG method. Indeed, as a consequence of this property,
obtuse triangles to form larger equilateral triangles as one can factor Ψ into a product of three spin states on
shown in Fig. 5c. In this way, the TRG method builds the three sides of the triangle,
up larger and larger triangles.
X
D
However, there is one subtlety. As we build up larger Ψ(σa , σb , σc ) ≈ Tijk ΨiA (σa )ΨjB (σb )ΨkC (σc ) (9)
and larger triangles, the corresponding tensors will have ijk=1
indices
√
with larger and larger ranges, increasing from 2L/l
3L/l
to 2 to 23L/l and so on. Yet, the TRG method with high accuracy. Such a factorization can be obtained
insists on approximating these tensors by a tensor with a by choosing ΨiA , ΨiB , ΨiC to be the ith largest eigenstates
fixed range D. How can this approximation possibly be of the density matrices of sides A, B, C, and letting Tijk
accurate? be the matrix elements of the tensor Ψ(σa , σb , σc ) be-
To answer this question, we must explain the physical tween these states. By the discussion above, the error ǫ
meaning of Ψαβγ . Let Ψ be a tensor obtained from the of this representation is (1) independent of the size L of
partition function of a very large triangle. Writing out the triangle, and (2) decreases rapidly with increasing D
the labels α, β, γ explicitly, we can write Ψ as a function (ǫ ∼ exp(−const · (log D)2 )).
 4

Thus, even though the exact tensor Ψαβγ has an expo- α = e−2βJ gives rise to the correct Boltzmann weight.
nentially large range, one can make a change of basis so
that in that basis, Ψ can be accurately approximated by a Applying the TRG method to the above tensor, we
tensor Tijk whose indices have a fixed range D. The TRG compute the free energy per unit spin, F = − N1β log(Z)
method can be thought of as a numerical technique for
in the thermodynamic limit N → ∞. The magnetization
(approximately) constructing this tensor Tijk for larger
M can be obtained by taking numerical derivatives of F
and larger triangles. The fixed point T ∗ is the value of
(though we need to use a more complicated tensor T to
this tensor in the limit of an infinitely large triangle.
represent an Ising model with an external magnetic field
The above analysis was based on the assumption that
H). Increasing D, the computation rapidly converges
the classical model was not critical. If instead the clas-
to the exact result [10] except in an increasingly
√ narrow
sical model is critical, the associated quantum states are
interval around the critical point αc = 1/ 3 (Fig. 6).
gapless ground states. Gapless ground states are more
entangled then their gapped counterparts. The entan-
glement entropy S = −Tr(ρ log ρ) of a region of size x in This interval becomes so narrow that one can study
a system of size L grows logarithmically with the region the critical point itself. For example, the magnetization
size, S ∼ log x. [9] This means that the factorization (9) curve for D = 34 predicts an αc within 10−4 of the exact
will always break down when the triangle is sufficiently result. One can even estimate the critical exponent β
large. Thus, in principle the TRG method - like DMRG from the scaling behavior of the magnetization. We find
[2]- breaks down at criticality. β = 0.12, not far from the exact value β = 1/8. How-
A simple example: In this section we demonstrate the ever, as explained earlier, the TRG method (like DMRG)
method with a simple P example:Pthe triangular lattice is best suited to studying systems off criticality. An inter-
Ising model: Z = {σ} exp(βJ hiji σi σj ). Note that esting question for further research is whether the TRG
the Ising model partition function can be written as a method can be modified so that it’s (almost exponential)
sum over domain wall configurations where the domain accuracy is uniform, both away from and near critical
walls live on the bonds of the honeycomb lattice. This points.
domain wall model can be easily realized by a tensor net-
work with D = 2. We think of the state i = 1 as denoting The authors would like to thank Michael Freedman and
”no domain wall” and i = 2 as denoting ”domain wall.” Frank Verstraete for motivating this problem and Nihat
Then the tensor with nonzero components Berker for useful discussions. This work was supported
by the Harvard Society of Fellows (M.L.), and by NSF
T 111 = 1, T 122 = T 212 = T 221 = α, (10) grant DMR-0517222 (C.P.N.).

[1] S. R. White, Phys. Rev. Lett. 69, 2863 (1992). (2006).

[2] U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005). [8] K. Okunishi, Y. Hieida, and Y. Akutsu, Phys. Rev. E 59,
[3] F. Verstraete and J. Cirac, cond-mat/0407066 (2004). R6227 (1999).
[4] K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975). [9] G. Vidal, J. Latorre, E. Rico, and A. Kitaev, Phys. Rev.
[5] T. Nishino, J. Phys. Soc. Jpn. 64, 3598 (1995). Lett. 90, 227902 (2003).
[6] V. Murg, F. Verstraete, and J. Cirac, Phys. Rev. Lett. [10] R. Houtappel, Physica 16, 425 (1950).
95, 057206 (2005).
[7] Y. Shi, L. Duan, and G. Vidal, Phys. Rev. A 74, 022320