How local is the information in MPS/PEPS tensor networks?

Two dimensional tensor networks such as projected entangled pairs states (PEPS) are generally hard to contract. This is arguably the main reason why variational tensor network methods in 2D are still not as successful as in 1D. However, this is not necessarily the case if the tensor network represents a gapped ground state of a local Hamiltonian; such states are subject to many constraints and contain much more structure. In this paper we introduce a new approach for approximating the expectation value of a local observable in ground states of local Hamiltonians that are represented as PEPS tensor-networks. Instead of contracting the full tensor-network, we try to estimate the expectation value using only a local patch of the tensor-network around the observable. Surprisingly, we demonstrate that this is often easier to do when the system is frustrated. In such case, the spanning vectors of the local patch are subject to non-trivial constraints that can be utilized via a semi-definite program to calculate rigorous lower- and upper-bounds on the expectation value. We test our approach in 1D systems, where we show how the expectation value can be calculated up to at least 3 or 4 digits of precision, even when the patch radius is smaller than the correlation length.


I. INTRODUCTION
Variational tensor-network methods 1 provide a promising way for understanding the low-temperature physics of many-body condensed matter systems. In particular, they seem suitable for studying the ground states of highly frustrated systems, where the sign problem limits many of the quantum Monte Carlo approaches. The best-known and by far the most successful tensor-network method is the Density Matrix Renormalization Group (DMRG) algorithm 2,3 . It has revolutionized our ability to numerically probe strongly correlated systems in 1D, often providing results up to machine precision. Today, it is generally believed to be the optimal method for finding ground states of 1D lattice models.
DMRG can be viewed as a variational algorithm for minimizing the energy of the system over the manifold of Matrix Product States (MPS) 4,5 . These are special types of tensor-network with a linear, 1D structure. In 2D and beyond, several generalizations of this approach are possible. Arguably, the most natural generalization is Projected Entangled Pairs State (PEPS) tensor network 6,7 , which was introduced by Verstraete and Cirac in 2004 6 , but was also used earlier under different names such as "vertex matrix product ansatz" in Ref. 8, "tensors product form ansatz" (TPFA) in Ref. 9, and "tensor product state" (TPS) in Ref. 10. It is the main tensor-network that we consider in this paper. Other constructions include Tree Tensor Networks 11 , Multi-scale Entanglement Renormalization ansatz (MERA) 12 , String bond states 13 , and the recently introduced Projected Entangled Simplex States (PESS) 14 , to name a few. These tensor-networks have proven a vital theoretical tool for understanding the physics of 2D lattice systems and in particular their en-tanglement structure. However, as a numerical method for studying highly frustrated 2D quantum systems, they still face substantial challenges which limit their applicability. In most cases, the best results are still obtained either by DMRG, in which a 1D MPS wraps around the 2D surface, or by quantum Monte Carlo methods.
There are several reasons for this qualitative difference between 1D and 2D systems. The most important one is the computational cost of contracting 2D tensor networks. While this cost scales linearly in the 1D case, it is exponential for 2D and above. Formally, this is reflected in the fact that contracting a PEPS is #P-hard 15 , which is at least NP-hard. To overcome this exponential hurdle, many approximation schemes have been devised. For example, in the original PEPS paper, the network is contracted column by column from left to right, by treating the tensors of a column as matrix product operators (MPO) that act on a MPS that represents the contracted part of the network. Throughout the contraction, the bond dimension of the MPS is truncated to some prescribed D , which introduces some errors (for details, see Refs. 6 and 7). Other approximate methods include, for example, the Corner Transfer Matrix (CTM) [16][17][18][19] , coarse-graining by tensor Renormalization 20-23 , the single layer method 24 and other variants.
While all of the above methods are physically motivated, none of them are rigorous. And to some extent they all produce uncontrolled approximations, even when dealing with the ground state itself. Moreover, while their computational cost is now polynomial in the bond dimension and the particle number, it still scales badly, which limits their practical use to small systems/resolutions (state of the art now days is around 15 × 15 sites with D = 6 25 ). This has led some re-searchers to develop the popular simulation framework, known as the "simple update" method, in which one completely abandons the contraction of the 2D network during the variational procedure 26 , essentially approximating the environment of a local tensor by a product state. While this allows for much higher bond dimensions and is often successful for translationally invariant systems, it may produce poor results for systems that approach critically; see, for example, the analysis in Refs. 25 and 27. In this paper we introduce a new approach for approximating the expectation values of local observables in a 2D PEPS tensor-network. Instead of contracting the full tensor network (or approximating such full contraction), we aim at approximating the expectation value using only the information inside a local patch of the tensor-network around the local observable. Clearly, this approach will fail for general PEPS since they often contain highly nonlocal correlations. However, as we shall demonstrate, such an approach can produce non-trivial results when applied to PEPS that are ground states of gapped local Hamiltonains. Indeed, it is known that such states exhibit strong properties of locality, such as exponential decay of correlations [28][29][30] and local reversibility 31 , and are therefore subject to many constraints to which arbitrary PEPS are not. Moreover, we have the local Hamiltonian at our hands, which can be used to evaluate the local expectation value. Our goal in this paper is to show that for these states, very good approximations for the local expectation values can be derived only from a local patch.
We identify two possible mechanisms for such approximations, which form the basis for two numerical algorithms. In both algorithms, the entire calculation is local and can therefore be done, in principle, efficiently. Moreover, both algorithms provide rigorous upper-and lower-bounds on the expectation value. While we usually cannot give rigorous bound on the distance between these bounds, we demonstrate numerically that this distance -and hence the error in our approximation -can be surprisingly small.
The first algorithm, which we call the 'basic algorithm', is expected to give good results in the case of frustrationfree gapped systems. The second one, which we call the 'commutator gauge optimization' (CGO) algorithm, works only for frustrated systems by utilizing the additional constraints in these systems. It does not rely directly on the existence of a gap, and may work even when considering patches of the PEPS that are much smaller than the correlation length. We numerically demonstrate both algorithms using MPS in 1D systems. Regarding the title of this paper, our findings suggest that when the PEPS/MPS tensor-network represents a frustrated ground state, the information about the local expectation value is largely encoded locally in the neighborhood of the observable.
We would like to stress from the start that the algorithms we present are not practical, and can only be used in 1D in a reasonable time. Their main goal, which is the main goal of this paper, is to illustrate the locality of information in MPS/PEPS representations of gapped ground states. Nevertheless, we strongly believe that the mechanisms behind these algorithms can be further exploited and turned into practical heuristic algorithms. We leave this interesting research direction for further work.
The organization of this paper is as follows. In Sec. II we formally define the problem we wish to solve and the assumptions we are using. In Sec. III we introduce the basic algorithm to solve the problem, and in Sec. IV we introduce the CGO algorithm. In Sec. V we present the results of the 1D numerical tests, and in Sec. VI we discuss the results and offer our conclusions.

