Differentiable Programming Tensor Networks

TL;DR

This paper introduces a differentiable programming framework for tensor network algorithms, enabling automatic differentiation and efficient gradient computations, which simplifies development and enhances optimization in tensor network applications.

Contribution

It presents a novel approach to differentiable tensor network programming, including techniques for stable differentiation and backpropagation through complex tensor operations.

Findings

Computed the specific heat of the Ising model via second derivatives.

Achieved state-of-the-art variational energy and magnetization for quantum antiferromagnetic Heisenberg model.

Demonstrated efficient gradient-based optimization in tensor network algorithms.

Abstract

Differentiable programming is a fresh programming paradigm which composes parameterized algorithmic components and trains them using automatic differentiation (AD). The concept emerges from deep learning but is not only limited to training neural networks. We present theory and practice of programming tensor network algorithms in a fully differentiable way. By formulating the tensor network algorithm as a computation graph, one can compute higher order derivatives of the program accurately and efficiently using AD. We present essential techniques to differentiate through the tensor networks contractions, including stable AD for tensor decomposition and efficient backpropagation through fixed point iterations. As a demonstration, we compute the specific heat of the Ising model directly by taking the second order derivative of the free energy obtained in the tensor renormalization group calculation. Next, we perform gradient based variational optimization of infinite projected entangled pair states for quantum antiferromagnetic Heisenberg model and obtain start-of-the-art variational energy and magnetization with moderate efforts. Differentiable programming removes laborious human efforts in deriving and implementing analytical gradients for tensor network programs, which opens the door to more innovations in tensor network algorithms and applications.

Equations

$$\frac{\partial \mathcal{L}}{\partial \theta }=\frac{\partial \mathcal{L}}{\partial T^n}\frac{\partial T^n}{\partial T^{n-1}}\cdots \frac{\partial T^2}{\partial T^1}\frac{\partial T^1}{\partial \theta }.$$

$$\overline{T^i}=\sum _{j:\mathrm{child}\mathrm{of}i}\overline{T^j}\frac{\partial T^j}{\partial T^i}.$$

$$\overline{A}=U\left[\overline{D}+F\odot (U^T\overline{U}-{\overline{U}}^{T}U)/2\right]U^T,$$

$$\overline{A}$$

$$+U\overline{D}{V}^T+(I-U{U}^T)\overline{U}{D}^{-1}{V}^T+U{D}^{-1}{\overline{V}}^T(I-V{V}^T),$$

$$\overline{A}=\left[\overline{Q}+Q\text{copyltu}(M)\right]R^{-T},$$

$$\overline{A}=\left(\left[(\overline{Q}+\overline{V}{Y}^T)+Q\text{copyltu}(M)\right]U^{-T},Q\overline{V}\right),$$

$$\overline{\theta }=\overline{T^{\ast }}\frac{\partial T^{\ast }}{\partial \theta }$$

$$=\sum _{n=0}^{\infty }\overline{T^{\ast }}{\left[\frac{\partial f(T^{\ast },\theta )}{\partial T^{\ast }}\right]}^n\frac{\partial f(T^{\ast },\theta )}{\partial \theta }.$$

$$\displaystyle Z=\raisebox{-34.44434pt}{\includegraphics[scale={0.3},trim=28.45274pt 28.45274pt 0.0pt 0.0pt,clip]{trg2d.pdf}}.$$

$$\displaystyle T_{uldr}=\raisebox{-12.91663pt}{\includegraphics[scale={0.3},trim=48.36958pt 28.45274pt 0.0pt 0.0pt,clip]{tensor.pdf}}=\frac{\sqrt{\lambda_{u}\lambda_{l}\lambda_{d}\lambda_{r}}}{2}\,\,\delta_{\mathrm{mod}(u+l-d-r,2)},$$

$$H=\sum _{⟨i,j⟩}{S}_i^x{S}_j^x+S_i^y{S}_j^y+S_i^z{S}_j^z.$$

$$\displaystyle A^{s}_{uldr}=\raisebox{-12.91663pt}{\includegraphics[scale={0.3},trim=28.45274pt 28.45274pt 0.0pt 0.0pt,clip]{singlelayer.pdf}},$$

$$\displaystyle T_{uldr}=\raisebox{-12.91663pt}{\includegraphics[scale={0.3},trim=48.36958pt 28.45274pt 0.0pt 0.0pt,clip]{tensor2.pdf}}=\raisebox{-12.91663pt}{\includegraphics[scale={0.3},trim=56.9055pt 28.45274pt 0.0pt 0.0pt,clip]{doublelayer.pdf}}.$$

$$\displaystyle\mathcal{L}={\raisebox{-15.0694pt}{\includegraphics[scale={0.25},trim=51.21504pt 85.35826pt 0.0pt 0.0pt,clip]{energy.pdf}}}\Bigg{/}{\raisebox{-15.0694pt}{\includegraphics[scale={0.25},trim=51.21504pt 85.35826pt 0.0pt 0.0pt,clip]{normalize.pdf}}},$$

Full text

I Introduction

