TensorNetwork on TensorFlow: Entanglement Renormalization for quantum critical lattice models

TL;DR

This paper implements a tensor network optimization algorithm using TensorFlow to study quantum critical models, achieving significant speed-ups on GPUs for ground state approximations and conformal data extraction.

Contribution

It introduces a GPU-accelerated tensor network optimization method for MERA using TensorFlow, enabling efficient analysis of quantum critical systems.

Findings

GPU implementation speeds up optimization by up to 200 times

Successfully approximates ground states of critical quantum models

Extracts conformal data from optimized tensor networks

Abstract

We use TensorNetwork [C. Roberts et al., arXiv: 1905.01330], a recently developed API for performing tensor network contractions using accelerated backends such as TensorFlow, to implement an optimization algorithm for the Multi-scale Entanglement Renormalization Ansatz (MERA). We use the MERA to approximate the ground state wave function of the infinite, one-dimensional transverse field Ising model at criticality, and extract conformal data from the optimized ansatz. Comparing run times of the optimization on CPUs vs. GPU, we report a very significant speed-up, up to a factor of 200, of the optimization algorithm when run on a GPU.

Equations

$$H=\sum _{n\in \text{mathbbm}Z}-X_n{X}_{n+1}-h{Z}_n\equiv \sum _{n\in \text{mathbbm}Z}{h}_{n,n+1}$$

$$\sum _{\alpha \beta }[w{]}_{\alpha \beta }^{\gamma }[\overline{w}{]}_{\alpha \beta }^{{\gamma }^{\prime }}={\delta }_{\gamma {\gamma }^{\prime }}$$