II. STATEMENT OF THE PROBLEM
Our construction involves PEPS and MPS tensor networks. For the exact definition and review of these tensor networks (and tensor networks in general), we refer the reader to Refs. 1 and 7.
We are given a local Hamiltonian H = i h i that is defined on a system of N × N spins of local dimension d that sit on a 2D rectangular lattice, which can be either open or with periodic boundary conditions. The local terms h i have O(1) bounded norms, and are assumed to be working on nearest-neighbors only. We further assume that H is gapped, with a unique ground state |Ω that is exactly described by a PEPS with bond dimension D = O(1). While in practice this is rarely the case, and D may have to be exponentially large in order for the PEPS to exactly describe the ground state, it is expected that a D = O(1) bond dimension should give very good approximations to a gapped ground state. For the sake of clarity we first assume that the D = O(1) description is exact and return to discuss this assumption in Sec. IV D.
Given a local observable B that acts on, say, 2 neighboring spins on the lattice, our task is to approximate B := Ω|B|Ω using only a local patch of the PEPS around B. Specifically, the local patch is defined by a ball L of radius around B. L contains all the sites on the lattice that can be connected to the support of B using at most steps on the lattice. We let L c denote the complement region of L that contains all the spins outside of L. An example of the local ball L for = 3 is given in Fig. 1.
Without loss of generality, we will assume that B is normalized so that B = 1; all of our results can be applied to the general case with a simple rescaling.
Consider then the PEPS representation of |Ω , and let α = (α 1 , α 2 , . . . , ) be the set of virtual indices that connect the spins in L to those outside of it. As illustrated in Fig. 2, it defines the following decomposition of |Ω : (1) Above, {|O α }, {|I α } are states that are defined by the PEPS tensor-network outside and inside L respectively. The sets of vectors {|O α }, {|I α } are not necessarily orthogonal between themselves, nor are they normalized, but we assume that both {|I α } and {|O α } are linearly independent among themselves. This is a stronger condition than the injectivity condition for the PEPS, in which only {|I α } have to be linearly independent 32 , but as injectivity, we expect it to hold for generic PEPS. We define V L := span{|I α } and let ρ L denote the reduced density matrix of |Ω on L, which lives in the subspace V L . Note that from our assumptions on {|I α } and {|O α } it follows that V L = Imρ L . We let q := dim V L , which is equal to the size of the range in which the composite index α = (α 1 , α 2 , . . .) runs. Notice that q = D O(|∂L|) , reflecting the fact that |Ω satisfies an area law.
Clearly, if we knew ρ L we could calculate B from B = Ω|B|Ω = Tr L (ρ L B). However, directly estimating ρ L involves the contraction of the full network, which is exactly what we are trying to avoid. Instead, since the radius is assumed to be small, we can easily calculate V L := span{|I α }. Can we estimate B using only V L and the fact that |Ω is the ground state of H? Formally, we ask: Problem 1 Given a local observable B, a ball L of radius around it, and the corresponding subspace V L , find a range [b min , b max ], as narrow as possible, such that B ∈ [b min , b max ].
We offer two algorithms to tackle this problem. The first, which we call the 'basic algorithm', is expected to work well when the system is gapped, and is suitable mostly for frustration-free systems. The second, which we call 'the commutators gauge optimization' algorithm, is much more powerful and applies only for frustrated systems. It relies on the assumption that |Ω with D = O(1) is described by a PEPS and therefore satisfies an area law (this, in turn, can also be attributed to the existence of a gap). We begin with the basic algorithm.