Tensor networks are prominent approaches for studying classical statistical physics and quantum many-body physics problems Orús (2014); Haegeman and Verstraete (2017); Orús (2018). In recent years, its application has expanded rapidly to diverse regions include simulating and designing of quantum circuits Markov and Shi (2008); Arad and Landau (2010); Kim and Swingle (2017); Huggins et al. (2019), quantum error correction Ferris and Poulin (2014); Bravyi et al. (2014), machine learning Stoudenmire and Schwab (2016); Stoudenmire (2018); Han et al. (2018); Cheng et al. (2019); Stokes and Terilla (2019), language modeling Gallego and Orús (2017); Pestun and Vlassopoulos (2017), quantum field theory Verstraete and Cirac (2010); Haegeman et al. (2013); Hu et al. (2018); Tilloy and Cirac (2019) and holography duality Swingle (2012); Hayden et al. (2016).

One of the central problems relevant to many research directions is the optimization of tensor networks in a general setting. Despite highly successful optimization schemes for one dimensional matrix product states White (1992); Vidal (2006); Schollwöck (2011); Stoudenmire and White (2011); Landau et al. (2015); Zauner-Stauber et al. (2018), optimizing tensor networks in two or higher dimensions has been a challenging topic. The hardness is partly due to the high computational cost of tensor contractions, and partly due to the lack of an efficient optimization scheme in the high dimensional situation.

The difficulty is particularly pressing in optimizing tensor network states for infinite translational invariant quantum systems. In these cases, the same tensor affects the variational energy in multiple ways, which results in a highly nonlinear optimization problem. Optimization schemes based on approximate imaginary time projection Jiang et al. (2008); Jordan et al. (2008); Corboz et al. (2010); Zhao et al. (2010); Phien et al. (2015) have been struggling to deal with nonlocal dependence in the objective function. Reference Corboz (2016); Vanderstraeten et al. (2016) apply gradient based optimization and show it significantly improves the results. However, it is cumbersome and error prone to derive the gradients of tensor network states analytically, which involves multiple infinite series of tensors even for a simple physical Hamiltonian. This fact has limited broad adoption of the gradient based optimization of tensor network states to more complex systems. Alternative approaches, such as computing the gradient using numerical derivative has limited accuracy and efficiency, therefore only applies to cases with few variational parameters Poilblanc and Mambrini (2017); Chen et al. (2018). While deriving the gradient manually using the chain rule is only manageable for purposely designed simple tensor network structures Wang et al. (2011).

Differentiable programming provides an elegant and efficient solution to these problems by composing the whole tensor network program in a fully differentiable manner. In this paper, we present essential automatic differentiation techniques which compute (higher order) derivatives of tensor network programs efficiently to numeric precision. This progress opens the door to gradient-based (and even Hessian-based) optimization of tensor network states in a general setting. Moreover, computing (higher order) gradients of the output of a tensor network algorithm offer a straightforward approach to compute physical quantities such as the specific heat and magnetic susceptibilities. The differentiable programming approach is agnostic to the detailed lattice geometry, Hamiltonian, and tensor network contraction schemes. Therefore, the approach is general enough to support a wide class of tensor network applications.

We will focus on applications which involve two dimensional infinite tensor networks where the differentiable programming techniques offer significant benefits compared to the conventional approaches. We show that after solving major technical challenges such as numerical stable differentiation through singular value decomposition (SVD) and memory efficient implementation for fixed point iterations, one can obtain the state-of-the-art results in variational optimization of tensor network states.

The organization of this paper is as follows. In section II we introduce automatic differentiation in the context of tensor network algorithms and formulate tensor network contractions as computation graphs. In section III we present the key techniques for stable and scalable differentiable programming of tensor network algorithms. And section IV we demonstrate the ability of the approach with applications to classical Ising model and quantum Heisenberg model on the infinite square lattice. Finally, we outlook for future research directions opened by this work in Sec. V. Our code implementation is publicly available at Git .

II General Theory

Automatic differentiation through a computation graph is a unified framework that covers training neural networks for machine learning, optimizing tensor networks for quantum physics, and many more. We first review the core idea of automatic differentiation and then explain its application to various tensor network contraction algorithms formulated in terms of computation graphs.

II.1 Automatic differentiation

Automatic differentiation mechanically computes derivatives of computation process expressed in terms of computer programs Bartholomew-Biggs et al. (2000). Unlike numerical differentiation, automatic differentiation computes the value of derivatives to the machine precision. The performance of automatic differentiation has a general theoretical guarantee, which does not exceed the algorithmic complexity of the original program Baur and Strassen (1983); Griewank (1989). Automatic differentiation is the computational engine of modern deep learning applications Rumelhart et al. (1986); Baydin et al. (2015). Moreover, automatic differentiation also finds applications in quantum optimal control Leung et al. (2017) and quantum chemistry calculations such as computing forces Sorella and Capriotti (2010) and optimizing basis parameters Tamayo-Mendoza et al. (2018).

Central to the automatic differentiation is the concept of the computation graph. A computation graph is a directed acyclic graph composed by elementary computation steps. The nodes of the graph represent data, which can be scalars, vectors, matrices, or tensors 111As noted in Laue et al. (2018) one will need tensor calculus as a precise language to present automatic differentiation, even the object being differentiated is a matrix.. The graph connectivity indicates the dependence of the data flow in the computation process. The simplest computation graph is a chain shown in Fig 1(a). Starting from, say vector valued, input parameters $\theta$ one computes a series of intermediate results until reaching the final output $\mathcal{L}$, which we assume to be a scalar. The so called forward evaluation simply traverses the chain graph in sequential order $\theta\rightarrow T^{1}\rightarrow\cdots\rightarrow T^{n}\rightarrow\mathcal{L}$.

