Conversion of projected entangled pair states into a canonical form

TL;DR

This paper introduces an algorithm to convert PEPS into a canonical form using a variational gauging approach, enhancing stability and efficiency for tensor network calculations in quantum many-body physics.

Contribution

It presents a novel variational gauging algorithm for PEPS canonicalization, improving convergence and stability in tensor network computations.

Findings

Effective in norm calculations for quantum models

Enables energy optimization within PEPS canonical form

Demonstrates rapid convergence with a new initialization scheme

Abstract

We propose an algorithm to convert a projected entangled pair state (PEPS) into a canonical form, analogous to the well-known canonical form of a matrix product state. Our approach is based on a variational gauging ansatz for the QR tensor decomposition of PEPS columns into a matrix product operator and a finite depth circuit of unitaries and isometries. We describe a practical initialization scheme that leads to rapid convergence in the QR optimization. We explore the performance and stability of the variational gauging algorithm in norm calculations for the transverse-field Ising and Heisenberg models on a square lattice. We also demonstrate energy optimization within the PEPS canonical form for the transverse-field Ising and Heisenberg models. We expect this canonical form to open up improved analytical and numerical approaches for PEPS.

Equations

$$|\Psi ⟩=\sum _{\{s_i\}}\mathcal{F}(A_1^{s_1},A_2^{s_2},\cdots ,A_{l_x\times l_y}^{s_{l_x\times l_y}})|s_1,s_2\cdots ,s_{l_x\times l_y}⟩$$

$${\mathcal{G}}_i=\mathcal{F}\left(\prod _{j\ne i}{E}_j\right)=\text{mathbbm}1,$$

$$H_{\text{ITF}}=-\sum _{⟨ij⟩}{\sigma }_z^i{\sigma }_z^j-\lambda \sum _i{\sigma }_x,$$

$$H_{\text{Heisenberg}}=\sum _{⟨i,j⟩}{\text{S}}_i\cdot {\text{S}}_j,$$

Full text

I Introduction

Tensor network states (TNS) White (1992, 1993); Vidal (2007); Changlani et al. (2009); Mezzacapo et al. (2009) are widely used as variational wave functions to approximate low-energy states of quantum many-body systems Verstraete et al. (2008); Orus (2014). Their power arises from their ability to efficiently capture global behaviors of quantum correlations in the system, as described by entanglement area laws Vidal et al. (2003); Calabrese and Cardy (2004). As a consequence, the global wave function is encoded in local tensors with finite bond dimension. A concrete example is the matrix product state (MPS) Fannes et al. (1992); Östlund and Rommer (1995); Fannes et al. (1994); Schollwöck (2011), a class of tensor-network states that capture the area law in 1D, and which underlie the success of the density matrix renormalization group (DMRG) White (1992, 1993).

The local tensors in a TNS are not uniquely defined and contain redundant parameters known as a local gauge. In MPS, such gauges can be fixed by bringing the MPS into a canonical form where all tensors but one are isometric Schollwöck (2011). The canonical form is simple to compute through QR decompositions, and has many applications, such as in defining optimal local truncations Schollwöck (2005, 2011), the DMRG algorithm, constructing the tangent space of excitations Dorando et al. (2009); Haegeman et al. (2011); Nakatani et al. (2014); Vanderstraeten et al. (2015), and providing a framework to characterize phases Pollmann and Turner (2012); Huang et al. (2014).

Projected entangled pair states (PEPS) Nishino and Okunishi (1996); Verstraete and Cirac (2004); Verstraete et al. (2006); Orus (2014) are higher-dimensional generalizations of MPS with analogous area laws. The PEPS has widely been used as a variational ansatz to explore physical properties of quantum many-body systemsMurg et al. (2009); Wang et al. (2016); Zheng et al. (2017); Liu et al. (2018); Haghshenas et al. (2018); Haghshenas and Sheng (2018); O’Rourke et al. (2018); Haghshenas et al. (2019). It has already been observed that partially fixing the gauge of local tensors can dramatically improve the efficiency and stability of PEPS algorithms Lubasch et al. (2014a); Phien et al. (2015a, b); Evenbly (2018). However, unlike in MPS, computing a fully canonical form for a PEPS remains a challenge.