$$\sum _{\alpha \beta }[u{]}_{\alpha \beta }^{\mu \nu }[\overline{u}{]}_{\alpha \beta }^{{\mu }^{\prime }{\nu }^{\prime }}={\delta }_{\mu {\mu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}$$

$$h_{n,n+1,n+2}^{\tau =1}$$

$$H$$

$$\mathcal{O}\to 2^{-{\Delta }_{\mathcal{O}}}\mathcal{O}$$

Full text

I Introduction

Tensor networks are efficient and very powerful parametrizations of restricted classes of vectors in high-dimensional vector spaces. They were originally developed to simulate many-body quantum systems in condensed matter physics. Tensor networks like the Matrix Product State (MPS) Fannes et al. (1992); White (1992); Bridgeman and Chubb (2017), the Multi-Scale Entanglement Renormalization Ansatz (MERA) Vidal (2007, 2006) or Projected Entangled Pair States (PEPS) Verstraete et al. (2006); Verstraete and Cirac (2004) are nowadays routinely used to approximate ground states of quantum many-body Hamiltonians in one (1d) and two (2d) spatial dimensions or to simulate real-time dynamics of many-body quantum states. Over the past decade it has become evident that tensor networks are also useful well outside their original scope, and they have for example found applications in quantum-chemistry White and Martin (1999); Chan et al. (2007); Szalay et al. (2015); Krumnow et al. (2016); White and Stoudenmire (2019), real material calculations Bauernfeind et al. (2017); Ganahl et al. (2015, 2014), quantum field theory Verstraete and Cirac (2010); Haegeman et al. (2013); Ganahl et al. (2017), machine learning Stoudenmire and Schwab (2016); Stoudenmire (2018); Evenbly (2019) and even quantum gravity Swingle (2012); Bény (2013); Czech et al. (2016); Bao et al. (2017); Milsted and Vidal (2018).

TensorFlow Abadi et al. (2015) is a free, open source software library for data flow and differentiable programming, developed by the Google Brain team, that can be used for a range of tasks including machine learning applications such as neural networks. Recently, the open source library TensorNetwork Roberts et al. (2019, 2015) has been released to facilitate the development and implementation of tensor network algorithms in TensorFlow.

In previous work we have provided benchmark results for the use of TensorNetwork for tensor network computations in condensed matter physics Milsted et al. (2019) and machine learning Efthymiou et al. (2019). In Milsted et al. (2019) we have used TensorNetwork to implement and benchmark an optimization algorithm for so-called Tree Tensor Networks (TTNs) to approximate the ground state of the Ising model on a torus. Using TensorNetwork’s TensorFlow backend to run the optimization on accelerated hardware, we obtained speed-ups of a factor of 100 over optimization on CPUs. Similar speed-ups are reported in Efthymiou et al. (2019) for applications of TensorNetwork to classification tasks in the area of machine learning.

In this manuscript we use TensorNetwork to implement an optimization algorithm for MERA. The purpose of this paper is to (1) provide sample code which shows how to use TensorNetwork to implement a MERA optimization. This is of interest to tensor network practitioners who want to get started quickly with TensorNetwork, as well as for newcomers who want to understand basic concepts of MERA optimizations; (2) to benchmark runtimes for MERA on CPU and GPU, and analyze in detail how the computational cost is distributed among different computational tasks.

To motivate the MERA, let us briefly review the basics of binary TTNs, shown in Fig.(1) (a). In a TTN all tensors (shown as blue triangles) are structurally identical (apart from the top-most tensor, shown in red). The TTN shown in Fig.(1) (a) can be contracted and optimized efficiently. Similar to a TTN, the MERA (shown in Fig.(1) (b)) is a tensor network designed to approximate ground states of many-body Hamiltonians. The MERA is a powerful generalization of a TTN which introduces a second set of structurally different tensors, called disentanglers (shown as green squares in Fig.(1) (b)) into the network. Abundant numerical evidence has long established that the MERA is well suited to approximating ground states of critical quantum systems Evenbly and Vidal (2009), more so than for example the simple TTN in Fig.(1) (a). In this paper we review a standard MERA energy minimization algorithm from Refs. Evenbly and Vidal (2009); Pfeifer et al. (2009); Evenbly and Vidal (2011), which results in a (scale invariant) MERA representation of an infinite, critical quantum spin chain. While MERA is a more powerful tensor network than a TTN, the addition of disentanglers significantly increases the computational complexity for fixed bond dimension $\chi$ (the leading cost changes from $\chi^{4}$ to $\chi^{9}$).

It is expected that, similar to TTNs Milsted et al. (2019), MERA algorithms can be considerably sped up by running them on accelerated hardware like GPUs or TPUs. Using TensorNetwork with TensorFlow as backend, we confirm this expectation and show that running a MERA optimization algorithm on state-of-the-art accelerated hardware can yield speed-ups of up to a factor of 200 as compared to running it on a CPU. The code for all calculations can be downloaded from Roberts et al. (2015).

II Review of the Multi-scale Entanglement Renormalization Ansatz

The purpose of this section is to introduce the reader to notation and give a high-level introduction to MERA. For a more comprehensive review, we refer the reader to Evenbly and Vidal (2009). The MERA is a class of variational tensor network wave functions which can be used to efficiently approximate ground states of critical quantum many-body Hamiltonians. In the following we will use the scale invariant MERA Pfeifer et al. (2009) (see Sec.(II.1) and Fig.(2)), to approximate the ground state of critical quantum lattice models in the thermodynamic limit on a one dimensional (1d) lattice $\mathcal{L}$. As a benchmark example, we will consider the quantum Ising model with Hamiltonian

[TABLE]

with X=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right) and Z=\left(\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right) Pauli spin operators, i.e. each site on the lattice $\mathcal{L}$ hosts a spin $1/2$ degree of freedom. The Hamiltonian has a phase transition at a magnetic field $h=1$ from a paramagnetic ($h>1$) to a ferromagnetic ($h<1$) phase. At the critical point $h=1$ the spectrum of $H$ becomes gapless and ground state correlations decay algebraically (that is, as a power law $x^{-2\Delta}$ with distance $x$). One of the key advantages of MERA over other variational classes of wave functions is its ability to exactly capture algebraic correlations at arbitrary distances. This makes it an ideal ansatz for ground states of critical quantum systems.

II.1 Scale invariant MERA

The scale invariant MERA is a tensor network representation of a many-body wave function $\ket{\Psi}$ on an infinite lattice $\mathcal{L}$. Fig.(2) shows an example of a scale invariant MERA. The physical degrees of freedom of this lattice are denoted by $i_{n}$, where $n\in\mathbbm{Z}$ denotes the lattice site, and $i_{n}\in\{1,2\}$ for a two level quantum spin degree of freedom, such as in the transverse field Ising model. In this paper we will focus on the so-called binary MERA. It consists of two types of tensors (see Fig.(1) (b) and (c)) called isometries $w$ and disentanglers $u$. Isometries $[w]_{\alpha\beta}^{\gamma}$ are rank-3 tensors of dimension $\chi\times\chi\times\chi$ which obey the isometric constraint