To compute the gradient of the objective function with respect to input parameters, one can exploit the chain rule

[TABLE]

Since we consider the case where the input dimension is larger than the output dimension, it is more efficient to evaluate the gradient in (1) by multiplying terms from left to the right using a series of vector-Jacobian products. In terms of the computation graph shown in Fig. 1(a), one traverses the graph backward and propagates the gradient signal from the output back to the input. Computing the derivative this way is called reverse mode automatic differentiation. This approach, commonly referred to as the backpropagation algorithm Rumelhart et al. (1986), is arguably the most successful method for training deep neural networks.

It is instructive to introduce the adjoint variable $\overline{T}={\partial\mathcal{L}}/{\partial T}$ to denote the gradient of the final output $\mathcal{L}$ with respect to the variable $T$. One sees that the reverse mode automatic differentiation propagates the adjoint from $\overline{T^{n}}=\overline{\mathcal{L}}\,\frac{\partial\mathcal{L}}{\partial T^{n}}$ with $\overline{\mathcal{L}}=1$ all the way back to $\overline{T^{i}}=\overline{T^{i+1}}\frac{\partial T^{i+1}}{\partial T^{i}}$ with $i=n-1,\ldots,1$, and finally computes $\overline{\theta}=\overline{T^{1}}\,\frac{\partial T^{1}}{\partial\theta}$. In each step, one propagates the adjoint backward via a local vector-Jacobian product.

The adjoint backpropagation picture generalizes well to more complex computation graphs. For example, the data node $T^{1}$ in Fig. 1(b) affects the final output via two different downstream computation paths. In the backward pass, one needs to accumulate all contributions from its child nodes for its adjoint. In general, the backpropagation rule reads

[TABLE]

The reverse mode automatic differentiation algorithm can be understood as a message passing process on the computation graph. After a topological sort of the computation graph defined by the forward pass, one visits the graph backward from the output node with adjoint $\overline{\mathcal{L}}=1$. Each node collects information from its child nodes to compute its own adjoint, then passes this information to its parents. Thus, one can compute the gradient with respect to all parameters in one forward and one backward pass. Typically one caches necessary information in the forward pass for efficient evaluation of the vector-Jabobian product in the backward pass.

The building blocks of a differentiable program are called primitives. The primitives can be elementary operations such as addition, multiplication, and math functions Griewank and Walther (2008). Each primitive has an associated backward function in addition to the ordinary forward function. The backward function backpropagates the adjoints according to the vector-Jacobian product Eq. (2). Note that one does not need to explicitly instantiate or store the full Jacobian matrices. Moreover, one can group many elementary computation steps together as a primitive. For example, the linear algebra operations such as matrix multiplication can be regarded as a primitive. In this way, the forward pass of these customized primitive can be operated as a black box. Designing the differentiable program in such a modular way allows one to control the level of granularity of automatic differentiation. There are several advantages of crafting customized primitives for domain specific problems. First, this can reduce the computation and memory cost. Second, in some cases, it is numerically more stable by grouping several steps together. Third, one can wrap function calls to external libraries into primitives, without the need to track each individual computation step.

Modern machine learning frameworks support automatic differentiation via various mechanisms. For example, TensorFlow Abadi et al. (2016) explicitly constructs computation graphs using a customized language, autograd Maclaurin et al. (2015), PyTorch Paszke et al. (2017) and Jax Johnson et al. (2018) track the program execution order at run time, and Zygote Innes (2018) performs source code transformation at the compiler level. These frameworks allow one to differentiate through function calls, control flows, loops and many other computation instructions. Building on these infrastructures, differentiable programming is emerging as a new programming paradigm that emphasizes assembling differentiable components and learning the whole program via end-to-end optimization Baydin et al. (2015). Letting machines deal with automatic differentiation mechanically has greatly reduced laborious human efforts and reallocated human attention to design more sophisticated and creative deep learning algorithms.

Finally, we note that it is also possible to evaluate Eq. (1) from right to left, which corresponds to the forward mode automatic differentiation. The operation of the forward mode automatic differentiation is akin to the perturbation theory. One can compute the objective function and the gradient in a single forward pass without storing any intermediate results. However, the forward mode automatic differentiation is not favorable for computation graphs whose input dimension is much larger than the output dimension Baydin et al. (2015). Therefore, the majority of deep learning work employs the reverse mode automatic differentiation. For the same reason, we focus on the reverse mode automatic differentiation for tensor network programs.

II.2 Computation graphs of tensor network contractions

A tensor network program maps input tensors to output tensors, which we assume to be a scalar. Depending on the contexts, the input tenors may represent classical partition function or quantum wavefunction, and the outputs can be various physical quantities of interest. It is conceptually straightforward to apply automatic differentiation to a tensor network program by expressing the computation, and in particular, the tensor network contractions as a computation graph.

As a pedagogic example, consider the partition function of infinite one dimensional Ising model $Z=\lim_{N\rightarrow\infty}{\rm Tr}(T^{N})$ where T=\left(\begin{array}[]{cc}e^{\beta}&e^{-\beta}\\ e^{-\beta}&e^{\beta}\end{array}\right) is a second rank tensor (matrix) representing the Boltzmann weight. One can numerically access the partition function in the thermodynamic limit by repeatedly squaring the matrix $T$ and tracing the final results. The computation graph shows the simple structure as shown in Fig. 1(a). Differentiating with respect to such computational process involves backpropagating the adjoint through matrix trace and multiplication operations, which is straightforward.

