# A scaling hypothesis for projected entangled-pair states

Bram Vanhecke<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Ghent, Krijgslaan 281, 9000 Gent, Belgium  Juraj Hasik
<EMAIL_ADDRESS>Laboratoire de Physique Théorique UMR5152,
C.N.R.S. and Université de Toulouse,
118 rte de Narbonne, 31062 Toulouse, France  Frank Verstraete  Department of
Physics and Astronomy, University of Ghent, Krijgslaan 281, 9000 Gent, Belgium
Laurens Vanderstraeten  Department of Physics and Astronomy, University of
Ghent, Krijgslaan 281, 9000 Gent, Belgium

###### Abstract

We introduce a new paradigm for scaling simulations with projected entangled-
pair states (PEPS) for critical strongly-correlated systems, allowing for
reliable extrapolations of PEPS data with relatively small bond dimensions
$D$. The key ingredient consists of using the effective correlation length
$\xi$ for inducing a collapse of data points, $f(D,\chi)=f(\xi(D,\chi))$, for
arbitrary values of $D$ and the environment bond dimension $\chi$. As such we
circumvent the need for extrapolations in $\chi$ and can use many distinct
data points for a fixed value of $D$. Here, we need that the PEPS has been
optimized using a fixed-$\chi$ gradient method, which can be achieved using a
novel tensor-network algorithm for finding fixed points of 2-D transfer
matrices, or by using the formalism of backwards differentiation. We test our
hypothesis on the critical 3-D dimer model, the 3-D classical Ising model, and
the 2-D quantum Heisenberg model.

Although variational algorithms over the class of projected entangled-pair
states (PEPS) Verstraete and Cirac (2004); Cirac  _et al._ (2020) are by now
well established for simulating ground states of 2-D quantum spin systems
Jordan  _et al._ (2008); Vanderstraeten  _et al._ (2016); Corboz (2016) and
free energies of 3-D statistical mechanical systems Nishino  _et al._ (2001);
Vanderstraeten  _et al._ (2018), the accurate simulation of critical points
and gapless phases in (2+1)-D remains elusive due to the prohibitive
computational cost of simulating PEPS with a large bond dimension. Indeed,
PEPS are the natural higher-dimensional generalization of matrix product
states (MPS), but, unlike for MPS, expectation values cannot be computed
exactly as that would require the contraction of an infinite two-dimensional
tensor network. Instead, such networks are typically contracted by
constructing effective environments – either as a boundary MPS Verstraete and
Cirac (2004); Fishman  _et al._ (2018) or with corner-transfer matrices
Nishino and Okunishi (1996); Orús and Vidal (2009); Corboz  _et al._ (2011) –
both of which carry a bond dimension $\chi$ controlling the error in the
environment.

The investigation of critical systems by PEPS has been severely hampered by
the presence of these two control parameters. In the case of MPS, it has been
realized that working at finite bond dimension induces a relevant perturbation
in the system, much in the same way as finite sizes enter in exact
diagonalization or Monte-Carlo techniques. As an analogue to finite-size
scaling Brézin (1988); Cardy (1984); Fisher and Barber (1972), the theory of
finite-entanglement scaling Nishino  _et al._ (1996); Pollmann  _et al._
(2009); Tagliacozzo  _et al._ (2008); Pirvu  _et al._ (2012); Vanhecke  _et
al._ (2019a) identifies an induced length scale due to the entanglement
compression caused by finite $\chi$. This has allowed MPS simulations to reach
a high precision in the determination of critical points and exponents, both
for 1-D quantum systems and problems in 2-D statistical mechanics.