Here, we introduce a gauging variational ansatz that efficiently brings a PEPS wave function into a full canonical form in direct analogy with that of an MPS. To do so, we re-express the columns of the PEPS as a QR tensor product, where $Q$ is an isometric column tensor network and $R$ is a matrix product operator (MPO). We show that $Q$ can be compactly parametrized by a finite-depth circuit of block isometries and unitaries that can be determined by variational optimization. After transforming all columns but one (a central column) to be isometric, we obtain the (column) canonical form of the PEPS, where part of the entanglement in the environment is transferred to the central column. We explore the stability and performance of the QR decomposition and PEPS canonical representation in calculating the norm in the 2D transverse-field Ising (ITF) and Heisenberg models on a square lattice. We also analyze the behavior of imaginary-time energy optimization in the canonical PEPS form in the context of the ground-state of the ITF and Heisenberg models.

The paper is organized as follows. The basic concepts of the PEPS ansatz are introduced in Sec. II. We first discuss the canonicalization procedure in the context of MPS in Sec. III as a basis to describe our approach to addressing this problem for PEPS. We then study the cost, accuracy, and stability of our gauging variational ansatz in calculations on the ITF and Heisenberg models. We also compare the results of direct energy optimization in the canonical form to results from standard PEPS optimization algorithms. Finally, we summarize our findings in Sec. IV and discuss future research directions.

II PEPS definition and background

A PEPS is a TNS defined by a set of local tensors $\{A^{s_{i}}_{i}\}$ connected by virtual bonds along the grid of the physical lattice. The bond dimension of the virtual bonds is denoted $D$, which controls the number of parameters (or, more physically, the amount of entanglement in the wavefunction) and hence the accuracy of the ansatz. The physical indices $s_{i}$ encode the local physical Hilbert space of dimension $d$. A PEPS wave function $|\Psi\rangle$ on the $l_{x}\times l_{y}=4\times 4$ square lattice with open boundary conditions is depicted in Fig. 1(a),

[TABLE]

where $\mathcal{F}$ denotes tensor contraction of the virtual bonds. The tensors are all colored differently to indicate that we do not assume translational invariance in the tensor network.

The tensor contraction in Eq. (1) is invariant under insertion of a gauge matrix and its inverse $G$, $G^{-1}$ between two tensors (along with a virtual bond). In an MPS, the canonical form at site $i$ is defined as the choice of gauges such that the environment tensor $\mathcal{G}_{i}$ (constructed by partial norm-contraction over all sites except $i$) is the identity tensor

[TABLE]

with $E_{j}=\sum_{s_{j}}A_{j}^{s_{j}\dagger}A_{j}^{s_{j}}$; $\mathbbm{1}$ denotes the tensor $\delta_{i_{1}i_{1}^{\prime}}\delta_{i_{2}i_{2}^{\prime}}\ldots$ where $i_{1}i_{2}\ldots$, $i_{1}^{\prime}i_{2}^{\prime}\ldots$ index the virtual bonds of $A^{s_{i}\dagger}_{i}$, $A^{s_{i}}_{i}$ respectively. By ensuring that the PEPS tensors satisfy Eq. (2), it also defines an analogous canonical form for a PEPS, depicted in Fig. 1(b). In the case of an MPS, we can convert an arbitrary MPS into canonical form by sequential QR (LQ) decompositions of tensors to the left (right) of site $i$, $A_{j}^{s_{j}}\to Q_{j}^{s_{j}}R_{j}$ ($A_{j}^{s_{j}}\to L_{j}Q_{j}^{s_{j}}$) where $Q_{j}$ is orthogonal in the sense $\sum_{s_{j}}{Q_{j}^{s_{j}{\dagger}}}Q_{j}^{s_{j}}=\mathbbm{1}$ (for LQ, $\sum_{s_{j}}{Q_{j}^{s_{j}}}{Q_{j}^{s_{j}{\dagger}}}=\mathbbm{1}$). For simplicity, we henceforth do not distinguish between QR and LQ, with the choice implicit from the diagrammatic representation. $R_{j}$ is then absorbed into the adjacent tensor for the subsequent QR decomposition until the full canonical form is reached Schollwöck (2011).

III PEPS canonical form and column QR ansatz