At this point, it is also worth distinguishing between the exact and approximate tensor network contraction schemes. Tensor networks on tree like graph can be contracted exactly and efficiently. While other exact approaches, such as the one used for counting and simulating quantum computing Kourtis et al. (2018); Villalonga et al. (2018), in general exhibit exponentially scaling with the problem size. Nevertheless, it is straightforward to apply automatic differentiation to these exact algorithms since they mostly involve tensor reshape and contractions.

We will focus on the less trivial cases of differentiating through approximate tensor contractions, which typically involve truncated tensor factorization or variational approximations. They cover important tensor network applications which show great advantages over other numerical methods Orús (2014); Haegeman and Verstraete (2017); Orús (2018). In particular, we are interested in contracting infinite tensor network, where the fundamental data structure is the bulk tensor. The contraction schemes loosely fall into three categories, the ones based on coarse graining transformations Levin and Nave (2007); Xie et al. (2008); Gu and Wen (2009); Xie et al. (2012, 2014); Evenbly and Vidal (2015); Yang et al. (2017), the ones based on the corner transfer matrix Nishino and Okunishi (1996); Orús (2012); Corboz et al. (2014), and the ones based on matrix product states Vidal (2006); Orus and Vidal (2008); Haegeman and Verstraete (2017). Since the last two contraction schemes are closely related Haegeman and Verstraete (2017), in the following we will focus on automatic differentiation for the tensor renormalization group (Sec. II.2.1) and corner transfer matrix renormalization group approaches (Sec. II.2.2) respectively.

II.2.1 Tensor renormalization group

Tensor renormalization group (TRG) contracts the tensor network by factorizing and blocking the bulk tensors iteratively Levin and Nave (2007). Figure 2(a) shows one step of the TRG iteration as the computation graph, which includes \raisebox{-.9pt} {1}⃝ Split the bulk tensor in two ways using SVD, where we have truncated the singular values and vectors to a prescribed bond dimension $\chi$. \raisebox{-.9pt} {2}⃝ Assemble the four 3-leg tensors generated in the last step into a 4-leg tensor. After this contraction, we obtain a new tensor for the next iteration.

The TRG method grows the lattice size exponentially fast. So one quickly reaches the thermodynamic limit after a few tens of iterations. Note that for numerical stability one needs to rescale the tensor elements after each iteration. The computational cost of TRG method scales $\mathcal{O}(\chi^{6})$ and the memory cost scales as $\mathcal{O}(\chi^{4})$. After unrolling the iterations, the computation graph of the TRG method is similar to the simple chain graph shown in Fig. 1(a). Within each iteration step, the basic operations are tensor index permutation, truncated SVD and tensor contractions. Since each of these operations is differentiable, one can backpropagate through the TRG procedure to compute the derivative of a downstream objective function with respect to the input tensor.

II.2.2 Corner transfer matrix renormalization group

The computation graph of the corner transfer matrix renormalization group (CTMRG) Nishino and Okunishi (1996) has a more interesting topology. The goal of CTMRG calculation is to obtain converged corner and edge tensors which represent the environment degrees of freedom of the bulk tensor.

In cases where the bulk tensor has the full symmetry of the square lattice, the step of one CTMRG iteration is shown in Fig. 2(b). \raisebox{-.9pt} {1}⃝ Contract the bulk tensor with the corner and edge tensors to form a 4-leg tensor. \raisebox{-.9pt} {2}⃝ Perform truncated SVD to the 4-leg tensor, keeping the singular dimensions up to the cut off $\chi$. Keep the truncated singular matrix as the isometric projector. \raisebox{-.9pt} {3}⃝ Apply the isometry to the 4-leg tensor from the first step to find a new corner tensor. \raisebox{-.9pt} {4}⃝ Apply the same isometry to find a new edge tensor for the next step. And iterate this procedure until convergence. One sees that the same bulk tensor with bond dimension $d$ appears in each step of the CTMRG iteration. Due to this reason, the converged environment tensors will depend on the bulk tensor in a complicated way.

Unlike the TRG method Levin and Nave (2007), the CTMRG approach grows the system size linearly. So one may need to iterate a bit more steps to reach convergences in CTMRG. On the other hand, the computational complexity $\mathcal{O}(d^{3}\chi^{3})$ and memory cost $\mathcal{O}(d^{2}\chi^{2})$ of CTMRG are smaller than the ones of TRG in terms of the cutoff bond dimension.

III Technical Ingredients

To compute gradients of a tensor network program using reverse mode automatic differentiation, one needs to trace the composition of the primitive functions and propagate the adjoint information backward on the computation graph. Thankfully, modern differentiable programming frameworks Maclaurin et al. (2015); Abadi et al. (2016); Paszke et al. (2017); Johnson et al. (2018); Innes (2018) have taken care of tracing and backpropagation for their basics data structure, differentiable tensors, automatically.

What one needs to focus on is to identify suitable primitives of tensor network programs and define their vector-Jacobian products for backpropagation. The key components of tensor network algorithms are the matrix and tensor algebras. And there are established results on backward through these operations Giles (2008); Townsend (2016); Seeger et al. (2017). First of all, it is straightforward to wrap all BLAS routines as primitives with customized backward functions. Next, although being less trivial, it is also possible to derive backward rules for many LAPACK routines such as the eigensolver, SVD, and QR factorization Giles (2008). By treating these linear algebra operations as primitives, one can compose a differentiable program with efficient implementations of matrix libraries.