Recently, there have been two crucial developments in realizing the same
program for PEPS simulations. First, variational algorithms were devised that
allow to find an optimal PEPS approximation for a given model Hamiltonian
Vanderstraeten  _et al._ (2016); Corboz (2016) or transfer matrix Nishino  _et
al._ (2001); Vanderstraeten  _et al._ (2018); Vanhecke  _et al._ (2019b) at a
certain value of the bond dimension $D$ – in contrast to the simple-update or
full-update schemes that were shown Vanderstraeten  _et al._ (2016) to yield
suboptimal PEPS results. Second, the correlation length of the optimal PEPS at
a certain bond dimension, $D$, was identified as a suitable length scale for
extrapolating variational energies and other observables Rader and Läuchli
(2018); Corboz  _et al._ (2018); Czarnik and Corboz (2019). In both these
advances, the role of the environment bond dimension $\chi$ was not given any
meaningful content: one simply assumes that all optimizations are done for
large enough values of $\chi$, and correlation lengths are then extrapolated
to get rid of any $\chi$ dependence. This makes PEPS simulations very costly,
as the $\chi$-limit is often very tough if not impossible to reach, especially
for large $D$. Moreover, as an implication, only a very limited number of data
points can contribute to the extrapolation – i.e., a single point for every
value of $D$ accessible by computational resources. Contrary to MPS, which can
be optimized up to bond dimensions of order $O(10^{5})$, the gradient-based
PEPS optimizations have so far been limited to $D<10$.

In this paper, we place finite-$\chi$ and finite-$D$ data on an equal footing
in PEPS extrapolation schemes. In order for this to be possible, the optimized
PEPS must be an extreme point of the cost function for a given value of $D$
and $\chi$. That is, we need to solve the optimization problem

$A=\operatorname*{arg\ min}_{A}f_{\\{D,\chi\\}}(A,\bar{A})$ (1)

for tensor(s) $A$ forming the PEPS ansatz where the cost function is the value
of the energy (quantum) or free energy (classical) as a result from the
approximate PEPS contraction _at a fixed value of $\chi$_. In the case of two-
dimensional PEPO transfer matrices, we show that this can be achieved by
gradient optimization, since the $\chi$-dependent gradient of the cost
function is obtained by a simple contraction (see Appendix). For quantum
Hamiltonians, one can construct a complex time evolution PEPOVanhecke  _et
al._ (2019b) and do the same, or the optimization can be done by supplementing
existing contraction algorithms with a backwards differentiation step, as e.g.
done using automatic differentiation Liao  _et al._ (2019); Hasik  _et al._
(2020).

The crucial contribution of this paper is the following hypothesis: close to a
critical point the dependence of expectation values $\braket{O}$ of an
optimized PEPS wavefunction – optimized in the sense of Eq. 1 – as a function
of $D$ and $\chi$ exhibits a data collapse when expressed through the induced
correlation length $\xi$ of the environment –the one used during optimization–
with bond dimension $\chi$:

$\braket{O}(D,\chi)=\braket{O}(\xi(D,\chi)).$ (2)

The consequences of the above hypothesis are twofold. On the one hand, it
implies that even optimized PEPS with relatively small $\chi$ can be used for
accurate extrapolations, provided that they are in the scaling regime. On the
other hand, it allows for a substantial increase in the amount of data
relevant for extrapolation as much more data points can be constructed by
varying $\chi$ for each $D$ considered. We want to stress that this hypothesis
only holds for PEPS algorithms that yield extremal points for the cost
function in Eq. 1, a condition that is not satisfied for simple and full-
update algorithms and for earlier variational algorithms Vanderstraeten  _et
al._ (2016); Corboz (2016), but which is satisfied for the exact
$(D,\chi)$-gradient algorithms such as the one presented in the Appendix and
the ones based on backwards differentiation.