To similarly canonicalize a PEPS, we sequentially decompose the PEPS columns, denoted $M$ (composed of tensors $\{m_{i}\}$), as a QR tensor contraction, where the column tensor network $Q$ is isometric, satisfying $Q^{\dagger}Q=\mathbbm{1}$, see Fig. 1(c, d). The gauge column tensor network $R$ (composed of tensors $\{r_{i}\}$) is an MPO acting on the horizontal virtual bonds. Once all columns (around a central column) are decomposed to be isometric, the central column $C$ can be viewed as an MPS by grouping the horizontal bonds with the physical bonds. This central column can then be canonicalized around a chosen site using the MPS canonicalization algorithm above, to yield a complete PEPS canonicalization (Fig. 1(e)). Note that the PEPS canonicalization condition around a site (Fig. 1(b)) does not itself specify that columns to the left and right of the central column separately contract to the identity; the conditions we impose are thus sufficient and convenient when canonicalizing a PEPS, but are more constrained than the necessary conditions for Fig. 1(b).

To explicitly carry out the QR decomposition, we first rewrite $M$ and thus $Q$ as MPOs by fusing physical bonds with the left virtual bonds (Fig. 1(f))—for the equivalent LQ decomposition, the physical bonds should be fused with the right virtual bonds, as in MPS. Then, to explicitly enforce the isometric constraint on $Q$, we write it as a finite depth-$n$ circuit of block-size $l$ isometries and unitaries $\{u_{i}\}$, where the isometries appear in the edge layer of the circuit (Fig. 1(g): the isometries have thick bonds (dimension $Dd$) and thin bonds (dimension $D$) while the unitaries only have thick bonds (dimension $Dd$))Pollmann et al. (2016); Wahl et al. (2017). The tensors in the first layer are chosen to be unitary and those in the remaining layers are isometries. The layer depth and block size control the distribution of entanglement between $Q$ and $R$. In practice, to obtain a faithful QR decomposition we have found it sufficient to use $n=2$ (a single layer of unitaries and isometries), increasing $l$ if necessary. In addition, we set the vertical bond dimension of $R$ equal to that of $M$.

To determine the tensors in the QR ansatz, we minimize the distance (cost function) ${F}=\parallel M-Q(\{u_{i}\})R(\{r_{i}\})\parallel$ with respect to $\{u_{i},r_{i}\}$ ($\parallel\cdot\parallel$ is the Hilbert-Schmidt norm) using standard tensor network techniques Evenbly and Vidal (2009); Kao et al. (2015); Evenbly (2017). We optimize the tensors one at a time and sweep until convergence. The cost function depends quadratically on $\{r_{i}\}$, i.e. ${F}=r_{i}^{\dagger}Nr_{i}-Sr_{i}+\textrm{const}$, which is explicitly minimized by solving the linear equation $Nr^{\prime}_{i}=S$. To update the isometric/unitary tensors $\{u_{i}\}$, we observe that the cost function only depends linearly on them due to cancellations, i.e. ${F}=u_{i}^{\dagger}Y+\text{const}$, thus the optimal solution is given by $u^{\prime}_{i}=-VU^{\dagger}$, where $V$, $U$ appear in the singular value decomposition $Y=UsV^{\dagger}$ Evenbly and Vidal (2009). The tensor network diagrams of $N$, $S$ and $Y$ appear in Fig. 2(a, b, c) respectively.

To accelerate the QR optimization (and to avoid local minima) we start with a good initial guess for $\{u_{i},r_{i}\}$. We have used two techniques. The first uses a local projective truncation on the tensors $\{m_{i}\}$ to initialize $\{u_{i}\}$ and $\{r_{i}\}$. To this end, we contract an approximate resolution of the identity (i.e. $\mathbbm{1}_{D}\approx p_{i}p_{i}^{{\dagger}}$ where $p_{i}$ is a local isometry) and a unitary $v_{i}$ into two adjacent tensors $m_{i},m_{i+1}$ and optimize $p_{i}$, $p_{i+1}$ and $v_{i}$ to minimize the local cost function shown in Fig. 2(d). Once we have the optimized tensors $p_{i},p_{i+1}$ and $v_{i}$, we construct a guess for $r_{i},r_{i+1}$ and $u_{i}$ as shown in Fig. 2(d). This initialization is purely local but in practice, we find that it performs well.

