Tensor network states and algorithms in the presence of a global U(1) symmetry

Tensor network decompositions offer an efficient description of certain many-body states of a lattice system and are the basis of a wealth of numerical simulation algorithms. In a recent paper [Phys. Rev. A 82, 050301 (2010)] we discussed how to incorporate a global internal symmetry, given by a compact, completely reducible group $\mathcal{G}$, into tensor network decompositions and algorithms. Here we specialize to the case of Abelian groups and, for concreteness, to a U(1) symmetry, associated, e.g., with particle number conservation. We consider tensor networks made of tensors that are invariant (or covariant) under the symmetry, and explain how to decompose and manipulate such tensors in order to exploit their symmetry. In numerical calculations, the use of U(1)-symmetric tensors allows selection of a specific number of particles, ensures the exact preservation of particle number, and significantly reduces computational costs. We illustrate all these points in the context of the multiscale entanglement renormalization Ansatz.

Figure 1

(a) Graphical representation of a tensor $\mathrm{T̂}$ of rank $k$ and components ${\mathrm{T̂}}_{i_1{i}_2\dots i_k}$. The tensor is represented by a shape (circle) with $k$ emerging lines corresponding to the $k$ indices $i_1,i_2,\dots ,i_k$. Notice that the indices emerge in counterclockwise order. (b) Graphical representation of tensors with rank $k=0,1,$ and $2$, corresponding to a complex number $c\in \mathbb{C}$, a vector $|v\rangle \in {\mathbb{C}}^{\left|i\right|}$, and a matrix $\mathrm{M̂}\in {\mathbb{C}}^{\left|i_1\right|\times \left|i_2\right|}$, respectively.

Figure 2

Transformations of a tensor: (a) Permutation of indices $b$ and $c$. (b) Fusion of indices $b$ and $c$ into $d=b\times c$; splitting of index $d=b\times c$ into $b$ and $c$.

Figure 3

(a) Graphical representation of the matrix multiplication of two matrices $\mathrm{R̂}$ and $\mathrm{Ŝ}$ into a new matrix $\mathrm{T̂}$ [Eq. (4)]. (b) Graphical representation of an example of the contraction of two tensors $\mathrm{R̂}$ and $\mathrm{Ŝ}$ into a new tensor $\mathrm{T̂}$ [Eq. (5)].

Figure 4

Graphical representations of the five elementary steps (1)–(5) into which one can decompose the contraction of the tensors of Eq. (5).

Figure 5

(a) Factorization of a matrix $\mathrm{T̂}$ according to a singular value decomposition (11). (b) Factorization of a rank-4 tensor $\mathrm{T̂}$ according to one of several possible singular value decompositions.

Figure 6

(a) Example of a tensor network $\mathcal{N}$. (b) Tensor $\mathrm{T̂}$ of which the tensor network $\mathcal{N}$ could be a representation. (c) Tensor $\mathrm{T̂}$ can be obtained from $\mathcal{N}$ through a sequence of contractions of pairs of tensors. Shading indicates the two tensors to be multiplied together at each step.

Figure 7

Examples of tensor network states for 1D systems: (a) matrix product state (MPS), (b) tree tensor network (TTN), (c) multiscale entanglement renormalization Ansatz (MERA). Examples of tensor network states for 2D systems: (d) projected entangled-pair state PEPS, (e) 2D TTN (2D MERA not depicted).

Figure 8

(a) Tensor $\mathrm{T̂}$ with one incoming index and two outgoing indices, denoted by incoming and outgoing arrows, respectively [Eq. (16)]. (b) A tensor network $\mathcal{N}$ with directed links can be interpreted as a linear map between incoming and outgoing spaces (of the incoming and outgoing indices) obtained by composing the linear maps associated with each of the tensors in $\mathcal{N}$.

Figure 9

(a) Constraint fulfilled by a U(1)-invariant vector. The only allowed particle number on the single index is $n=0$. (b) Constraint fulfilled by a U(1)-invariant matrix. It follows from Schur’s lemma that the matrix is block diagonal in particle number. (c) Constraint fulfilled by a rank-3 tensor with one incoming index and two outgoing indices.

Figure 10