Figure 1: Residual entropy per site for the dimer covering problem on the
cubic lattice. We have optimized PEPS tensors for $D=(2,3,4,5,6)$ and
$\chi=(20,26,34,44,58,75,97,126,164,213,276)$; the same MPS bond dimension was
used for the contraction of the double- and triple-layer (see appendix). The
correlation length is extracted from the double-layer boundary MPS. A fit
(black) on the $D>3$ data reveals a clear $x^{4}$ power law (see inset) for
which the origin is to us an enigma. Figure 2: Magnetization of the classical
3-D Ising model. We have optimized PEPS tensors for $D=(2,3,4)$ and
$\chi=(10,15,20,25,30,40,50)$ for 100 equally spaced temperatures around the
critical temperature. The same MPS bond dimension was used for the contraction
of the double- and triple-layer (see appendix). To plot the data, we shift the
temperature by $T_{c}=4.511528$ (obtained from Monte-Carlo simulations Talapov
and Blöte (1996)) and rescaled both axes with appropriate powers
$\beta=0.326419$ and $\nu=0.629971$ of $\xi$ (obtained from the conformal
bootstrap method Simmons-Duffin (2015)), with $\xi$ the correlation length of
the double-layer MPS environment. In the inset we provide a partial data
collapse (also showing the abundance of data points), showing the crossing of
the rescaled magnetization at the critical temperature.