III. THE BASIC ALGORITHM
Let P V be the projector into the subspace V L , and define The basic algorithm then simply calculates the largest and the smallest eigenvalues of B L in the V L subspace and use these as lower-and upper-bounds to B .
Intuitively, when the system is gapped, we expect the range [b 1 , b q ] to narrow down as → ∞ for the following reason. Every α = (α 1 , . . . , α q ) can be seen as a specific boundary condition for a restricted system on the ball L, and as the system is gapped, we expect the effect of these boundary conditions to be negligible. Hence, we expect that for every |I α , we have I α |B|I α / I α 2 → B as → ∞, and similarly for any linear combination of the |I α vectors. In particular this means that the eigenvalues of B L are expected to converge to B . The following lemma shows that this is indeed the case when the eigenvalues of ρ L are not too small.
Lemma 1 Let λ 2 min be the minimal eigenvalue of ρ L , and let b 1 , . . . , b q be the eigenvalues of B L = P V BP V in the subspace V L . Then for every β = 1, . . . , q, where ξ > 0 is the correlation length of |Ω .
Proof: The proof is a simple application of the basic idea in the Choi-Jamio lkowski isomorphism, together with the exponential decay of correlation property of gapped ground states 28 . Let |b 1 , |b 2 , . . . |b q be the eigenvectors of B L in V L that correspond to the eigenvalues b 1 , b 2 , . . . , b q . They constitute an orthonormal basis of V L . Let |Ω = α λ α |Ô α ⊗ |Î α be the Schmidt decomposition of |Ω with respect to L and L c such that Both {|b β } and {|Î α } are orthonormal bases of V L , and therefore they are connected by a unitary transformation: |b β = α U βα |Î α . We use it to define a set of (un-normalized) states on the spins of L c : In addition, we define the operator C β := |c β c β | ⊗ 1 L , where 1 L is the identity operator on the spins of L. Using these definitions, the following identities are easy to verify: Next, by using the exponential decay of correlations 28 and the uniqueness of the ground state, we deduce that Let us analyze the above inequality term by term. First, by Eq. (6) and the fact that C β and B commute, A similar argument shows that Ω|C β |Ω = 1. Substituting this into the LHS of (7), we get |b β − B | ≤ B · C β ·e − /ξ . Finally, using Eq. (5) and the fact that γ |U βγ | 2 = 1 (since U βγ are the entries of a unitary matrix), proves the lemma.
Lemma 1 shows that if the smallest eigenvalue of ρ L is lowerbounded by e −c /ξ for some c < 1, the range [b 1 , b q ] should converge to the point B exponentially fast in . This is certainly expected in most 1D ground states that are approximated by an MPS of a constant bond dimension D. However, it can hardly happen in 2D, where even if the state is approximated by a PEPS with D = O(1), the dimension of V L grows at least as D O( ) , and consequently, the smallest eigenvalue of ρ L is expected to fall off exponentially like D −O( ) . In 3D things are even worse, since dim V L grows like D O( 2 ) .
When the underlying Hamiltonian is frustration-free, V L is a subspace of the groundspace of H L , the part of H that is supported on the spins of L. In such case, it is easy to see that if the Hamiltonian satisfies the local topological order (LTQO) property 33 (see also Ref. 34), then necessarily |b α − B | decays as a function of (the rate of the decay depends on the particular definition of LTQO). Related to that, Lemma 1 has much in common with Theorem 11 of Ref. 34, which shows that parent Hamiltonians of translationally invariant, injective MPS satisfy LTQO.
We conclude this section with a somewhat stronger result than Lemma 1 that holds for frustration-free systems: Lemma 2 When |Ω is the unique ground state of a frustration-free Hamiltonian with an O(1) spectral gap, then where |δ ≤ e −O( ) . In other words, |Ω is an approximate eigenvector of B L .
The proof of this lemma is given in the Appendix. Formally, it is at least as strong as Lemma 1 because it can be used to derive Lemma 1 when the system is frustrationfree. Indeed, using the notation of the proof of Lemma 1, the lemma can be proved by multiplying inequality (8) by C β and then by c β | ⊗ β β |. But it also seems stronger since it tells us something about the eigenvalues of B L even when λ 2 min is much smaller then e − /ξ . To see this, note that (8) implies that If {|b β } is the eigenbasis of B L with eigenvalues b 1 , b 2 , . . ., then the above inequality can be written as As { b β |ρ L |b β } is the probability distribution that corresponds to the measurement of B L , we can use Markov inequality to deduce that if we measure B L on |Ω , then with probability of at least 1 − e −O( ) we will obtain an eigenvalue of B L that is e −O( ) close to B . This holds regardless of the minimal eigenvalue of ρ L . Therefore in cases where the weight of all the eigenvalues of B L in ρ L is of the same order (which can be much much smaller than e −O( ) in 3D), then an exponentially large fraction of them must be exponentially close to B . This does not prove that the range [b min , b max ] rapidly shrinks, but it supports the intuition that this should generally be the case.
We conclude this section noting that for frustrated systems a similar result to Lemma 2 can be proved using the same techniques that are used to prove the exponential decay of correlations in the frustrated case 28 : a combination of Lieb-Robinson bounds and a suitable filtering function. In fact, a very similar result was already proven by Hastings in Ref. 35 for the special case when B is a local Hamiltonian term. However, instead of pursuing that direction, we turn to a different algorithm, which turns out to be much more powerful for our problem.