[TABLE]

while disentanglers $[u]_{\alpha\beta}^{\mu\nu}$ are rank-4 tensors of dimension $\chi\times\chi\times\chi\times\chi$ obeying

[TABLE]

see Fig.(1) (d). Square brackets are used to denote tensor elements, overlines denote complex conjugation, and $\delta_{\mu\nu}$ is the Kronecker delta. To simplify explanations, we assume that all isometries $w$ and disentanglers $u$ have identical dimensions $\chi$. A MERA is organized into layers $\tau=1,2,3,\dots\;$. In Fig.(2) we show the first four layers of a scale invariant MERA. Each layer corresponds to a different scale Evenbly and Vidal (2009) of the network, and contains a row of isometries $w_{\tau n}$ and a row of disentanglers $u_{\tau n}$, where $\tau$ and $n$ label layer (or scale) and position, respectively. In the following we impose translational invariance by choosing $w_{\tau n}$ and $u_{\tau n}$ to be independent of $n$ within each layer $\tau$, and will henceforth omit the subscript $n$. The physical degrees of freedom $i_{n}$ reside at the bottom of the network. The scale invariant MERA consists of an infinite number of layers $\tau$. The layers are divided into transitional layers with $\tau<\tau^{*}$ (where $\tau^{*}$ is a hyper parameter of the network) and scale invariant layers $\tau\geq\tau^{*}$. For instance, in Fig.(2) we have $\tau^{*}=3$, indicating that there are two transitional layers of tensors for $\tau=1,2$ before the scale invariant layers for $\tau\geq 3$. For the scale invariant layers, all tensors $u_{\tau},w_{\tau}$ are identical and we will occasionally denote them by $w$ and $u$ without subscript. Tensors in the transitional layers are different for each layer. That is, if $\tau^{*}=3$ then the entire MERA is specified by the six tensors $u_{1},w_{1},u_{2},w_{2},u_{3},w_{3}$.

II.2 Causal cone, local observables, scaling operators

In the following we explain how to efficiently compute observables for a given MERA. As a concrete example, we will focus on the expectation value of the Hamiltonian Eq.(1). In a binary MERA, Hamiltonians appear naturally as sums of three body terms. It is therefore convenient to regroup the Hamiltonian Eq.(1) into a sum of three body terms:

[TABLE]

We refer to $h_{n,n+1,n+2}^{\tau=1}$ as the Hamiltonian density in the first layer. To calculate the expectation value $\bra{\Psi}h_{n,n+1,n+2}^{\tau=1}\ket{\Psi}$ of a local Hamiltonian term, we need to contract an infinite network of isometries and disentanglers. Fig.(3) shows a diagrammatic representation of $\bra{\Psi}h_{n,n+1,n+2}^{\tau=1}\ket{\Psi}$. One of the key features of the MERA is the presence of a causal cone Evenbly and Vidal (2009, 2011); Evenbly (2015) which allows for an efficient contraction of the network. The causal cone for a binary MERA is shown in Fig.(3) as highlighted region. Using the isometric and unitary conditions Eq.(2) and Eq.(3) (see also Fig.(1) (c) and (d)), all contractions outside the causal cone can be carried out trivially, and the only contractions that need to be performed explicitly are over tensors inside the causal cone. This is illustrated in Fig.(3). The contraction of the network shown in the top-right corner of Fig.(3) is done in several steps: the first step is the contraction of the infinite number of tensors inside the scale invariant layers of the causal cone. This is achieved by calculating the dominant eigenvector $\rho_{\tau^{*}}$ of the descending super operator $\mathcal{D}$ (defined in Fig.(4) (b)) by iterative methods Evenbly and Vidal (2009). In the second step, $\rho_{\tau^{*}}$ is then descended through the transitional layers of the MERA down to the bottom of the network using the descending super-operator of the transitional layers. Finally, at the bottom of the network, $\rho_{\tau=1}$ is contracted with the operator $h_{n,n+1,n+2}^{\tau=1}$ to give its expected value. For the scale invariant binary MERA, all contractions can be carried out at a leading computational cost of $\chi^{9}$.