We first illustrate this hypothesis by scaling the free energy and order
parameters of statistical mechanical problems in three dimensions, where we
optimize PEPS fixed points of the 2-D transfer matrix Vanderstraeten  _et al._
(2018) with a gradient optimization at fixed values of $D$ and $\chi$ (see
Appendix). First, we consider the partition function of dimer coverings on a
cubic lattice. This problem was first tackled with PEPS in Ref. Vanderstraeten
_et al._ , 2018, but here we optimize at fixed values of $D$ and $\chi$
allowing for better convergence of the PEPS optimization in addition to
enabling us to test our PEPS scaling scheme. The data is condensed on a single
plot in Fig. 1. We observe that, indeed, the PEPS data for different
$(D,\chi)$ all lie on the same curve, even for relatively small values of
$\chi$. An extrapolation in $1/\xi^{4}$ (motivated by the inset) gives rise to
an extremely precise value for the residual entropy, $S=-0.449886267(25)$, and
has to be compared with the Monte-Carlo result Beichl and Sullivan (1999) of
$S=-0.4466\pm 0.0006$ – the accurate determination of absolute quantities with
Monte Carlo are well known to be very challenging. Next, we study the 3-D
classical Ising model near criticality, for which we perform a collapse of all
our $(D,\chi$ data for the magnetization, shown in Fig. 2. To perform the data
collapse, we use the critical exponents calculated by bootstrap Simmons-Duffin
(2015) and find the temperature $T_{c}$ that optimizes the collapse:
$T_{c}=4.51170$; this should be compared to $T_{c}=4.511528(6)$ from Monte-
Carlo simulations Talapov and Blöte (1996).

Figure 3: Data for the 2-D Heisenberg model for $D=(6,7,8)$ and selected
$\chi$ from 17 up to 200 for $D=6,7$ and 147 for $D=8$. Top: energy vs.
inverse correlation length of the environment MPS. Bottom-left: energy vs. the
cube of the inverse correlation length. Bottom-right: magnetization vs that
same inverse correlation length. Fits are performed on all the data and are
shown in black.

Finally we consider the square-lattice quantum Heisenberg model. Here, the
PEPS wavefunction for the ground state is completely described by a single
real tensor $A$ with $C_{4v}$ point group symmetry Hasik  _et al._ (2020), so
we can use the original straightforward and highly-efficient formulation of
the corner-transfer matrix algorithm Nishino and Okunishi (1996). In order to
evaluate the gradient of the variational energy, the whole process of
evaluating the energy expectation value is differentiated by backwards
differentiation Liao  _et al._ (2019); Hasik  _et al._ (2020); Hasik and Mbeng
(2020). Again, we optimize PEPS tensors for different $(D,\chi)$, and plot the
energy and magnetization as a function of inverse correlation length in Fig.
3. The staggered magnetization can be extrapolated as a linear function of
$1/\xi$, as was noted in Ref. Rader and Läuchli, 2018, and is here motivated
by a plot in Fig. 3; we find the extrapolated value $m_{s}=0.30771(31)$. The
energy can be extrapolated in $1/\xi^{3}$, as shown in the left-bottom plot of
Fig. 3, to $e=-0.6694401(10)$. These values should respectively be compared to
the Monte-Carlo Sandvik (1997); Sandvik and Evertz (2010) values
$e=-0.669437(5)$ and $m_{s}=0.30743(1)$. For the Heisenberg model, we impose
$U(1)$ symmetry on the PEPS (not on the environment) and select the class that
previously proved to be most effective Hasik  _et al._ (2020).

Notice that the $D=2$ data (blue indicators) for the dimer model are clearly
not in the scaling regime. We have included them as an illustration of the
fact that high enough bond dimensions must be chosen to reach the scaling
regime. For all the above results we implement some filtering of the data,
imposing that the magnetization decreases monotonically in $D$ and $\chi$
separately 111Note that such a filtering on the cost function (residual
entropy for the dimers, free energy for the 3D Ising model, and energy for the
Heisenberg model) is not allowed, as it is not guaranteed to be monotonic in
$\chi$ or $D$, and there are indeed data points with a lower value of the cost
function than the extrapolated value.. This allows us to identify data points
that seemed completely converged, based on the norm of the gradient, but are
in fact stuck at a local minimum and also do not satisfy the scaling
hypothesis. This can be done more thoroughly for the 3-D Ising model, as we
can also impose monotonicity in temperature in addition to $D$ and $\chi$.
Finally, we always use the largest length scale in the boundary MPS. In the
case of Heisenberg model it coincides with the spin-spin correlation length,
but in other models it might prove more fruitful to use some smaller
correlation length associated to a more relevant correlator; this can be
achieved by using symmetric tensor networks or direct fitting of the relevant
correlation function.

In conclusion, we have formulated a scaling hypothesis for PEPS and identified
an induced correlation length $\xi(D,\chi)$ as the only relevant parameter
when the PEPS is optimized using variational gradient methods at fixed bond
dimensions $D$ and $\chi$. This allows for simulating quantum critical points
in 2+1D with much greater accuracy and smaller cost. We conjecture that a
similar scaling hypothesis holds in any dimension, therefore opening up the
possibility of simulating quantum critical points in 3-D and higher.

_Acknowledgments._ — The authors would like to thank Fabien Alet and Andreas
Läuchli for inspiring discussions. This work was supported by the Research
Foundation Flanders (G0E1820N) and the ERC grant QUTE (647905). Part of the
computations were carried out on the HPC resources of CALMIP supercomputing
center under the allocation 2020-P1231. JH is supported by the TNSTRONG
ANR-16-CE30-0025 and the TNTOP ANR-18-CE30-0026-01 grants awarded by the
French Research Council.

## References

  * Verstraete and Cirac (2004) F. Verstraete and J. I. Cirac, “Renormalization algorithms for Quantum-Many Body Systems in two and higher dimensions,” cond-mat/0407066 (2004).
  * Cirac  _et al._ (2020) J. I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete, “Matrix product states and projected entangled pair states: Concepts, symmetries, and theorems,” arXiv:2011.12127 (2020).
  * Jordan  _et al._ (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, “Classical Simulation of Infinite-Size Quantum Lattice Systems in Two Spatial Dimensions,” Phys. Rev. Lett. 101, 250602 (2008).
  * Vanderstraeten  _et al._ (2016) L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, “Gradient methods for variational optimization of projected entangled-pair states,” Phys. Rev. B 94, 155123 (2016).
  * Corboz (2016) P. Corboz, “Variational optimization with infinite projected entangled-pair states,” Phys. Rev. B 94, 035133 (2016).
  * Nishino  _et al._ (2001) T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima, Y. Akutsu, and A. Gendiar, “Two-dimensional tensor product variational formulation,” Progress of theoretical physics 105, 409–417 (2001).
  * Vanderstraeten  _et al._ (2018) L. Vanderstraeten, B. Vanhecke, and F. Verstraete, “Residual entropies for three-dimensional frustrated spin systems with tensor networks,” Phys. Rev. E 98, 042145 (2018).
  * Fishman  _et al._ (2018) M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, “Faster methods for contracting infinite two-dimensional tensor networks,” Phys. Rev. B 98, 235148 (2018).
  * Nishino and Okunishi (1996) T. Nishino and K. Okunishi, “Corner Transfer Matrix Renormalization Group Method,” Journal of the Physical Society of Japan 65, 891 (1996).
  * Orús and Vidal (2009) R. Orús and G. Vidal, “Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction,” Phys. Rev. B 80, 094403 (2009).
  * Corboz  _et al._ (2011) P. Corboz, S. R. White, G. Vidal, and M. Troyer, “Stripes in the two-dimensional t-J model with infinite projected entangled-pair states,” Phys. Rev. B 84, 041108 (2011).
  * Brézin (1988) E. Brézin, “An investigation of finite size scaling,” in  _Finite-Size Scaling_, Current Physics–Sources and Comments, Vol. 2, edited by J. L. Cardy (Elsevier, 1988) pp. 12 – 19.
  * Cardy (1984) J. L. Cardy, “Conformal invariance and universality in finite-size scaling,” Journal of Physics A: Mathematical and General 17, L385–L387 (1984).
  * Fisher and Barber (1972) M. E. Fisher and M. N. Barber, “Scaling theory for finite-size effects in the critical region,” Phys. Rev. Lett. 28, 1516–1519 (1972).
  * Nishino  _et al._ (1996) T. Nishino, K. Okunishi, and M. Kikuchi, “Numerical renormalization group at criticality,” Physics Letters A 213, 69–72 (1996).
  * Pollmann  _et al._ (2009) F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore, “Theory of Finite-Entanglement Scaling at One-Dimensional Quantum Critical Points,” Phys. Rev. Lett. 102, 255701 (2009).
  * Tagliacozzo  _et al._ (2008) L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I. Latorre, “Scaling of entanglement support for matrix product states,” Phys. Rev. B 78, 024410 (2008).
  * Pirvu  _et al._ (2012) B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo, “Matrix product states for critical spin chains: Finite-size versus finite-entanglement scaling,” Phys. Rev. B 86, 075117 (2012).
  * Vanhecke  _et al._ (2019a) B. Vanhecke, J. Haegeman, K. Van Acoleyen, L. Vanderstraeten, and F. Verstraete, “Scaling hypothesis for matrix product states,” Phys. Rev. Lett. 123, 250604 (2019a).
  * Vanhecke  _et al._ (2019b) B. Vanhecke, L. Vanderstraeten, and F. Verstraete, “Symmetric cluster expansions with tensor networks,” arXiv:1912.10512 (2019b).
  * Rader and Läuchli (2018) M. Rader and A. M. Läuchli, “Finite Correlation Length Scaling in Lorentz-Invariant Gapless iPEPS Wave Functions,” Phys. Rev. X 8, 031030 (2018).
  * Corboz  _et al._ (2018) P. Corboz, P. Czarnik, G. Kapteijns, and L. Tagliacozzo, “Finite Correlation Length Scaling with Infinite Projected Entangled-Pair States,” Phys. Rev. X 8, 031031 (2018).
  * Czarnik and Corboz (2019) P. Czarnik and P. Corboz, “Finite correlation length scaling with infinite projected entangled pair states at finite temperature,” Phys. Rev. B 99, 245107 (2019).
  * Liao  _et al._ (2019) H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, “Differentiable programming tensor networks,” Phys. Rev. X 9, 031041 (2019).
  * Hasik  _et al._ (2020) J. Hasik, D. Poilblanc, and F. Becca, “Investigation of the Néel phase of the frustrated Heisenberg antiferromagnet by differentiable symmetric tensor networks,” arXiv:2009.02313 (2020).
  * Talapov and Blöte (1996) A. L. Talapov and H. W. J. Blöte, “The magnetization of the 3D Ising model,” Journal of Physics A Mathematical General 29, 5727–5733 (1996).
  * Simmons-Duffin (2015) D. Simmons-Duffin, “A semidefinite program solver for the conformal bootstrap,” Journal of High Energy Physics 2015, 174 (2015).
  * Beichl and Sullivan (1999) I. Beichl and F. Sullivan, “Approximating the permanent via importance sampling with application to the dimer covering problem,” Journal of Computational Physics 149, 128 – 147 (1999).
  * Hasik and Mbeng (2020) J. Hasik and G. B. Mbeng, “peps-torch: A differentiable tensor network library for two-dimensional lattice models,”  (2020).
  * Sandvik (1997) A. W. Sandvik, “Finite-size scaling of the ground-state parameters of the two-dimensional heisenberg model,” Phys. Rev. B 56, 11678–11690 (1997).
  * Sandvik and Evertz (2010) A. W. Sandvik and H. G. Evertz, “Loop updates for variational and projector quantum monte carlo simulations in the valence-bond basis,” Phys. Rev. B 82, 024407 (2010).
  * Note (1) Note that such a filtering on the cost function (residual entropy for the dimers, free energy for the 3D Ising model, and energy for the Heisenberg model) is not allowed, as it is not guaranteed to be monotonic in $\chi$ or $D$, and there are indeed data points with a lower value of the cost function than the extrapolated value.
  * Vanderstraeten  _et al._ (2019) L. Vanderstraeten, J. Haegeman, and F. Verstraete, “Tangent-space methods for uniform matrix product states,” SciPost Phys. Lect. Notes , 7 (2019).

## Appendix: Differentiating variational cost functions with PEPS

In the main body of the text, we have explained that it is crucial to
differentiate the PEPS cost function evaluated at a fixed value of environment
bond dimension $\chi$. In the context of two-dimensional quantum systems,
differentiating the $\chi$-dependent cost function can be done through the
technique of automatic differentiation, as shown in Refs. Liao  _et al._ ,
2019 and Hasik  _et al._ , 2020. Here, we show that the same procedure is
possible for PEPS approximations of transfer-matrix fixed points, which arise
in the tensor-network representation of 3-D partition functions Vanderstraeten
_et al._ (2018) or the imaginary-time evolution operator for 2-D quantum
Hamiltonians Vanhecke  _et al._ (2019b), under certain symmetry conditions on
the involved tensors.

We start from a uniform projected entangled-pair operator (PEPO) for an
infinite system, described by a 6-leg tensor $T$ as

$O(T)=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=7]{peps1.pdf}}}\quad,$
(3)