IV. THE COMMUTATOR GAUGE OPTIMIZATION
In this section we introduce the Commutator Gauge Optimization algorithm (CGO for short), which is appli-cable for frustrated systems. As we shall see, it can be viewed as pair of primal-dual SDP optimization problem. We start with the formulation of the primal optimization problems.

A. Primal problem
Our starting point is the simple observation of the basic algorithm, that the expectation value B must be inside eigenvalues range [b min , b max ] of the operator B L := P V BP V . The additional idea is that this must hold for any other operator K for which K = B . Therefore, if we can find an operator K such that B = K , then necessarily B must also be inside the range [k min , k max ] of the minimal and maximal eigenvalues of P V KP V . By optimizing over a subset of such operators, we may significantly narrow down the range in which B is found. We therefore look for operators, K min , K max on L such that Ω|K min |Ω = Ω|K max |Ω = Ω|B|Ω , but at the time, P L K max P L has the smallest possible maximal eigenvalue, and P L K min P L has the largest possible minimal eigenvalue. How can we find such operators? A very simple yet powerful trick is to look for operators of the form: where A min , A max are anti-Hermitian such that K min , K max are Hermitian. For any eigenstate | of H and for any operator A, it is easy to verify that |[H, A]| = 0 and therefore Eq. (10) holds. In that sense K min and K max can be viewed as a different "gauges" of B, which have the same expectation value with respect to |Ω .
By taking the support of A min , A max to be small enough, we can guarantee that [H, A min ] and [H, A max ] are supported in L. Formally, we partition the spins of L into two disjoint subsets, L = ∂L∪L 0 . Here, ∂L contains the spins in L that are coupled to spins in L c via one or more local terms in H, and L 0 contains the rest of the spins in L. An illustration of this decomposition is given in Fig. 3. Letting H L denote the sum of all the local h i terms whose support is inside L, we conclude that when the support of A is inside radius around B, calculate and k (P ) Above, the optimizations are over all anti-Hermitian operators A that are supported on L 0 .
Two notes are in order. First, the presence of P V in the above equations is crucial, as it is the only place where information about the ground state is entering the problem. Without it, values from other eigenstates of H might worsen our estimate. It is therefore clear why some sort of an area-law is required for CGO to work; if the support of ρ L (i.e., the subspace V ) is the entire available Hilbert space, there would be nothing in V L to tell us that we are dealing with the ground state and not, say, some other eigenstate of H. Second, when H is frustration free, P V H L = H L P V = 0 and so P V (B + [H L , A])P V = P V BP V for all A; the CGO algorithm then simply becomes the basic algorithm.
Before discussing how the optimization Problem 2 can be done in practice, let us introduce its dual.

