Algorithms for finite projected entangled pair states

Projected entangled pair states (PEPS) are a promising ansatz for the study of strongly correlated quantum many-body systems in two dimensions. However, due to their high computational cost, developing and improving PEPS algorithms is necessary to make the ansatz widely usable in practice. Here we analyze several algorithmic aspects of the method. On the one hand, we quantify the connection between the correlation length of the PEPS and the accuracy of its approximate contraction and discuss how purifications can be used in the latter. On the other hand, we present algorithmic improvements for the update of the tensor that introduce drastic gains in the numerical conditioning and the efficiency of the algorithms. Finally, the state-of-the-art general PEPS code is benchmarked with the Heisenberg and quantum Ising models on lattices of up to $21\times 21$ sites.

Figure 1

A PEPS on a $5\times 5$ square lattice. Physical indices are depicted pointing down, and virtual ones connect tensors in the plane.

Figure 2

Original contraction [1, 16] of the norm TN for a PEPS with bond dimension $D$. The product of a bulk row with a boundary MPO of bond dimension $D^{\prime }$ is approximated by a new boundary MPO of the same dimension, at a cost $\mathcal{O}\left(d{D}^6{D}^{\prime 2}\right)+\mathcal{O}\left(D^4{D}^{\prime 3}\right)$.

Figure 3

Cluster environment of a bulk row for cluster size $\delta =1$. The contraction outside the cluster is achieved by means of a separable positive boundary MPO in the way explained in Ref. [21] with computational cost $\mathcal{O}\left(d{D}^5\right)$. Then the cluster is contracted with the help of a general boundary MPO of bond dimension $D^{\prime }$, which, in this case, can be found with $\mathcal{O}\left(d{D}^5{D}^{\prime 2}\right)$ operations.

Figure 4

Relation between cluster error and correlation function in the Ising model on a $21\times 21$ lattice. (a) Cluster error (main plot) and correlation function (inset) for observable ${\sigma }^Z$, for $D=2$ (open symbols), 3 (filled symbols), and $B=2.0$ (triangles), $2.5$ (squares), $2.8$ (circles). (b) Characteristic cluster size, ${\delta }_0$, vs correlation length, $\zeta$, for several values of $B\in [2,4]$, for observable ${\sigma }^Z$ with $D=2$ (plusses), 3 (crosses), and for ${\sigma }^X$ with $D=2$ (triangles), 3 (circles).

Figure 5

Norm environment for a single site of a $5\times 5$ PEPS. Because each PEPS tensor $A$ from Fig. 1 is contracted with its complex conjugate $A^{*}$ over the physical index, an exact contraction of this norm TN would give a positive Hermitian norm matrix $N$. For large PEPS, an exact contraction of $N$ is not feasible, and the original contraction approximation [1, 16] based on general boundary MPO, shown in Fig. 2, does not keep the positivity.

Figure 6

Positive contraction of the norm TN for a PEPS with bond dimension $D$. The product of a bulk row with a boundary purification of virtual bond dimension $D^{\prime \prime }$ and purification bond dimension $d^{\prime }$ is approximated by a new boundary purification of the same dimensions. Each update of a tensor $A$ constitutes a nonlinear problem, and solving the linearized equations costs $\mathcal{O}\left(d{D}^6{D}^{\prime \prime 4}\right)+\mathcal{O}\left(D^4{D}^{\prime \prime 6}\right)+\mathcal{O}\left(d^{\prime 3}{D}^3{D}^{\prime \prime 6}\right)$ where $d^{\prime }\le D{D}^{\prime \prime 2}$.

Figure 7

Relative error of the norm contraction using boundary purification MPO, for SU Ising ground-state approximations. (a) $B=3.0$ and $D=2$ on lattices of size $21\times 21$ (main plot) and $11\times 11$ (inset). (b) $B=1.0$ and $D=4$ on a $11\times 11$ lattice. For reference, we show the error of the SL method with maximum $d^{\prime }=D{D}^{\prime \prime 2}$ (open circles) and of the original algorithm [1, 16] (filled circles). Our purification contraction was performed with $d^{\prime }=1$ (triangles), 2 (squares), and 3 (diamonds), and $D^{\prime \prime }=\sqrt{D^{\prime }}$.

Figure 8