which represents a layer of a three-dimensional tensor network. Evaluating
this partition function amounts to finding the leading eigenvector
$\ket{\Psi}$ of this transfer matrix, such that the free energy density is
given by

$f=-\lim_{N\to\infty}\frac{1}{N}\log\left(\frac{\bra{\Psi}O(T)\ket{\Psi}}{\braket{\Psi}{\Psi}}\right),$
(4)

where $N$ is the diverging system size – we will show that we can give meaning
to the density in the thermodynamic limit directly, i.e. without taking the
$N\to\infty$ limit. Assuming this PEPO is Hermitian, we will approximate the
fixed point as a uniform PEPS, parametrized by a single complex-valued tensor
$A$,

$\ket{\Psi(A)}=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=1]{peps1.pdf}}}\quad,$
(5)

and solve the optimization problem

$A=\operatorname*{arg\ min}_{A}f(A,\bar{A}),$ (6)

with

$f(A,\bar{A})=-\frac{1}{N}\log\left(\frac{\bra{\Psi(\bar{A})}T\ket{\Psi(A)}}{\braket{\Psi(\bar{A})}{\Psi(A)}}\right)$
(7)

to find an optimal tensor $A$. Here, $\bar{A}$ is the complex conjugate of the
tensor $A$. In order to solve this optimization problem, we have to evaluate
numerator and denominator in the above expression for a given PEPS, which
boils down to contracting an infinite two-layer and three-layer tensor
network. Here we use boundary MPS methods for the double-layer transfer matrix