There are, however, a few practical obstacles to stable and scalable implementation of differentiable tensor network programs. First, the backward for the eigensolver and SVD may face numerical instability with degeneracy in the eigenvalues or singular values. Second, the reverse mode automatic differentiation may incur large memory consumption, which prevents one from reaching the same bond dimension of an ordinary tensor network program. We present solutions to these problems in below.

III.1 Stable backward through linear algebra operations

We present several key results on matrix derivatives involving linear algebra operations that are relevant to tensor network algorithms. Recall the modular nature of reverse mode automatic differentiation, one just needs to specify the local backward function to integrate these components into a differentiable program. We will comment on their connections to physics literature and pay special attention to stable numerical implementations Git . For more information, one can refer to Giles (2008); Townsend (2016); Seeger et al. (2017).

III.1.1 Symmetric eigensolver

The forward pass reads $A=UDU^{T}$, where the diagonal matrix $D$ is a diagonal matrix of eigenvalues $d_{i}$ and each column of the orthogonal matrix $U$ is a corresponding eigenvector. In the computation graph the node $A$ has two child nodes $U$ and $D$.

In the backward pass, given the adjoint $\overline{U}$ and $\overline{D}$, we have Giles (2008); Seeger et al. (2017)

[TABLE]

where $F_{ij}=(d_{j}-d_{i})^{-1}$ if $i\neq j$ and zero otherwise. The symbol $\odot$ denotes an element-wise Hadamard product. One can readily check that the gradient is also a symmetric matrix.

Equation (3) can be regarded as "reverse" perturbation theory. When the downstream calculation does not depend on the eigenstate, i.e. $\overline{U}=0$, the backward equation is related to the celebrated Hellmann-Feynman theorem  Feynman (1939) which connects the perturbation to the Hamiltonian and its eigenvalues. Ref. Fujita et al. (2018) applied this special case of Eq. (3) for inverse Hamiltonian design based on energy spectra.

The appearance of the eigenvalue difference in the denominator of $F$ is a reminder of the first order nondegenerate perturbation theory. Reference Tamayo-Mendoza et al. (2018) concerned about the stability of the backward through the eigensolver, thus turned to less efficient forward mode automatic differentiation for variational optimization of Hartree-Fock basis. Actually in many physical problems the final object function depends on part of the eigenvalues and eigenstates in a gauge independent way, e.g., a quadratic form of occupied eigenstates. In these cases, only the eigenvalue difference between the occupied and unoccupied states will appear in the denominator of $F$, which is a familiar patten in the linear response theory Martin (2004). Therefore, degenerated eigenvalues would not necessarily cause problem for these physical applications 222The cases where one has exact degeneracy between occupied and unoccupied eigenstates is a singular case. One should include or exclude all degenerated states in the occupation.. In practice, we found that by using a Lorentzian broadening with $1/x\rightarrow x/(x^{2}+\varepsilon)$ with $\varepsilon=10^{-12}$, one can stabilize the calculation at the cost of introducing a small error in the gradient, see also Seeger et al. (2017).

III.1.2 Singular value decomposition

A ubiquitous operation in tensor network algorithms is the matrix SVD, which is used for canonicalization and factorization of tensor networks Schollwöck (2011); Orús (2014); Haegeman and Verstraete (2017). The forward pass reads $A=UDV^{T}$, where $A$ is of the size $(m,n)$, and $U,V^{T}$ has the size $(m,k)$ and $(k,n)$ respectively, and $k=\min(m,n)$. $D$ is a diagonal matrix contains singular values $d_{i}$. In the reverse mode automatic differentiation, given the adjoints $\overline{U},\overline{D}$ and $\overline{V}$, one can obtain Townsend (2016)

[TABLE]

where $[F_{\pm}]_{ij}=\frac{1}{d_{j}-d_{i}}\pm\frac{1}{d_{j}+d_{i}}$ for $i\neq j$ and zero otherwise. To prevent the numerical issue in case of degenerate singular values, we use the same Lorentzian broadening as Sec. III.1.1 for the first term, which works well in our experience. In practice, for variational tensor network calculation starting from random tensors, the chance of having exact degenerate eigenvalues is small. And even if this happens, applying the rounding is a reasonable solution. While for the cases of degeneracy due to intrinsic reasons Gu and Wen (2009); Pollmann et al. (2010); Xie et al. (2014); Efrati et al. (2014), one will still obtain the correct gradient as long as the end-to-end gradient is well defined. Lastly, inverting the vanishing singular values in Eq. (4) is not a concern since the corresponding space is usually truncated.

III.1.3 QR factorization

QR factorization is often used for canonicalization of tensor networks Schollwöck (2011); Orús (2014); Haegeman and Verstraete (2017). In the forward pass, one factorizes $A=QR$, where $Q^{T}Q=I$ and $R$ is an upper triangular matrix 333One can also require the diagonal of the $R$ matrix to be positive to fix the redundant gauge degrees of freedom. This can be easily implemented on top of the existing QR libraries with backward support.. Depending on the dimensions $(m,n)$ of the matrix $A$ there are two cases for the backward function.

For input shape of $A$ matrix $m\geq n$, $R$ is a $n\times n$ matrix. The backward pass reads Seeger et al. (2017)

[TABLE]