Two hallmarks of quantum criticality are the phenomena of scale invariance and universality. Scale invariance refers to the fact that at a critical point, a physical system looks similar when observed at different length scales. Universality refers to the fact that microscopically different systems may have identical properties at their respective critical points. More generally, critical systems can be divided into universality classes. The universal physical properties of systems within one such class can be described by so-called scaling operators. Scaling operators are operators which transform covariantly under changes of scale. For example, under a change of scale by a factor of two, a scaling operator $\mathcal{O}$ transforms as Pfeifer et al. (2009)

[TABLE]

where $\Delta_{\mathcal{O}}$ is called the scaling dimension of operator $\mathcal{O}$. In a so-called conformal field theories Di Francesco et al. (1997), scaling operators are organized into distinct conformal towers. Operators at the bottom of a tower are called primaries, and those higher up are called descendants. There is one tower for each primary operator. The number of primary operators is a property of the physical system. For the Ising model at criticality, the low energy spectrum of Eq.(1) is described by the so-called Ising conformal field theory (Ising CFT) Di Francesco et al. (1997) with central charge $c=\frac{1}{2}$ and three primary operators $\mathbbm{1},\sigma$ and $\epsilon$, with scaling dimensions $\Delta_{\mathbbm{1}}=0,\Delta_{\sigma}=\frac{1}{8}$ and $\Delta_{\epsilon}=1$, respectively.

One of the highlights of MERA is its realization of a discrete scale transformation by a factor of two for the binary MERA (we refer the reader to Evenbly and Vidal (2009); Pfeifer et al. (2009); Evenbly and Vidal (2011); Evenbly (2015) for details). From an optimized MERA wave function, lattice versions of scaling operators and their corresponding scaling dimensions can be approximately obtained from diagonalizing the scale invariant ascending super-operator $\mathcal{A}$ of Fig.(4) (a).

II.3 Optimization

We use a standard variational optimization algorithm to optimize the scale invariant MERA to approximate the ground state of Hamiltonian $H$ in Eq.(1). The algorithm iteratively updates isometries and disentanglers to lower the expectation value of the energy, which is calculated as outlined in the previous section. We call one round of updates a sweep. One sweep consists of several steps, which we will describe in the following. In the first step, the isometries and disentanglers in the scale invariant layers $\tau\geq\tau^{*}$ are updated. In the second step, isometries and disentanglers in the transitional layers $\tau<\tau^{*}$ are updated layer by layer (we refer the reader to Evenbly and Vidal (2009, 2011); Evenbly (2015) for a detailed discussion of optimization algorithms). To update isometries $w_{\tau}$ and disentanglers $u_{\tau}$ in a layer $\tau$, we need to calculate the environments of the corresponding tensors. The environments of $w_{\tau}$ and $u_{\tau}$ are rank-3 and rank-4 tensors, respectively. Exemplary contributions to the environments of the isometry $w_{\tau}$ and the disentangler $u_{\tau}$ are shown in Fig.(4) (c) and (d). To calculate the environments in layer $\tau$, we use the descending super-operator $\mathcal{D}$ to descend $\rho_{\tau^{*}}$ to layer $\tau+1$, and the ascending super-operator $\mathcal{A}$ to ascend the local Hamiltonian density $h_{n,n+1,n+2}^{\tau=1}$ up to layer $\tau$ . The full environments consist of sums of several contributions of similar form (see Fig.(A1)). Finally, the updates to $u_{\tau}$ and $w_{\tau}$ are obtained from an SVD of the environments. This is shown graphically in Fig.(5). The sweeps are repeated until sufficient convergence is reached.

III Benchmark results

In the following we present benchmark results for the MERA approximation of the ground state of the critical transverse field Ising model at magnetic field $h=1$, Eq.(1). We have performed calculations for bond dimensions $\chi=4,6,8,10,12,14$ and $16$.

III.1 Ground state energy, scaling dimensions