$T_{2}=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=4]{peps1.pdf}}}\quad,$
(8)

with

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=5]{peps1.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=6]{peps1.pdf}}}\quad,$
(9)

and triple-layer transfer matrix

$T_{3}=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=12]{peps1.pdf}}}\quad,$
(10)

with

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=11]{peps1.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=10]{peps1.pdf}}}\quad,$
(11)

Here we assume that the tensor $T$ has the symmetry

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=8]{peps1.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=9]{peps1.pdf}}}\quad.$
(12)

with a given unitary matrix $X$, and that we have chosen the PEPS tensor with
the symmetry constraint

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=2]{peps1.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=3]{peps1.pdf}}}\quad$
(13)

with $X$ the same unitary matrix. These symmetry constraints can be
generalised a bit, but the important thing is that the double-layer and
triple-layer transfer matrix are both Hermitian operators. In the examples of
the dimer and Ising model the $X$ matrix is a simple unit matrix.

### Vumps equations

We characterize the boundary MPS as the fixed point of the vumps algorithm;
for more details on uniform MPS and the derivation of the fixed-point
equations, see Ref. Vanderstraeten  _et al._ , 2019. We first list the
equations for a normalized MPS in the center gauge, described by four tensors
$\\{A_{l},A_{r},C,A_{c}\\}$,

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=11]{vumps.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=12]{vumps.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=13]{vumps.pdf}}}\quad,$
(14)