where $M=R\overline{R}^{T}-\overline{Q}^{T}Q$ and the copyltu function generates a symmetric matrix by copying the lower triangle of the input matrix to its upper triangle, $[\texttt{copyltu}(M)]_{ij}=M_{\max(i,j),\min(i,j)}$. The multiplication to $R^{-T}$ can be dealt with by solving a linear system with a triangular coefficient matrix.

For the case of $m<n$, $Q$ is of the size $(m,m)$, and $R$ is a $m\times n$ matrix. We denote $A=(X,Y)$ and $R=(U,V)$, where $X$ and $U$ are full rank square matrices of size $(m,m)$. This decomposition can be separated into two steps, first $X=QU$ uniquely determines $Q$ and $U$, then we calculate $V=Q^{T}Y$. Applying the chain rule, the backward rule gives

[TABLE]

where $M=U\overline{U}^{T}-(\overline{Q}+\overline{V}Y^{T})^{T}Q$.

III.2 Memory efficient reverse mode automatic differentiation with checkpointing function

A straightforward implementation of the reverse mode automatic differentiation for tensor networks has a large memory overhead. This is because one needs to store intermediate results in the forward pass for evaluating the vector-Jacobian products in the backward pass. The number of stored variables is related to the level of granularity in the implementation of the automatic differentiation. In any case, the memory consumption of reverse mode automatic differentiation will be proportional to the depth of the computation graph. This is particularly worrying for tensor networks with large bond dimensions and a large number of renormalization iterations.

The solution to the memory issue of reverse mode automatic differentiation is a well known technique called checkpointing Griewank and Walther (2008); Chen et al. (2016). The idea is to trade the computational time with memory usage. Taking the chain computation graph in Eq. (1) as an example, one can store the tensor every a few steps in the forward process. And in the backward pass, one recomputes intermediate tensors whenever needed by running a small segment of the computation graph forwardly. In this way, one can greatly reduce the memory usage with no more than twice of the computational effort.

In a deeper understanding, checkpointing amounts to define customized primitives which encapsulates a large part of the computation graph. These primitives have their own special backward rules which locally runs the forward pass again and then backpropagates the adjoint. Therefore, in the forward pass one does not need to cache internal states of these checkpointing primitives.

The checkpointing is a general strategy that is applied to the computation graph of any topological structure. When applied to tensor network algorithms, it is natural to regard the renormalization steps shown in Fig. 2 as checkpointing primitives. In this way, one avoids storing some large intermediate tensors in the forward pass.

III.3 Backward through fixed point iteration

Fixed point iteration is a recurring pattern in tensor network algorithms. For example, one iterates the function $T^{i+1}=f(T^{i},\theta)$ until reaching a converged tensor $T^{\ast}$ and uses it for downstream calculations. To compute the gradient with respect to the parameter $\theta$, one can certainly unroll the iteration to a deep computation graph and directly apply the reverse mode automatic differentiation. However, this approach has the drawback of consuming large memory if it takes long iterations to find the fixed point.

One can solve this problem by using the implicit function theorem on the fixed point equation  Christianson (1994). Taking the derivative on both sides of $T^{\ast}=f(T^{\ast},\theta)$, we have

[TABLE]

The second line expands the matrix inversion in the square bracket as a geometric series. Therefore, to backpropagate through a fixed point iteration, the basic operation is just the vector-Jacobian products involving the single step iteration function. And in the backward function one performs iteration to accumulate the adjoint $\overline{\theta}$ until reaching its convergence. The geometric series show the same convergence rate to the fixed point as the forward iteration Christianson (1994).

Many of the tensor network contraction schemes, including the CTMRG method reviewed in Sec. II.2.2, fall into the framework of fixed point iterations. Thus, one can use Eq. (7) for backward through CTMRG calculation, where the iteration is over the RG step shown in Fig. 2(b). We note that the analytical gradient of infinite tensor network contraction derived in Refs. Corboz (2016); Vanderstraeten et al. (2016) contains similar pattern, which is a summation of geometric series.

Similar to the checkpoint technique of Sec. III.2, Eq. (7) also reduces the memory usage in the reverse mode automatic differentiation since one does not need to store a long chain of intermediate results in the forward iteration. Moreover, since the downstream objective function is independent of how the fixed point tensor is obtained, one can also exploit an accelerated iteration scheme Fishman et al. (2018) in the forward process 444A related example is backward through an iterative solver for the dominant eigenstates $x^{\ast}$ of a matrix $A$. In the forward process, one may use the Krylov space method such as Lanczos or Arnoldi iterations. While for the backward pass one can leverage the fact that the dominant eigenvector is the fixed point of the power iteration $f(x,A)=Ax/||Ax||$. Therefore, one can use Eq. (7) to propagate $\overline{x^{\ast}}$ back to $\overline{A}$. Note that the forward and the backward iteration functions are decoupled, thus the backward function does not need to be the same as the forward pass, as long as it ensures $x^{\ast}$ is a fixed point.. There is, however, a caveat when applying Eq. (7) to differentiating tensor network RG algorithms. One may need to pay special attention to the redundant global gauge in the RG iterations to ensure the fixed point equation indeed holds.

III.4 Higher order derivatives

Since the gradient can also be expressed as a computation graph, one can compute the seconder order derivatives by applying automatic differentiation to the graph again. In this way, one can in principle compute arbitrary higher order derivatives of a program using automatic differentiation Griewank and Walther (2008). Deep learning frameworks Maclaurin et al. (2015); Abadi et al. (2016); Paszke et al. (2017); Johnson et al. (2018); Innes (2018) have out-of-the-box support for computing higher order derivatives.