A second strategy is based on an accurate estimate of the $\{r_{i}\}$ tensors using a Schulz iteration for the matrix square-root Higham (1997). Note that $R$ is formally the square root of $M^{\dagger}M$, due to the isometric property of $Q$. We thus rewrite $M^{\dagger}M$ as an MPO as in Fig. 2(e). Then, starting from $Y_{0}=M^{\dagger}M$ and $Z_{0}=\mathbbm{1}$, the coupled Schulz iteration, $Y_{k+1}=\frac{1}{2}Y_{k}(3\mathbbm{1}-Z_{k}Y_{k})$, $Z_{k+1}=\frac{1}{2}(3\mathbbm{1}-Z_{k}Y_{k})Z_{k}$, gives $Y_{k\rightarrow\infty}=R$ and $Z_{k\rightarrow\infty}=R^{-1}$. The vertical bond dimension of $Y_{k}$ and $Z_{k}$ increases with each MPO multiplication (Fig. 2(f)) thus we perform MPO compression after each iteration (viewing the MPO as an MPS). The vertical bond dimension of the final $Y_{k}$ ($R$) is compressed back to the original bond dimension of $M$. Also, since $Z_{k}$ approximates $R^{-1}$ which may have arbitrarily large norm, we regularize the iteration using $M^{\dagger}M\rightarrow M^{\dagger}M+\delta I$, where $\delta$ is a small number ($\sim 10^{-6}$). The Schulz iteration converges rapidly (see SM) and we use this accurately estimated $R$ to initialize the optimization of the tensors in $Q$ with respect to the cost function $F$. Although computing the Schulz iteration is more expensive than the local initialization, we expect it to be better when canonicalizing PEPS with more entanglement.

III.1 PEPS canonicalization sweep and truncations

To canonicalize all columns $M^{[1]}M^{[2]}\ldots M^{[l_{x}]}$ in the PEPS, we sweep over all the columns in a prescribed order (say from left to right) and compute the QR decomposition to each. After column $M^{[1]}$ has been converted to QR form, we then absorb the $R$ gauge into the neighboring $M^{[2]}$ column, creating a combined column $C^{[2]}$ with an increased vertical bond dimension of $D^{2}$. To avoid increasing the vertical bond dimension of subsequent columns, we compress $C^{[2]}$ to a smaller vertical bond dimension $D_{c}<D^{2}$. We can perform this column truncation as an MPS truncation with enlarged physical bond dimension $D^{2}d$. The role played by $D_{c}$ is somewhat related to the $\chi$ in PEPS contraction algorithms Corboz et al. (2010a, b); Lubasch et al. (2014a), but here $D_{c}$ is an (auxiliary) bond dimension for a single PEPS layer, rather than for a double layer. A more relevant comparison is therefore to the $\chi$ used in single-layer PEPS algorithms Pižorn et al. (2011); Xie et al. (2017), which is argued to be $\propto D$, thus leading us to conjecture that (asymptotically) $D_{c}\propto D$. Continuing, we perform the QR decomposition on the truncated $C^{[2]}$ column (with vertical bond dimension $D_{c}$), absorb $R$ into $M^{[3]}$ and proceed as before over the remaining columns, to finally produce a PEPS in the canonical form as $Q^{[1]}Q^{[2]}\ldots Q^{[l_{x-1}]}C^{[l_{x}]}$. Note that if MPS truncation is used to compress $C$, then the final central column will be in canonical form around a single site.

Canonicalization redistributes entanglement in the PEPS, thus the canonicalized PEPS has different bond dimensions than the original PEPS. If we use the $n=2$, $l=2$ ansatz for $Q$, then when viewed as an MPO the $Q$ columns have a vertical bond dimension of $\mathcal{O}(D^{2})$, while the central column $C$ has a vertical bond dimension of $D_{c}$. The formally large bond dimension of $Q$ is primarily an artifact of expressing the isometric constraint in terms of gates. Thus it is computationally most efficient to use the structure of $Q$ (i.e. viewing the column of isometries and column of unitaries separately) in the tensor network contractions.

Given a canonical PEPS, the canonicalization sweep can be used to convert between canonical forms (where we move the central column) which is important in algorithms such as energy optimization. The only difference then from the canonicalization sweep discussed above is that the $Q$ columns have a larger vertical bond dimension than $M$. Thus, when absorbing $R$ into a neighboring $Q$ column (in the $n=2$, $l=2$ ansatz), we create a central column $C$ of vertical bond dimension $D^{2}D_{c}$, which we subsequently compress to $D_{c}$.

From the above, we see that in computing the canonical form, and in moving the central column, there are two potential sources of error that must be controlled. The first is the QR approximation error, controlled by the finite-depth/block-size ($n$, $l$) of $\{u_{i}\}$ and the vertical MPO bond dimension $D$ of $R$ (which we fix). The second is the absorption error, that arises from the truncation of the central column’s vertical bond dimension to $D_{c}$.