(a) U(1)-covariant vector $|\Psi \rangle$, with some nonvanishing particle number $n\ne 0$. Under the action of U(1) on its index, the covariant vector $|\Psi \rangle$ acquires a phase $e^{-\mathrm{i}n\phi }$ [Eq. (79)]. (b) The U(1)-covariant vector $|\Psi \rangle$, with components $({\Psi }_n){}_{t_n}$, can be represented by a U(1)-invariant matrix $\mathrm{T̂}$ with components ${\mathrm{T̂}}_{i_{1}i}=({\Psi }_n){}_{i_1}$, where $i$ is a trivial index ($\left|i\right|=1$) with charge $n$ and is decorated by the opposite arrow to $i_1$.

Figure 11

A tensor network $\mathcal{N}$ made of U(1)-invariant tensors represents a U(1)-invariant tensor $\mathrm{T̂}$. This is seen by means of two equalities. The first equality is obtained by inserting resolutions of the identity $\mathbb{I}={\mathrm{Ŵ}}_{\phi }{\mathrm{Ŵ}}_{\phi }^{†}$ on each index connecting two tensors in $\mathcal{N}$. The second equality follows from the fact that each tensor in $\mathcal{N}$ is U(1) invariant.

Figure 12

Computation times (in seconds) required to permute and fuse two indices of a rank-4 tensor $\mathrm{T̂}$, as a function of the size of the indices. All four indices of $\mathrm{T̂}$ have the same size, $5d$, and therefore the tensor contains $\left|\mathrm{T̂}\right|=5^4{d}^4$ coefficients. The figures compare the time required to perform these operations using a regular tensor and a U(1)-invariant tensor, where in the second case each index contains five different values of the particle number $n$ (each with degeneracy $d$) and the canonical form of Eqs. (83) and (84) is used. The upper figure shows the time required to permute two indices: For large $d$, exploiting the symmetry of a U(1)-invariant tensor by using the canonical form results in shorter computation times. The lower figure shows the time required to fuse two adjacent indices. In this case, maintaining the canonical form requires more computation time. Notice that in both figures the asymptotic cost scales as $O(d^4)$, or the size of $\mathrm{T̂}$, since this is the number of coefficients which need to be rearranged. We note that the fixed-cost overheads associated with symmetric manipulations could potentially vary substantially with choice of programming language, compiler, and machine architecture. The results given here show the performance of the authors’ matlab implementation of U(1) symmetry.

Figure 13

Computation times (in seconds) required to multiply two matrices (upper panel) and to perform a singular value decomposition (lower panel) as a function of the size of the indices. Matrices of size $5d\times 5d$ are considered. The figures compare the time required to perform these operations using regular matrices and U(1)-invariant matrices, where for the U(1)-invariant matrices each index contains five different values of the particle number $n$, each with degeneracy $d$, and the canonical form of Eqs. (92) and (93) is used. That is, each matrix decomposes into five blocks of size $d\times d$. For large $d$, exploiting the block-diagonal form of U(1)-invariant matrices results in shorter computation time both for multiplication and for singular value decomposition. The asymptotic cost scales with $d$ as $O(d^3)$, while the size of the matrices grows as $O(d^2)$. We note that the fixed-cost overheads associated with symmetric manipulations could potentially vary substantially with choice of programming language, compiler, and machine architecture. The results given here show the performance of the authors’ matlab implementation of U(1) symmetry.

Figure 14

MERA for a system of $L=2\times 3^2=18$ sites, made of two layers of disentanglers $\mathrm{û}$ and isometries $\mathrm{ŵ}$, and a top tensor $\mathrm{t̂}$.

Figure 15

Error in ground-state energy $\Delta E$ as a function of $\chi$ for the $\mathit{XX}$ and Heisenberg models with $2L=108$ spins and periodic boundary conditions, in the particle number sector $N=L$ (or $S_z=0$). The error is calculated with respect to the exact solutions, and is seen to decay exponentially with $\chi$.

Figure 16

Decay of energy gaps $\Delta$ with system size $L$ in the $\mathit{XX}$ model. The upper line corresponds to the energy gap ${\Delta }_L$ between the ground state and the first excited state in the $N=L$ particle number (or $S_z=0$) sector. The lower line corresponds to the energy gap ${\Delta }_{L+1}$ between the ground states of the $N=L$ and $N=L+1$ particle number sectors.

Figure 17