The ability to compute higher derivatives supports Hessian based optimization of tensor network states, such as the Newton method Nocedal and Wright (2006), for tensor network states. However, computing and inverting the full Hessian matrix explicitly could be prohibitively expensive and unnecessary. One can efficiently compute the Hessian-vector product via $\sum_{j}\frac{\partial^{2}\mathcal{L}}{\partial\theta_{i}\partial\theta_{j}}x_{j}=\frac{\partial}{\partial\theta_{i}}\left(\sum_{j}\frac{\partial\mathcal{L}}{\partial\theta_{j}}x_{j}\right)$ without constructing the Hessian matrix explicitly Pearlmutter (1994). This is sufficient for iterative linear equation solvers used for the Newton method.

IV Applications

We present two applications to demonstrate the versatility of differentiable programming tensor network approach for statistical physics and quantum many-body problems. Our public available code implementation Git employs PyTorch pyt with a customized linear algebra automatic differentiation library for improved numerical stability, see discussions in Sec. III.1. While we note that one is readily to reproduce the results with other modern deep learning frameworks such as autograd aut , TensorFlow ten , Jax Jax , and Zygote Zyg frameworks.

IV.1 Higher order derivative of the free energy

Consider an Ising model on the square lattice with inverse temperature $\beta$, its partition function can be expressed as a two dimensional tensor network with bond dimension $D=2$

[TABLE]

The bulk tensor is Chen et al. (2017)

[TABLE]

where $\lambda_{u}=e^{\beta}+(-1)^{u}e^{-\beta}$. We contract the infinite tensor network using the TRG approach discussed in Sec. II.2.1. We use a cut off bond dimension $\chi=30$ and iterate for $30$ TRG steps. Finally, we obtain the partition function Eq. (8) and the free energy by tracing out the bulk tensor.

Next, we compute the physical observables such as energy density and specific heat by directly taking derivatives of the free energy using automatic differentiation, as shown in Fig. 3. One notices that the energy density shows a kink and the specific heat exhibits a peak around the critical temperature $\beta_{c}=\ln(1+\sqrt{2})/2\approx 0.44068679$. Unlike numerical differentiation, these results are free from the finite difference error Xie et al. (2012); Li et al. (2012). Accurate computation of higher order derivatives of the tensor network algorithm will be useful to investigate thermal and quantum phase transitions. We note that it is in principle possible to obtain the specific heat by directly computing the energy variance Vanderstraeten et al. (2015, 2016), which, however, involves cumbersome summation of geometric series expressed in term of tensor networks.

There are alternative ways to compute the specific heat with automatic differentiation. For example, one can directly compute the energy via using the impurity tensor and then take the first order derivative to obtain the specific heat. Or, one can also use forward mode automatic differentiation since there is only one input parameter $\beta$ to be differentiated. We have purposely chosen the present approach to show off the power of differentiable programming with the reverse mode automatic differentiation technique. Backpropagating through the whole TRG procedure, and in particular the SVD, allows one to compute physical observables using higher order derivatives. It is remarkable that this works at all given many of the degenerate singular values due to the $Z_{2}$ symmetry of the Ising model Tamayo-Mendoza et al. (2018). To obtain correct physical results, it is crucial to implement the SVD backward function in a numerical stable way as explained in Sec. III.1.2.

IV.2 Gradient based optimization of iPEPS

We consider a variational study of the square lattice antiferromagnetic Heisenberg model with the Hamiltonian

[TABLE]

We consider an infinite projected entangled pair state (iPEPS) as the variational ansatz. The variational parameters are the elements in the iPEPS

[TABLE]

where $s$ denotes the physical indices, and the remaining indices $u,l,d,r$ are for virtual degrees of freedom of the bond dimension $D$. We initialize the tensor elements with random Gaussian variables. The overlap of the iPEPS forms a tensor network, where the bulk tensor is the double layer tensor with bond dimension $d=D^{2}$

[TABLE]

To contract the infinite tensor network formed by this bulk tensor we use the CTMRG method reviewed in Sec. II.2.2. We initialize the corner and edge tensors by partially tracing out legs from the bulk tensor, then perform the CTMRG iteration until we reach convergence in the corner and edge tensors. After contraction, we can evaluate the expected energy ${\langle\psi|H|\psi\rangle}/{\langle\psi|\psi\rangle}$. Due to the translational invariance of the problem, it is sufficient to consider the expected energy on a bond

[TABLE]

where the black rectangle in Eq. (13) is the Hamiltonian operator acting on a bond. We have performed a basis rotation to the Hamiltonian so that the ground state will have a single site unit cell. We use cutoff bond dimension $\chi=30,50,80,100,144,160$ for $D=2,3,\ldots,7$ respectively. Since the expected energy decreases with the cutoff dimension Poilblanc and Mambrini (2017); Xie et al. (2017), the approximated CTMRG contraction gives variational upper bound to the ground state energy. The expected energy Eq. (13) has both explicit and implicit dependence on the variational parameters in Eq. (11) via the corner and edge tensors.