III.2 Cost of conversion to canonical form

We now discuss the leading costs of the computational steps in the conversion to the canonical form, assuming the $n=2$, $l=2$ ansatz for Q and setting the vertical bond dimension of $R$ to always be equal to that of the column that is being decomposed.

QR optimization. For a column of vertical bond dimension $D$, the cost to determine the isometries/unitaries is $\mathcal{O}(l_{y}D^{6})$ (non-iterative) and $\mathcal{O}(l_{y}D^{4})$, using an iterative algorithm to solve the linear equation ${F}=u_{i}^{\dagger}Y+\text{const}$. The cost to determine $\{r_{i}\}$ tensors is $\mathcal{O}(l_{y}D^{12})$ (non-iterative) and $\mathcal{O}(l_{y}D^{8})$, using a minimization algorithm to solve ${F}=r_{i}^{\dagger}Nr_{i}-Sr_{i}+\textrm{const}$.

Absorption step. As discussed above, when converting a standard PEPS with columns $M^{[1]}M^{[2]}\ldots M^{[l_{x}]}$ that have vertical bond dimension $D$ into canonical form, the typical absorption step during the sweep creates a central column $C$ with an enlarged vertical bond dimension $DD_{c}$ (because $R$ has dimension $D_{c}$). Compressing this down to a vertical bond dimension $D_{c}$, using sequential SVD on the MPS bonds, costs $\mathcal{O}(l_{y}D^{5}D_{c}^{3})$. Alternatively, if we consider direct minimization $||\,|\phi\rangle-|\psi\rangle||$ (the cost typically reported in boundary PEPS algorithms) where $|\phi\rangle$ and $|\psi\rangle$ are MPSs with physical bond dimension $dD^{2}$ and virtual bond dimensions $DD_{c}$ and $D_{c}$ respectively, the cost is $\mathcal{O}(l_{y}D^{4}D_{c}^{4})$.

Alternatively when carrying out the absorption step for a PEPS already in canonical form, e.g. when moving $C$ from $l_{x}$ to $l_{x-1}$ in the PEPS with columns $Q^{[1]}Q^{[2]}\ldots C^{[l_{x}]}$, then the absorption step involves compressing $C^{[l_{x-1}]}$ from vertical bond dimension $D^{2}D_{c}\to D_{c}$. Using sequential SVD on the MPS bonds, the cost is $\mathcal{O}(l_{y}D^{8}D^{3}_{c})$, while direct minimization gives a cost of $\mathcal{O}(l_{y}D^{6}D^{3}_{c})+\mathcal{O}(l_{y}D^{4}D^{4}_{c})$).

However, in this case, since $C^{[l_{x-1}]}=Q^{[l_{x-1}]}R$, we can use the ansatz structure of $Q$ to reduce the cost of the truncation, by absorbing and truncating first the column of isometries, then the column of unitaries. (In both these truncations, the surrounding columns are canonical, and thus each can be performed as an MPS truncation). With this technique, the cost of the absorption step is reduced to $\mathcal{O}(l_{y}D^{4}D^{4}_{c})$ by using direct minimization. The cost is the same for truncating the column of isometries and for the column of unitaries. Using sequential SVD truncation, the leading term is $\mathcal{O}(l_{y}D^{6}D^{3}_{c})$.

III.3 Accuracy of QR ansatz

To assess the accuracy of the QR ansatz, we first study its performance for a single PEPS column. As our initial state, we use the (approximate) ground-state of the spin-$\frac{1}{2}$ ITF and Heisenberg models on the square lattice. The ITF model and the Heisenberg model are respectively defined by

[TABLE]

where ${\textbf{S}_{i}}\equiv(\sigma_{x},\sigma_{y},\sigma_{z})$ and $\sigma_{\alpha}$ are the Pauli matrices. The ITF model has a critical point at $\lambda_{c}\approx 3.05$. Our initial PEPS is constructed from the bulk tensors of an infinite PEPS ground-state Jordan et al. (2008) (optimized with a full-update scheme Corboz et al. (2010b); Phien et al. (2015a) and a $2\times 2$ unit cell) that is repeated periodically across the finite PEPS lattice.