In Fig.(6) we present our benchmark results for the ground state energy for different bond dimensions $\chi$ of the MERA. The exact ground state energy density is given by $e_{exact}=-\frac{4}{\pi}$. The relative error in the ground state energy decays quickly with increasing bond dimension. At our largest bond dimension $\chi=16$, the error is $\approx 5\times 10^{-10}$.

In Fig.(7) we show the lowest twelve scaling dimensions $\Delta$ obtained from the eigenvalues $\lambda=2^{-\Delta}$ of the ascending super-operator $\mathcal{A}$ (blue dots), together with the exact scaling dimensions of the Ising CFT (red lines), for a bond dimension of $\chi=16$. We observe a very good agreement with the theoretical results. Larger scaling dimensions are significantly less accurate.

III.2 Runtimes

One of the key advantages of the TensorNetwork API is the use of TensorFlow as a backend to perform tensor contractions. Thus, it can be easily run on accelerated hardware like GPUs or TPUs, with the benefit of considerable reduction of runtimes. In the following we present benchmark results comparing optimization runtimes on CPU and GPU. We ran our simulations on TensorFlow v1.13.1 built with the Intel math kernel library (MKL) on Google’s cloud compute engine. CPU benchmarks were run on an Intel Skylake architecture with 1, 16, 32, and 64 cores. The GPU benchmarks were run on an NVIDIA Tesla V100 graphics processor 111SVDs were computed on the host CPU, which was seen to be faster than performing the SVD on GPU. Since the SVDs contribute only a small fraction to the total runtime, the overall runtimes were essentially unaffected by this choice..

The most expensive operations of the optimization algorithm described above scale asymptotically as $O(\chi^{9})$. In Fig.(8) we compare average runtimes of one step of the optimization (averaged over twentyx iteration steps) run on CPUs and GPU. When run on one CPU (single thread), the asymptotic scaling is seen to be $\chi^{9}$, consistent with the theoretical expectation. For single-threaded operation and large bond dimensions, we observe a $200$ fold speed-up of the optimization when run on a GPU, as compared to CPU. For multi-threaded optimization on multiple CPUs, we observe a convergence of the speed-up with the number of CPUs at 32 threads. Compared to 32 CPUs, the GPU is still faster by a factor of $6$.

In Fig.(9) we show the runtimes of the individual steps of a MERA optimization. The upper panel shows runtimes on a GPU running with a single thread, the lower panel results for CPU running with a single thread. As expected, in either case the dominant cost is given by the computation of the steady-state reduced density matrix $\rho^{*}$ and the computation of the environments. Running on GPU gives speed-ups of these operations by a factor of  200. Note that the SVDs in both cases were performed on CPU. At the largest bond dimension $\chi=16$ the cost for calculating $\rho^{*}$ is still about 40 times larger than the SVD of the environment of $u$. This is in contrast to the case of the TTN Milsted et al. (2019), where the factor between the most expensive tensor contractions and the SVD is on the order of 3. For this reason, the MERA optimization algorithm exhibits a greater speed-up as compared to a TTN optimization when run on a single thread.

IV Conclusion

We have used TensorNetwork, an API for tensor contractions with a TensorFlow backend, to implement the scale invariant MERA optimization algorithm of Ref. Pfeifer et al. (2009); Evenbly and Vidal (2011). As a benchmark, we have used the algorithm to approximate the ground state of the critical transverse field Ising model in the thermodynamic limit, using MERA with bond dimensions of up to $\chi=16$. From the optimized MERA, we calculated the lowest twelve scaling dimensions of the Ising model, which are found to be in excellent agreement with their theoretically predicted values. When run on a GPU, we observe a $200$ fold speed-up as compared to 1 CPU core, and a speed-up of $6$ when compared to 32 CPUs.

Acknowledgements M. Ganahl, A. Milsted and G. Vidal thank X for their hospitality. X is formerly known as Google[x] and is part of the Alphabet family of companies, which includes Google, Verily, Waymo, and others (www.x.company). G. Vidal is a CIFAR fellow in the Quantum Information Science Program. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science.

V Appendix

V.1 Environment contributions

For completeness, in Fig.(A1) we show all contributions to the environments of isometry $w$ and disentangler $u$.