B. The dual problem
To introduce the dual problem, we begin with the constraint Ω|[H L , A]|Ω = 0, which holds for every A that is supported on the spins in L 0 . Tracing out the spins in L c , we get Using But since A is an arbitrary anti-Hermitian operator on L 0 , it implies that the expression it multiplies has to vanish. We therefore reach the following corollary: Corollary 3 Let ρ L be the reduced density matrix of an eigenstate of H in the region L. Then, The above local identity holds for all eigenstates of H and for all regions L, regardless of any assumption about a gap or the existence of a PEPS description. When H is frustration free and ρ L is the reduced density matrix of the ground state, it is satisfied trivially since [ρ L , H L ] = 0. However, when the system is frustrated, Eq. (15) provides many non-trivial constraints as it holds pointwise on the Hilbert space of L 0 . Corollary 3 allows us to formulate a dual to the optimization problem given in Problem 2. Instead of optimizing over anti-Hermitian operators A, we may optimize over 'legal' reduced density matrices, which satisfy Eq. (15).
One can prove directly that k max and that k (P ) min . We omit the proof since it will follow naturally from viewing these two optimizations as dual semi-definite programs (SDPs), as we shall now do.

C. The commutator gauge optimization as a SDP
The primal CGO problem and its dual can be cast as semi-definite programs. Let us demonstrate it for the upperbound of B ; the lowerbound follows similarly. To write (12) as an SDP, note that a number λ is an upperbound to all eigenvalues of an operator O iff λ1 − O 0, where 0 stands for non-negative matrix. Then choosing a basis A 1 , A 2 , . . . , A m for all the anti-Hermitian matrices in L 0 , the optimization over (12) can be written as Above, the optimization is done over the real numbers b, c 1 , c 2 , . . . , c m . The dual SDP is then maximize: Tr(ρ L P V BP V ) (20) subject to: ρ L 0, Tr(ρ L 1) = 1, (21) and Tr(ρ L P V [H L , A i ]P V ) = 0 ∀i which is identical to the optimization in (16). By the weak duality of semi-definite programming, we note that the value in (20) lowerbounds the one in (18), and that in many cases they are identical. In the next section, we present the results of some numerical tests we performed on 1D systems to check the performance of such SDPs.