with the conditions

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=19]{vumps.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=20]{vumps.pdf}}}\quad,\qquad\quad\vbox{\hbox{\includegraphics[scale={0.3},page=21]{vumps.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=22]{vumps.pdf}}}\quad$
(15)

and the normalization

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=5]{vumps.pdf}}}\quad=1,\qquad\quad\vbox{\hbox{\includegraphics[scale={0.3},page=14]{vumps.pdf}}}\quad=1.$
(16)

We characterize the fixed point of a given hermitian MPO by the following
fixed-point equations

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=1]{vumps.pdf}}}\quad=\lambda\quad\vbox{\hbox{\includegraphics[scale={0.3},page=2]{vumps.pdf}}}\quad$
(17)

and

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=15]{vumps.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=16]{vumps.pdf}}}\quad.$
(18)

Note that the second equation for $C$ follows from the one for $A_{c}$. The
left and right fixed points are determined by the following eigenvalue
equations

$\displaystyle\quad\vbox{\hbox{\includegraphics[scale={0.3},page=6]{vumps.pdf}}}\quad=\lambda\quad\vbox{\hbox{\includegraphics[scale={0.3},page=7]{vumps.pdf}}}\quad,$
(19)
$\displaystyle\quad\vbox{\hbox{\includegraphics[scale={0.3},page=8]{vumps.pdf}}}\quad=\lambda\quad\vbox{\hbox{\includegraphics[scale={0.3},page=9]{vumps.pdf}}}\quad,$
(20)

with the normalization

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=10]{vumps.pdf}}}\quad=1.$
(21)

Importantly, through the hermiticity of the transfer matrix, we also have the
fixed-point equations for $\bar{A}_{c}$ and $\bar{C}$,

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=3]{vumps.pdf}}}\quad=\lambda\quad\vbox{\hbox{\includegraphics[scale={0.3},page=4]{vumps.pdf}}}\quad,$
(22)