We measure the accuracy of the QR ansatz by the value of its optimization cost function $F$. Here, the parameter controlling the accuracy is the block size $l$ of the isometric/unitary circuit (the number of layers is kept as $n=2$, and the vertical bond dimension of $R$ is kept as $D$). In Fig. 3(a), we show the plot of the distance $F$ versus ITF magnetic field $\lambda$. As expected, when the system is close to criticality, the accuracy is reduced as the ground state becomes more entangled. Increasing $l$ increases the disentangling effect of the unitaries, and the accuracy increases rapidly (we conjecture exponentially with $l$), especially far from criticality.

Next, we investigate the QR accuracy as a function of system size $l_{y}$ for the ITF model at field strength $\lambda=3.5$. As shown in Fig. 3(b), the relative error in $F$ increases linearly with system size, although the slope shows a rapid decay with the isometry/unitary block size $l$. Thus, the variational gauging ansatz introduces a constant error per lattice site, consistent with a fidelity that goes like $e^{-\epsilon l_{y}}\sim 1-\epsilon l_{y}$.

We can also study the accuracy of the Schulz iteration. We show this by evaluating the distance $F_{SC}=\parallel M^{\dagger}M-R^{2}\parallel$ as a function of the Schulz iteration. The MPO compression is controlled by a truncated bond dimension $\chi$ (in the final iteration in the canonicalization algorithm, this is always set to $D$). In Fig. 3(c), we show how the accuracy of the Schulz iteration depends on $\chi$ for different initial bond dimensions $D$ for the ITF model at magnetic field $\lambda=2.0$. The regularization parameter is always set to $\delta\sim 10^{-6}$.

Finally, we give additional results for the QR optimization in the Heisenberg model. We plot the relative distance $F$ versus iteration number for a single bulk column with $l_{y}=10$ in Fig. 3(d), using block size $l=2$. We see that the relative error of the norm contraction $\Delta$ is similar to that of the ITF model at the critical point ($\sim 10^{-2}$) for a lattice of size $10\times 10$.

III.4 Accuracy of PEPS canonical form

We next investigate the accuracy and stability of the full PEPS canonical form constructed from a sweep of the QR approximation and absorption steps across the columns. We estimate the faithfulness of the canonical form from the norm contraction $\mathcal{N}=\langle\Psi|\Psi\rangle$. We compute the norm in the canonical form using only the central column $C$ since all other columns contract exactly to the identity. The relative error of norm contraction is then defined as $\Delta=\frac{\mathcal{N}-\mathcal{N}_{b}}{\mathcal{N}_{b}}$, where the reference value $\mathcal{N}_{b}$ is obtained using an accurate boundary contraction of the original (uncanonicalized) PEPS keeping a large boundary auxiliary bond dimension Verstraete and Cirac (2004); Verstraete et al. (2008); Lubasch et al. (2014b). In Figs. 4(a, b), we show a plot of the relative error $\Delta$ as a function of ITF magnetic field $\lambda$ (using the same approximate ground-state as above) and central column position. Similarly to the single-column results above, the accuracy of the full canonical form depends on the correlation length of the model, and the canonicalization error decreases rapidly (exponentially) as the block size $l$ is increased.

III.5 Norm stability of the canonical form

An important property to assess is the numerical stability of the PEPS in the canonical form. We can study this by examining the stability of the total norm $N$ with respect to small perturbations of the tensors. A numerically stable form is one where small perturbations to the tensors result in small perturbations in the norm. To study this property, we compare the stability of the norm in a standard PEPS versus in its canonical form. We add small perturbations to each tensor $A^{i}_{s_{i}}\rightarrow A^{i}_{s_{i}}+\alpha P$ in a canonical PEPS and a normal PEPS (the unperturbed canonical PEPS is obtained by first converting the normal PEPS to canonical form). The normalized tensors $P$ ($|P|=1$) are chosen to form a uniform random distribution in the interval $[0,1]$. In Fig. 5(a, b), we plot the norm $N$ as a function of perturbation strength $\alpha$ for the ITF ground-state PEPS and a random PEPS. We observe, in both cases, that the canonical PEPS (cPEPS) remains more stable, as the norm changes less with the perturbation as compared to the standard PEPS. This is a result of the fact that the $Q$ tensors are isometric and adding small perturbations to them only affects the isometric property on the unit scale.

III.6 Energy optimization in canonical form