D. Working with constant bond dimensions
Up to now, we have assumed that the ground state |Ω is given exactly as a PEPS with constant bond dimension D. However, in virtually all frustrated systems, a constant D can only give an approximation to the ground state. To distinguish the actual PEPS approximation from the exact ground state, we will denote it by |Ω p . Given that |Ω p is only an approximation to the ground state |Ω , we can only expect Ω p |[H, A]|Ω p to be close to zero, but not completely vanish. Similarly, Eq. (15) is expected to be only approximately satisfied. While this may look merely as an aesthetic defect, it raises a conceptual problem: if D is constant, the dimension of V L increases like D O( ) , while the number of constraints in Eq. (15) goes like d O( 2 ) . For large enough, yet constant , it is expected that Eq. (15) is over-determined, or simply unsatisfiable. From an SDP point of view, we expect the constraints in the primal and the dual problems to become linearly dependent, causing the dual problem not to have any feasible solution and the primal problem to yield infinities. A practical (though probably not optimal) solution to this is to use only a subset of all possible A i in the SDP procedure. Note that this only relaxes the program in Eq. (20) (and corresponding program for minimum eigenvalue) and hence increases the range [k min ]. One way to do it is to create an increasing list of random A i matrices and feed it to the SDP until they become linearly dependent. As long as the matrices are not linearly dependent and are of a comparable norm, we expect the solution of SDPs (18) and (20) with P V taken from |Ω p to produce results that are close to the exact case. The 1D numerical tests that we present in the following section support this intuition.