and

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=17]{vumps.pdf}}}\quad=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=18]{vumps.pdf}}}\quad.$
(23)

Given these fixed-point equations and normalizations, the transfer-matrix
eigenvalue per site of the MPO is given by

$\lambda=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=23]{vumps.pdf}}}\quad$
(24)

### Differentiating the fixed-point equations

All tensors
$\\{A_{l},A_{r},A_{c},C,A_{l}^{*},A_{r}^{*},A_{c}^{*},C^{*},G_{l},G_{r}\\}$
can be thought of being functionally dependent on $O$ through the above fixed-
point equations. If we change the transfer-matrix tensor as $O\to O+dO$, all
these tensors have a first-order change as well. We can find these first-order
tensors by differentiating the fixed-point equations. In the following we
denote the first-order tensors with small letters.

We first differentiate the normalization conditions, arriving at the
equalities

$\displaystyle\quad\vbox{\hbox{\includegraphics[scale={0.3},page=1]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=2]{diff.pdf}}}\quad=0,$
(25)
$\displaystyle\quad\vbox{\hbox{\includegraphics[scale={0.3},page=3]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=4]{diff.pdf}}}\quad=0.$
(26)

The normalization of $G_{l}$ and $G_{r}$ can be differentiated to give

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=5]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=6]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=7]{diff.pdf}}}\quad\\\
+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=8]{diff.pdf}}}\quad=0,$
(27)

which, after using the eigenvalue equation for $C$ and $C^{*}$ and the
equation 26, gives

$\quad\vbox{\hbox{\includegraphics[scale={0.3},page=5]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=6]{diff.pdf}}}\quad=0.$
(28)

Next, we find the first-order change in the transfer-matrix eigenvalue,

$d\lambda=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=9]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=10]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=11]{diff.pdf}}}\quad\\\
+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=12]{diff.pdf}}}\quad+\quad\vbox{\hbox{\includegraphics[scale={0.3},page=13]{diff.pdf}}}\quad.$
(29)

The four last terms can be simplified using the eigenvalue equations above
and, therefore, cancel, so that we find

$d\lambda=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=9]{diff.pdf}}}\quad.$
(30)

This implies that the first-order change in the eigenvalue of the transfer
matrix is given by the simple diagram where we differentiate the $O$ tensor.

### Optimizing the free energy

We now return to the objective function that we wish to optimize. Using the
boundary MPS for double- and triple-layer transfer matrix, characterized by
the vumps fixed-point equations, we approximate the free-energy density as

$f(A,\bar{A})=-\log\lambda(A,\bar{A}),$ (31)

and

$\lambda=\quad\vbox{\hbox{\includegraphics[scale={0.3},page=1]{opt1.pdf}}}\quad/\quad\vbox{\hbox{\includegraphics[scale={0.3},page=2]{opt1.pdf}}}\quad.$
(32)

The function $\lambda(A,\bar{A})$ is a real function of complex parameters.
The gradient of the objective function is given by

$g=-\frac{2}{\lambda(A,\bar{A})}\left(\quad\vbox{\hbox{\includegraphics[scale={0.3},page=2]{opt1.pdf}}}\quad\right)^{-1}\left(\quad\vbox{\hbox{\includegraphics[scale={0.3},page=3]{opt1.pdf}}}\quad-\lambda(A,\bar{A})\quad\vbox{\hbox{\includegraphics[scale={0.3},page=4]{opt1.pdf}}}\quad\right)$
(33)

with

$e_{2}=\partial_{\bar{A}}E_{2},\qquad e_{3}=\partial_{\bar{A}}E_{3}.$ (34)

Here, the differentials with respect to the environment tensors all vanish as
we have showed above. This means that the gradient $g$ is the exact gradient
for the objective function evaluated at a given value of the boundary MPS
$\chi$.