A natural application of the PEPS canonical form is to ground-state energy optimization, which mimics the use of the MPS canonical form in energy optimization. To show this, we perform imaginary time evolution on the ITF model, which we carry out with a sequence of gates $e^{-\tau h}$ on the horizontal and vertical bonds Verstraete et al. (2008); Lubasch et al. (2014a). Evolution on a column of vertical bonds is conveniently carried out on the central column $C$ of a canonical PEPS with bond dimension $D_{c}$, where it reduces to an MPS imaginary time evolution followed by an MPS truncation with an enlarged physical bond dimension $D^{2}d$. Evolution on a column of horizontal bonds can be carried out using a two-column canonical PEPS (analogous to the two-site MPS canonical form) where there are two central columns, and columns to the left and right of these two are isometric tensors $Q$, thus reducing the optimization problem to one of a PEPS with only two columns. In this case, rather than canonicalizing the remaining environment around the bond in the two-column PEPS, we contract it exactly, which is straightforward. Since there are only two columns, these can be reduced to an MPS with enlarged vertical bond dimension $D^{2}D_{c}$ and physical bond dimension $D^{2}d^{2}$. In Fig. 6(a, b) we show the energy as a function of the number of full imaginary time sweeps for the ITF model at field strength $\lambda=3.5$ compared to a near-exact DMRG result. Note that both $D$ (which controls the variational space of the standard PEPS) and $D_{c}$ (vertical bond dimension of the central column, which controls the accuracy of the absorption step in the canonicalization sweep) affect the final converged energy; in this setting, increasing $D_{c}$ has a larger effect than increasing $D$. The relative error of the energy per site reached for the largest bond dimension $D=4,D_{c}=8$ for both $4\times 4$ and $8\times 8$ lattice sizes is on the order of $10^{-4}$.

In Fig. 6(c, d), we benchmark the variational energy for the canonical PEPS and the standard PEPS for the ITF model at $\lambda=3.5$ and the Heisenberg model. For the standard PEPS we used the optimization algorithm from Ref. Lubasch et al., 2014a to obtain the ground state. Both algorithms are initialized by the same random PEPS. The maximum relative errors of the QR and absorption steps using the two layer, $l_{y}=2$ QR ansatz are on the order of $10^{-4}$ and $10^{-3}$ respectively in the ITF and Heisenberg simulation. This relative error bounds the ultimate accuracy in the energy that the canonical PEPS simulation can achieve. Thus we observe that standard PEPS optimization leads to lower energies than the canonical PEPS energy in both cases as expected.

IV Conclusions

In conclusion, we have described a procedure to convert a PEPS into a canonical form analogous to that of an MPS where all columns but one are isometric, by sequentially decomposing columns through a variational QR ansatz. We find that the canonicalization is stable and can be carried out with a small and controllable error. Canonicalization redistributes entanglement in the PEPS, resulting in a central column with increased bond dimension. Our procedure introduces the possibility to formulate canonical PEPS algorithms which make explicit use of an isometric environment, which we demonstrated in an imaginary time optimization of the ground-state energy. The canonical form is clearly numerically more stable, as we show in calculations of the stability of the norm with respect to perturbations of the tensors. However, a faithful comparison of the cost of algorithms using the canonical PEPS and standard PEPS requires considerably more analysis. This is because the canonicalization leads to a non-homogeneous PEPS with different bond dimensions on the vertical and horizontal bonds, quite different from the standard PEPS scenario, and in addition, the canonical PEPS introduces additional numerical parameters (to control the accuracy of the QR ansatz and $R$ absorption) which must be converged. To eliminate the additional errors of the QR decomposition, it may be more expedient to directly optimize the underlying network of isometries and unitaries suggested by the QR ansatz, without explicitly converting into the standard PEPS column-row form. Future investigations will focus on more detailed analysis of these and other algorithms as well as the general representational power of canonicalized PEPS.

Note: As this manuscript was prepared for submission, we were notified of Ref. Zaletel and Pollmann, 2019. In that work, the authors similarly pursue a full canonicalization of a PEPS, but use a different set of isometric conditions on the $Q$ tensors which are more constrained than the ones that we use. Further work is needed to understand the relationship between these techniques.

Acknowledgements.

Primary support for this work was from MURI FA9550-18-1-0095. Some of code used to test energy optimization was based on work supported by the US National Science Foundation (NSF) via grant CHE-1665333. MJO acknowledges a US NSF Graduate Research Fellowship via grant DEG-1745301. GKC acknowledges support from the Simons Foundation. We have used Uni10 Kao et al. (2015) as a middleware library to build up variational gauging ansatz.