V. NUMERICAL TESTS
Due to the reduced computational cost, we performed all of our tests on 1D systems whose ground states are described by MPS. Nevertheless, we strongly believe that our conclusions hold also for 2D systems, which we leave for future research.
We analyzed four well-known systems, which we label System-A, System-B, System-C, and System-D. All systems are defined over N = 100 spins with open boundary conditions. The underlying Hamiltonians are nearestneighbors with unique ground states and an O(1) spectral gap and correlation lengths which were calculated numerically. The full details are as follows: System-A -1D AKLT. This is a frustration-free system on spin one particles (d = 3) 36,37 given by The ground state is known analytically, and is described by the AKLT valence bond state, which can be written as a MPS with bond dimension D = 2. The spectral gap in the thermodynamic limit is ∆ A 0.35 38 . The correlation length is ξ A 0.9.
System-B -Transverse Ising model. This is the classical Ising model equipped with a transverse magnetic field in theẑ direction: At h = 1.0 the model experiences a phase-transition and the gap closes down. We used h = 1.1 for which the gap is ∆ B 0.07 and the correlation length is ξ B 8.8. System-C -Transverse XY model. This model resembles the transverse Ising model, but with an additional S y i S y i+1 interaction term: We used h = 1.1 and α = 0.5, for which the gap was evaluated numerically to be ∆ C 0.07 and the correlation length ξ C 4.3.
System-D -XY model with random field. Just like System-C, only that here the transverse field at spin i is given by h i , which is a uniformly distributed random number in the range [1.05, 1.15].
For h = 1.1, α = 0.5 and h i ∼ [1.05, 1.15], the gap was evaluated numerically to be ∆ D 0.06. The average correlation length was just like as in System-C. While system A is frustration-free, systems B,C,D are frustrated with substantially smaller spectral gaps and larger correlation lengths. Their ground states were found using a standard DMRG program with a constant bond dimension D. We first verified that the local expectation values do not change by more than 10 −5 when passing from D = 6 to D = 20. This indicated that D = 6 is a good enough bond dimension for our systems. However, in order to suppress as much as possible the finite truncation errors, we first computed the ground state as an MPS with D = 20 and then defined the subspace V L by truncating the bonds at the two cuts from index value 7 onward. All together, this defined a subspace V L of dimension q = 6 × 6 = 36. The essential difference between this procedure and directly using the D = 6 MPS is that now the 36 vectors {|I α } that span V L are resolved inside the ball L using D = 20 instead of D = 6.
In all systems we used 2-spin observables that acted on spin numbers 50,51 in the middle of the chain. For system-A we used three random projectors. Each projector was chosen by randomly picking a 2-dimensional subspace out of the 9-dimensional subspace of two spin-1 particles. In the B,C,D systems we used the three projectors: (i) P x P x projects to the product state of two spin up in thex direction, i.e., P x P x = |+ +| ⊗ |+ +|, (ii) P z P z = |0 0| ⊗ |0 0| is identical to P x P x but in theẑ direction, and (iii) a projector to a random (pure) state in the 4-dimensional space of the two spin 1/2 particles.
For each system and each observable we applied the basic and CGO algorithms with = 3, that corresponds to a ball L of 8 spins, and = 4, with 10 spins. Since System-A is frustration-free we only applied the basic algorithm to it.
In the frustrated systems our SDP calculations were performed using the open-source package SDPA 39,40 . As noted in Sec. IV D, when working with constant D, Eq. (15) is usually over-determined, and therefore one should not use a complete spanning basis of A i matrices. Following the idea from that section, we used a set of random anti-Hermitian matrices A i which were created as random MPO of bond dimension D = 20. Then we ran the SDP for the lower and upper bounds on B , gradually increasing the number of A i .
We used the following heuristics to decide when to stop adding more A i 's and read the final result. First, at every step, we calculated the difference between the upperbound and the lowerbound to B , and this served as a natural error scale. If the difference became negative (i.e., lowerbound was greater than the upperbound), we took the results of the previous step. Additionally, we stopped if the graph started to "wiggle": if at one step the upperbound or lowerbound were worsened with respect to the previous step by more than 1% of the er- ror length scale, we stopped and used the result of the previous step. In all runs (in systems B,C,D) this procedure resulted in an SDP estimate that is based on about 520-710 A i matrices.
A typical evolution of the upper/lower bounds for = 3, 4 as the number of A i increases is presented in Fig. 4. The graph shows the evolution of the bounds in the case of system D (the XY model with random transverse field) with the random projector.
The full numerical results for the four systems are summarized in Tables I,II,III,IV. Table I presents the = 3, 4 basic algorithm results for system A, the frustration-free AKLT model. The other three tables are for the frustrated models B,C,D, and they also include the CGO results. For each test we present the lowerbound and the upperbound in the form of a [k min , k max ] range, as well as the estimated error ∆ := 1 2 (k max − k min ). In the frustrated case, the = 3, 4 radii of the local environment are well below the range where decay of correlation should be felt. It is not surprising then that the basic algorithm in these cases performs very badly, and is unable to give a single meaningful digit in the approximation of B . At the same time, quite remarkably, the CGO algorithm manages to recover 3-4 digits (and often 4-5 digits) from B . Moreover, these results are always better than their equivalent basic results for the frustration-free AKLT model, despite having a much larger correlation length (ξ = 8.8 and ξ = 4.3 vs. ξ = 0.9 for the AKLT) and a much larger V L subspace (q = 36 vs. q = 4 for the AKLT).