Six environment tensors of a nearest-neighbor tensor pair. They form a periodic boundary MPO with virtual bond dimension $D^{\prime }$ and physical dimension $D$. It is constructed with $\mathcal{O}\left(d{D}^6{D}^{\prime 2}\right)+\mathcal{O}\left(D^4{D}^{\prime 3}\right)$ operations, resulting from the optimal search for a boundary MPO and the contraction of the environment up to the tensor pair.

Figure 9

A QR decomposition of the left full tensor $A_{\mathrm{L}}$ generates the left reduced tensor $a_{\mathrm{L}}$ as the $R$. Similarly, a LQ decomposition of the right full tensor $A_{\mathrm{R}}$ gives the right reduced tensor $a_{\mathrm{R}}$ as the $L$. The initial $d{D}^4$ variational parameters of the full tensor are decreased to the $d^2{D}^2$ variational parameters of the reduced tensor.

Figure 10

Environment tensor $N_{\mathrm{red}}$ of a reduced nearest-neighbor tensor pair and its closest positive semidefinite approximant $\tilde{X}{\tilde{X}}^{†}$ constructed as explained in the text. The contractions are characterized by the leading computational cost $\mathcal{O}\left(d^4{D}^4{D}^{\prime 2}\right)+\mathcal{O}\left(d^2{D}^6{D}^{\prime 2}\right)+\mathcal{O}\left(d^2{D}^4{D}^{\prime 3}\right)$ and the computation of the positive approximant requires additionally $\mathcal{O}\left(d^6{D}^6\right)$ operations.

Figure 11

Gauge fixing on the environment tensor of the reduced tensor pair, when the environment is nonseparable. (a) We perform a QR and LQ decomposition on $\tilde{X}$ from Fig. 10 independently of each other (notice that we have shortened here the horizontal open indices of $\tilde{X}$ compared to Fig. 10). (b) Contraction of $\tilde{X}$ with $L^{-1}$ and $R^{-1}$ gives the final square root of the environment tensor $X$. (c) In order to leave the state unchanged, the left and right reduced tensors $a_L$ and $a_R$ from Fig. 9 have to be contracted with the gauge transformations $L$ and $R$ as shown here, which gives the starting tensors ${\tilde{a}}_L$ and ${\tilde{a}}_R$ for the update explained in Fig. 12.

Figure 12

Initial and final steps of the reduced tensor update. (a) Initialization: before we update the reduced tensors by sweeping, we apply a SVD to their joint contraction with the Trotter gate, and keep only the $D$ largest singular values in $\Sigma$. Splitting $\Sigma$ gives the initial tensors $a_{\mathrm{L}}:=U\sqrt{\Sigma }$ and $a_{\mathrm{R}}:=\sqrt{\Sigma }V$ for the ALS procedure. (b) Final form: after convergence of the ALS sweeping, we put the two tensors on an equal footing.

Figure 13

Mean value of the relative change ${\epsilon }_d\left(u\right):=|d\left(u\right)-d(u-1)|/|d_{\mathrm{init}}|$ of the cost function $d$, Eq. (1), after consecutive update sweeps $u$ over a tensor pair computed with respect to the initial value of the cost function $d_{\mathrm{init}}$, for the $D=4$ reduced tensor update setting of Table 1. We compare the FU evolution without our gauge fixing using only the positive approximant (open symbols) to the same propagation with our gauge fixing (filled symbols), for a $11\times 11$ Ising model at $B=1.0$ (circles) and $3.0$ (squares), and for a $10\times 10$ Heisenberg model (triangles).

Figure 14

Relative energy error ${\epsilon }_E:=|E\left(D\right)-E_0|/|E_0|$, where $E_0$ denotes the exact ground-state energy from the ALPS library [34, 35, 36], of the SU (open symbols) and the FU (filled symbols) for different lattice sizes. In the case of the SU, we consider $N=10\times 10$ (triangles), $14\times 14$ (squares), $16\times 16$ (diamonds), and $20\times 20$ (circles). In the case of the FU, we consider $N=10\times 10$ (diamonds) and $14\times 14$ (circles).

Figure 15

Spin correlations $C\left(x\right):=|\langle {\vec{S}}_l·{\vec{S}}_{l+x}\rangle |$ (main) and relative error ${\epsilon }_C\left(x\right):=|C\left(x\right)-C_0\left(x\right)|/|C_0\left(x\right)|$ (inset), with the exact values $C_0\left(x\right)$ (thick line) from the ALPS library [34, 35, 36], for two sites separated by distance $x$ along the diagonal in the center of $10\times 10$ PEPS (a) and along the vertical in the center of $14\times 14$ PEPS (b). We consider PEPS Heisenberg ground-state approximations from the FU with $D=2$ (dash-double-dotted), 4 (dash-dotted), 5 (dashed), 6 (filled circles), and 7 (crosses), and from the SU with $D=4$ (squares), 6 (diamonds), and 8 (open circles).