Low-energy spectrum of ${\mathrm{Ĥ}}_{\mathit{XXX}}$ with $L=54$ sites ($=108$ spins). Depicted states have spins of zero ($\times$, blue loops), one (+, red loops), or two ($°$, green loop), and total number of particles ($N$) between 52 and 56. Note that the second and third spin-1 triplets are twofold degenerate.

Figure 18

Computation time (in seconds) for one iteration of the MERA energy minimization algorithm, as a function of the bond dimension $\chi$. For sufficiently large $\chi$, exploiting the U(1) symmetry leads to reductions in computation time. The horizontal line on this graph shows that this reduction in computation time equates to the ability to evaluate MERAs with a higher bond dimension $\chi$: For the same cost per iteration incurred when optimizing a standard MERA in matlab with bond dimension $\chi =20$, one may choose instead to optimize a U(1)-symmetric MERA with partial precomputation and $\chi =24$, or with full precomputation and $\chi =28$.

Figure 19

(a) Graphical representation of the fusion tensor ${\Upsilon }^{\mathrm{fuse}}$ and the splitting tensor ${\Upsilon }^{\mathrm{split}}$. (b) The tensors ${\Upsilon }^{\mathrm{fuse}}$ and ${\Upsilon }^{\mathrm{split}}$ are unitary, and thus yield the identity when contracted pairwise as shown. (c) A fusion tensor decomposed into two parts. The first part (indicated by a circle with an arrow) performs the tensor product of input irreps, $n_A{t}_A\times n_B{t}_B$. The result is an index that labels pairs $(n_A{t}_A,n_B{t}_B)$. The second part (indicated by a rectangle) is a permutation that associates each pair $(n_A{t}_A,n_B{t}_B)$ with a unique $(n_{\mathit{AB}}{t}_{n_{\mathit{AB}}})$, corresponding to a vector in the coupled basis of ${\mathbb{V}}^{(\mathit{AB})}$.

Figure 20

Binary tree decomposition of a symmetric tensor $\mathrm{T̂}$ having components ${\mathrm{T̂}}_{i_1{i}_2{i}_3{i}_4{i}_5{i}_6}$. The tree $\mathcal{T}$ is comprised of a matrix $\mathrm{M̂}$ as the root node, four splitting tensors as internal nodes, and $i_1,i_2,...,i_6$ as its leaf indices. No incoming or outgoing arrows are indicated on the indices in the figure, as the decomposition is valid for any such assignment of directional arrows.

Figure 21

Two possible tree decompositions of a rank-4 tensor $\mathrm{T̂}$. Different choices ${\mathcal{T}}_1,{\mathcal{T}}_2,\dots$ of tree decomposition for tensor $\mathrm{T̂}$ lead to different matrices ${\mathrm{M̂}}_1,{\mathrm{M̂}}_2,\dots$ for the same tensor.

Figure 22

Tree decompositions of tensor $\mathrm{T̂}$ are obtained by contracting the tensor with an appropriate resolution of the identity on its indices, selected according to the desired choice of the fusion tree $\mathcal{T}$. In each instance, evaluation of the contents of the shaded region yields the appropriate matrix $\mathrm{M̂}$.

Figure 23

A matrix ${\mathrm{M̂}}_1$ can be reorganized into another matrix ${\mathrm{M̂}}_2$ by means of fusion tensors, splitting tensors, and the permutation of indices. These operations define a one-to-one linear map $\Gamma$ that acts to reorganize the coefficients of ${\mathrm{M̂}}_1$. $\Gamma$ does not depend on the coefficients of ${\mathrm{M̂}}_1$, but solely on the sequence of operations performed.

Figure 24

Flow diagram for the execution of a predetermined number of iterations of a generic iterative tensor network algorithm, (a) without any precomputation and (b) with precomputation of the operations $\Gamma$.

Figure 25

(a) Permutations applied to one or more legs of a fusion or splitting tensor can be replaced by an appropriate permutation on the coupled index. This process can be used to replace all permutations applied on internal indices of a diagram such as Fig. 23 with net permutations on the indices of ${\mathrm{M̂}}_1$ and on the open indices of the network, as in shown in (b). The residual fusion and splitting operations, depicted as an arrow in a circle, simply perform the basic tensor product operation and its inverse, Eqs. (2) and (3), as described in Fig. 19c and Sec. 1 of this Appendix.