VI. CONCLUSIONS
We have presented a new approach to approximate the expectation value of a local operator given a PEPS tensor-network. Our initial observation is that while this is generally a computationally hard task (#P-hard), it is not necessarily the case for PEPS that represent gapped ground states, since these states have much more structure. Our approach circumvents the exponentially expensive contraction of the environment by estimating the local expectation value using only a local patch of the PEPS tensor-network around it. We have presented two algorithms to accomplish that. The basic algorithm simply diagonalizes the observable in the subspace V L and uses its extreme eigenvalues as upper/lower bounds for the expectation value. We have argued that this range is expected to converge as the radius of the local patch increases; it should in particular be useful for frustrationfree systems. The second algorithm, CGO, builds upon the basic algorithm, but uses optimization over commutators to narrow down the eigenvalue range significantly. For it to work, the system has to be frustrated.
As demonstrated by our numerical tests, both algorithms allow one to extract a significant amount of information about the local expectation value using only a local patch of the surrounding tensor network. Contrary to what one might have expected, it seems that in frustrated systems much more information can be extracted from a local patch. In hindsight, this is reasonable: in frustrated systems the ground state is not a common eigenstate of all local terms 41 ; only once we consider all local terms, we satisfy the eigenstate equation. This implies that the action of one local term must somehow cancel the action of other local terms, and this inter-dependence leads to many non-trivial constraints that can be exploited to recover ρ L from its underlying subspace (see Corollary 3).
Our work leaves many open questions for future research. From a numerical point of view, the most immediate question we would like to answer is how effective the CGO algorithm is in 2D. As it stands now, repeating the same procedure that we have used in the 1D case is impractical. Tentatively, for the CGO algorithm to work, one needs a high enough D that would yield good approximation to the groundstate of a frustrated system, and in addition the dimension of the V L subspace should be smaller than the physical dimension of the spins inside L -i.e., the area law should be non-trivially felt. Consider, for example, the rectangular grid in Fig. 1. There, we have dim V L = D 2(4 +3) while the physical dimension goes like d 2( +1) 2 . So if we consider a D = 4 PEPS over d = 2 spins, the physical dimension wins over dim V L at = 7, and reasonable results are expected at = 8 or beyond. This would mean one has to work with matrices of size 2 162 × 2 162 , an impossible task. It is therefore clear that the CGO algorithm cannot be used directly in 2D. Nevertheless, it might still be possible to apply it in approximate manner. For example, by projecting the system to a random subspace and using concentration results like the Johnson-Lindenstrauss lemma 42 , or perhaps, by exploiting the symmetries in the problem (like translation invariance) in a clever way.
It would be also interesting to see if some aspects of CGO can be used in existing PEPS algorithms. The constraints of Corollary 3 may be useful in improving the simple-update method; even if we deal with small regions L that do not provide a unique answer for ρ L , sat-isfying the constraints in Corollary 3 may improve upon the simple product-state approximation of the environment. In addition, they might be useful for increasing the numerical stability of current algorithms. Finally, from a numerical point of view it would be interesting to see if Corollary 3 could be used to verify that we have an eigenstate at our hands. This might be useful, for example, for settling the nature of the groundstate of the antiferromagnetic Heisenberg model on the kagome lattice [43][44][45][46][47] .
From a more theoretical point of view, it would be interesting to gain better understanding of the CGO algorithm, and in particular the implications of the constraints in Corollary 3. For example, can we find sufficient conditions to when the algorithm gives good approximations (i.e., narrow [k min , k max ] range)? From a computational complexity point of view, this might be a step in showing that PEPS of ground states can serve as an NP witness. A much more ambitious goal in that direction would be to show that the complexity of the local Hamiltonian problem for gapped Hamiltonians on a regular lattice with a unique ground state is inside NP.