We compute the gradient of Eq. (13) with respect to the single layer tensor Eq. (11) using automatic differentiation, which automatically resolves the intricate structure in the computation graph Fig. 2(b). The gradient computation takes time comparable to the forward evaluation of the expected energy. Then, we optimize the iPEPS using quasi-Newton L-BFGS algorithm Nocedal and Wright (2006) with the automatically computed gradient. One quickly reaches an optimum after a few hundred function and gradient evaluations. Figure 4(a) shows the relative error in energy compared to extrapolated quantum Monte Carlo (QMC) results Sandvik (2010) for various bond dimensions. The accuracy of the ground state energy is comparable to the state-of-the-art results Corboz (2016); Vanderstraeten et al. (2016), which were shown to be more accurate than imaginary time projection based simple and full update algorithms Jiang et al. (2008); Jordan et al. (2008); Corboz et al. (2010); Zhao et al. (2010); Phien et al. (2015) 555Although it was expected that the full update method will in principle reach the same accuracy as the variational method within the error of Trotter-Suzuki decomposition, it would require optimizing the iPEPS tensors in a globally optimal way.. Note that both Vanderstraeten et al. (2016) and our ansatz contain only half of the variational parameters of the one in Corboz (2016), so the energy results are slightly higher than Ref. Corboz (2016) at $D=2,3$. However, for larger bond dimensions $D=4,5,6,7$, our calculation reaches the lowest variational energy for the infinite square lattice Heisenberg model. Figure 4(b) shows the staggered magnetization measured on the optimized state, which approaches to the extrapolated QMC results at larger bond dimensions.

To obtain results for bond dimension $D>4$, we need to employ either the checkpoint technique III.2 or the fixed point iteration III.3 to keep the memory budget low enough to fit into a single Nvidia P100 GPU card with $12$G memory. It is rather encouraging that with moderate effort one can reach the state-of-the-art performance in variational optimizing iPEPS Corboz (2016); Vanderstraeten et al. (2016). The success is also a nontrivial demonstration that one can indeed stabilize reverse mode automatic differentiation for linear algebra operations appeared in scientific computation Tamayo-Mendoza et al. (2018).

We note that the present approach applies as well to finite systems or problems with larger unit cells, more complex Hamiltonians Corboz et al. (2011); Liao et al. (2017); Lee et al. (2018); Haghshenas et al. (2018), and more sophisticated contraction schemes with improved efficiency Xie et al. (2017), which is promising to deliver new physical results to quantum many-body problems.

V Discussions

Computing the gradient via automatic differentiation significantly boosts the power of existing tensor network algorithms. Researchers can focus on the core tensor network contraction algorithms without worrying about the tedious gradient calculations. The computational complexity of automatic differentiation is the same as the forward contraction of the tensor networks.

Besides greatly reducing human efforts, the automatic differentiation approach also computes a slightly different gradient than Refs. Corboz (2016); Vanderstraeten et al. (2016). The present approach computes numerical exact gradient of an approximated energy density via automatic differentiation. While Refs. Corboz (2016); Vanderstraeten et al. (2016) first derive analytical expression of the energy gradient as infinite tensor networks, and then contract these networks approximately to obtain approximated gradient. Thus, the two approaches perform differentiate the approximation and approximate the derivative respectively Baydin et al. (2015). Other than the general recommendation of Ref. Baydin et al. (2015), we find that differentiating through approximated tensor network contraction can be advantageous for infinite systems whose analytical derivative is complicated to derive and approximate.

In this paper, we have focused on the high level applications of automatic differentiation that differentiates through the whole contraction algorithms for optimizing tensor networks and computing physical observables. The same techniques are also applicable to low level cases such as finding optimal truncation bases or variational transformation of tensor networks Yang et al. (2017). Moreover, besides the optimization of the expected energy of quantum problems, the approach is also relevant to variational contraction of tensor networks Haegeman and Verstraete (2017); Vanderstraeten et al. (2018). We expect differentiable programming techniques will become an integrated part of the standard tensor network toolbox.

A bonus of implementing tensor network programs using deep learning frameworks Abadi et al. (2016); Paszke et al. (2017); Johnson et al. (2018); Innes (2018) is that one can readily enjoy the GPU acceleration. The calculations of this work were done with a single GPU card. Pushing this line of research further, we envision that it will be rewarding to deploy tensor network algorithms on emerging specialized hardware in a larger scale.

Finally, it is useful to comment on the difference of automatic differentiation for tensor networks and neural networks. Typical neural network architectures do not involve sophisticated linear algebra operations. However, with the development of tensorized neural networks Ionescu et al. (2015) and applications of various tensor networks to machine learning problems Novikov et al. (2015); Stoudenmire and Schwab (2016); Stoudenmire (2018); Han et al. (2018); Cheng et al. (2019); Stokes and Terilla (2019), the boundary between the two classes of networks is blurred. Thus, results presented this paper would also be relevant to tensor network machine learning applications when one moves to more sophisticated contraction schemes.

VI Acknowledgment

We are grateful to Pan Zhang, Zhi-Yuan Xie, Wei Li, Ling Wang, Dian Wu, Xiu-Zhe Luo, Shuo-Hui Li, Song Cheng, Jing Chen, Xi Dai, Anders Sandvik, and Frank Verstraete for useful discussions. We thank Philippe Corboz and Laurens Vanderstraeten for providing the reference data shown in Fig. 4. The authors are supported by the National Key Research and Development Project of China Grant No. 2017YFA0302901 and No. 2016YFA0302400, the National Natural Science Foundation of China Grant No. 11888101 and No. 11774398, and the Strategic Priority Research Program of Chinese Academy of Sciences Grant No. XDB28000000.