Figure 16

Correlation function $G\left(x\right):=\langle {\vec{S}}_l·{\vec{S}}_{l+x}\rangle -\langle {\vec{S}}_l\rangle ·\langle {\vec{S}}_{l+x}\rangle$ for two sites separated by distance $x$ along the diagonal in the center of $10\times 10$ (main) and $14\times 14$ (inset) PEPS Heisenberg ground-state approximations of bond dimension $D=4$. We compare SU (dotted), CU${}_1$ (dash-double-dotted), CU${}_2$ (dash-dotted), CU${}_3$ (dashed), and FU (solid).

Figure 17

Observables $\langle {\sigma }^Z\rangle$ (main) and $\langle {\sigma }^X\rangle$ (inset) evaluated in the center of $21\times 21$ PEPS Ising ground-state approximations from the FU with $D=2$ (plusses), 3 (crosses), and 4 (filled circles), where the $D=4$ are basically on top of the $D=3$ results. Open symbols show the SU at $B=3.0$, for $D=2$ (down-triangles), 4 (up-triangles), 5 (squares), 6 (diamonds), and 7 (open circles). We interpolate the FU $D=4$ $\langle {\sigma }^Z\rangle$ results between $B=2.85$ and $B_0:=3.000035$ with $|B-B_0{|}^{0.34}$.

Figure 18

Correlation function $G\left(x\right):=\langle {\sigma }_l^Z{\sigma }_{l+x}^Z\rangle -\langle {\sigma }_l^Z\rangle \langle {\sigma }_{l+x}^Z\rangle$ for two sites separated by distance $x$ along the diagonal in the center of PEPS ground-state approximations of the Ising model on a $21\times 21$ lattice with transverse field $B=3.0$ (main) and on a $11\times 11$ lattice with $B=2.8$ (inset). The results were obtained with the FU using bond dimension $D=2$ (plusses), 3 (crosses), and 4 (filled circles), and with the SU using $D=4$ (squares), 6 (diamonds), and 7 (open circles).

Figure 19

Gauge fixing on the environment tensors of the full tensor pair, when the environment is nonseparable. (a) The environment tensor $N_{\mathrm{L}}$ of the left full tensor. (b) We determine the positive approximant for the environment of the left full tensor via a diagonalization of the Hermitian approximant ${\tilde{N}}_{\mathrm{L}}:=(N_{\mathrm{L}}+N_{\mathrm{L}}^{†})/2=U\Sigma U^{†}$, in which, then, the negative eigenvalues are discarded in ${\Sigma }_{+}$, and, finally, the environment is written as ${\tilde{X}}_{\mathrm{L}}{\tilde{X}}_{\mathrm{L}}^{†}$ in terms of its square root ${\tilde{X}}_{\mathrm{L}}:=U\sqrt{{\Sigma }_{+}}$. (c) We perform three independent QR decompositions on ${\tilde{X}}_{\mathrm{L}}$. (d) After equally having carried out the previous steps (a) to (c) with the right full tensor, we have six different matrices $R$. Their inverses are contracted with the corresponding tensors of the boundary MPO.

Figure 20

Like in the reduced case of Fig. 12, the update of the full tensor pair also consists of the three stages initialization, optimization, and final form. The optimization is the standard ALS sweeping, in which each full tensor is gauged after its update in the standard way, i.e., the left tensor is QR decomposed along its right virtual bond and the right tensor is LQ decomposed along its left virtual bond. (a) Initialization I. Firstly, we contract the gauge transformations from Fig. 19 with the full tensors, and split off their reduced parts. (b) Initialization II. Secondly, we construct new reduced tensors from a SVD on the tensor pair and the Trotter gate, equally sharing the $D$ largest singular values between the left and right tensor. (c) Initialization III. We recover the full tensors $A_{\mathrm{L}}$ and $A_{\mathrm{R}}$, which are now the initial tensors for the optimization via ALS sweeping. (d) Final form. After convergence of the sweeping, we put the two tensors on the same footing.