Tangent-space methods for truncating uniform MPS

A central primitive in quantum tensor network simulations is the problem of approximating a matrix product state with one of a lower bond dimension. This problem forms the central bottleneck in algorithms for time evolution and for contracting projected entangled pair states. We formulate a tangent-space based variational algorithm to achieve this for uniform (infinite) matrix product states. The algorithm exhibits a favourable scaling of the computational cost, and we demonstrate its usefulness by several examples involving the multiplication of a matrix product state with a matrix product operator.

The density matrix renormalization group (DMRG) 1,2 and quantum tensor networks 3,4 provide algorithms for simulating ground states of strongly correlated quantum many body systems with a computational cost that scales linear in the system size, thereby overcoming the infamous exponential wall of the quantum many body problem. The physical parameter controlling the computational cost is the entanglement entropy, as directly reflected in the bond dimension χ of the corresponding matrix product states (MPS) 5 . However, there are many interesting physical problems for which this bond dimension can become prohibitively large, such as the problem of simulating time evolution of a quantum state out of equilibrium or of contracting a tensor network comprised of a projected entangled pair state (PEPS) with a large bond dimension. In both cases, the central problem is to approximate the product of an MPS and a matrix product operator (MPO) with a MPS of lower bond dimension. For both finite and infinite systems, a well known algorithm to achieve this is time-evolving-blockdecimation and variants thereof [6][7][8][9] . For finite systems, a considerable improvement over those algorithms can be obtained by adopting a variational algorithm which optimizes the fidelity by sweeping through the system while solving alternating linear problems 10,11 . The computational cost of the latter algorithm has a better scaling as it does not require to bring the joint MPS/MPO system in canonical form, and furthermore achieves a better overal fidelity due to its variational nature.
In this paper, we present the uniform and infinite version of that algorithm. It is based on ideas developed in the context of tangent space methods for uniform matrix product states 12,13 and the variational uniform matrix product state algorithm 14,15 . Our main motivation is the development of efficient MPS algorithms which can deal with time-evolution methods involving MPOs with large bond dimension and of efficient and wellconditioned ways of contracting PEPS 16 . It also overcomes a main limitation of algorithms based on the timedependent variational principle (TDVP) [17][18][19] , where it is difficult to build up entanglement starting from a lowentangled state by allowing large time steps.
The paper is organized as follows. In the first section we discuss how to approximate a given uniform MPS variationally with another with with lower bond dimension.
In a second section, we illustrate this algorithm with several relevant examples. Fixed-point equations.-We start from the diagrammatic expression of a uniform MPS in the thermodynamic limit, parametrized by a single tensor A (1) We will assume a trivial unit cell in this text for simplicity, the case of larger unit cells is treated straightforwardly. Using the gauge freedom of the MPS we can choose this tensor to be in the left canonical gauge A L or the right canonical gauge A R , with These gauge-fixed tensors are related by a matrix C as allowing us to bring the MPS into the so-called mixed gauge For a given MPS |Ψ(M described by a tensor M , we now wish to find an MPS |Ψ(A) such that the latter approximates the former in some optimal way. A natural choice for an optimality condition is that they should have a maximal fidelity, which leads us to a variational optimization problem for the tensor A A = arg max This cost function being a real-valued function of A and A * , the gradient is obtained by differentiating the cost function with respect to A * . An optimal point is reached when the gradient vanishes, The left-hand side of this equation can be interpreted as a tangent vector on the manifold of MPS 12,13 , and the optimality condition can be reformulated as where P A represents the projector on the space of tangent vectors to |ψ A . An explicit form of the tangent-space projector in the mixed-gauge is given by 13 Applying this operator to |ψ M , which we assume to be a uniform MPS parameterized by a single tensor M , we find that the optimality condition [Eq. (7)] implies that where A C and C are given by and with the fixed points G L and G R given by the eigenvalue equations Here, the factor λ appears as the 'fidelity per site', formally given in the thermodynamic limit as In order to find a fixed point, we can use an iterative scheme. One crucial step in each iteration will be the compute error Eq. (16) 7: until < η 8: return AL, AR, λ extraction of a new set of MPS tensors {A L , A R } from the A C and C that were obtained. A close-to-optimal solution of this problem is given by the prescription 13 and where all decompositions involve unique polar decompositions or their transposed. This approach is very similar to the one adopted in the standard variational uniform MPS (VUMPS) algorithm 15 . Once we have obtained a new set {A L , A R }, we can re-compute the fixed-point tensors G L and G R and the scheme can be reiterated. As a convergence measure we take the norm of the fixed-point equation in Eq. (7), which is given by where A C and C are given by Eqs. 10 and 11. A specific instance of the above scheme occurs when applying a uniform matrix product operator (MPO) to a given MPS, and approximating the resulting state as an MPS with a certain bond dimension. In that case the above fixed-point equations are given by and Our variational method can, therefore, be used for approximating an MPS-MPO state by an MPS with the original bond dimension of M . This is an operation that appears in many MPS methods (see further), and we can show that our approach scales more favourably as compared to the standard local-truncation approach. Indeed, supposing that both the original and new MPS have bond dimension χ and physical dimension d and the MPO has bond dimension D, the time-complexity of the above scheme is O(χ 3 Dd+χ 2 D 2 d 2 ), and the memory required scales as O(χ 2 Dd). We can compare this to the complexity of cutting the bond dimension by truncating local Schmidt values. The most costly operation required to cut the bond this way is following contraction: The time-complexity of this operation is O(χ 3 D 2 d + χ 2 D 3 d) and the memory required O(χ 2 dD 2 ). In addition, one typically performs a full singular-value decomposition of a square χD matrix, for which the time complexity scales as O(χ 3 D 3 ). This analysis shows that for MPOs of large virtual dimension D, the method we prescribe can be a significant, even crucial, improvement.
Truncating an MPS.-Let us first illustrate the method by truncating the bond dimension of a given MPS. The most commonly used technique for that purpose is the truncataction of the Schmidt values on all bonds simultaneously 7 . We compare the two techniques in Fig. 1 for an MPS of considerable dimension. We find that truncating all Schmidt values simultaneously performs fairly well across the board, but that our variational scheme still finds a slightly better state after convergence. This example shows that our fidelity optimization can be useful only if precision is of the utmost importance. Time evolution.-There are roughly two different classes of methods used to time-evolve an infinite MPS. The first class tries to directly transform the Schrödinger equation into a (non-linear) differential equation on the variational manifold. This is exactly the mechanism behind TDVP 17,19 , where the direction in which the state needs to change (the right hand side of the Schrödinger equation) is projected onto the tangent space of the MPS. The second class of methods instead starts from an approximation of the time evolution operator exp(−iHδ) for a certain time step δ. This approximation is provided in terms of a quantum circuit, or, more generally, an MPO, and can be obtained from e.g. a Suzuki-Trotter decomposition 6,20,21 or a cluster expansion 16,22 . The resulting MPO is then applied to the current state, encoded as MPS, followed by a bond truncation 1 . With a (low-order) Suzuki-Trotter decomposition, the MPO bond dimension can remain low, but feasible time steps δ are also very small. With the cluster expansion, it is easier to reach larger δ, at the cost of a higher MPO bond dimension. It is therefore infeasible to apply this MPO to an MPS and truncate directly according to the Schmidt values due to prohibitive memory constraints or time complexity considerations. In this case thus, our method is indispensable. We illustrate this usage by evolving the Néel state with the XXZ Hamiltonian., where S α i the spin-1/2 operators at site i and we choose ∆ = 1/2. This problem is closely related to the one considered in Ref. 23 asserting the supremacy of quantum simulators. We have exploited the U (1) symmetry of the system and used an MPS bond dimension of 994. The MPO bond dimension is 21, which enabled an accurate time step of up to dt = 1.2. In Fig. 2 we show the occupation number as a function of time, and benchmark it with a simulation with the TDVP algorithm with much smaller time steps. Power method for transfer matrices.-Let us now consider the calculation of an MPS fixed point of an MPO transfer matrix by way of the power method: In each iteration we apply the MPO and truncate the bond dimension, until the MPS converges to a fixed point. Power methods have been used for computing transfer fixed points where the local singular-value truncation was adopted in each iteration 9 , but here we use our variational truncation. In contrast to the former, the fixed point of our variational-truncation approach is, in fact, a variationally optimal MPS in the sense that it optimizes the leading eigenvalue for hermitian transfer matrices. Indeed, in the fixed point of this power method, the toplayer MPS in the fixed-point equations [Eqs. (17)- (19)] should be the same as the down-layer, and the equations reduce to the usual fixed-point equations of the VUMPS algorithm (which is variationally optimal for hermitian transfer matrices). Hence, both approaches share at least the same fixed point, which is not true if truncation based on singular values is performed.
For hermitian transfer matrices the performance of a power method is inferior to that of the Krylov-inspired VUMPS algorithm, but it is very useful in cases of spatial symmetry breaking where the fixed point alternates between different MPSs or for non-Hermitean MPOs. We illustrate this case by studying the MPO transfer matrix of the classical antiferromagnetic Ising model on the square lattice (Fig. 3). In the (low-temperature) symmetry-broken phase, we find that the power method alternates between two MPSs that are the same up to a one-site translation. We look at some convergence criteria and also compare to the sublattice rotated ferromagnetic fixed point, found using VUMPS. Additionally we find that this technique allows for slightly better convergence of the fixed point MPS than VUMPS, as can be be seen from the stagnation of the fourth (purple) curve and the continued convergence of the second (red) curve and in Fig. 3. Dynamical growing of bond dimension.-Our variationaltruncation approach is particularly useful as a way of enlarging the bond dimension of an MPS when simulating time evolution or computing fixed points of transfer matrices. With respect to the former, the most persistent critique to the TDVP algorithm revolves around the fact that it projects the time evolution on the manifold of MPS with a fixed bond dimension, and that it is impossible to grow the bond dimension during the evolution. Our variational algorithm is not confined to a manifold of fixed bond dimension, because we can choose the bond dimension at each time step. We believe that a 'hybrid' between TDVP and our current scheme can provide a good way of simulating time evolution variationally using MPS where the amount of entanglement increases through time.
For fixed points of transfer matrices we can exploit our fidelity optimization in a similar way. We imagine the situation in which we have found a fixed-point MPS of a certain bond dimension, and we wish to find a better MPS of larger bond dimension. We can now use the previous MPS to construct an initial guess, apply the transfer matrix to this MPS, and then truncate to an MPS of the desired bond dimension using the equations above [Eqs. (17)- (19)]. The resulting MPS is already a more accurate approximation of the desired state than the previous one, and thus makes an excellent initial guess for running a new fixed-point algorithm at this higher bond dimension. This is especially useful in the context of PEPS algorithms, where the fixed point calculation of the PEPS double layer is the main bottleneck.
Conclusions.-We have discussed a method for approximating a uniform and infinite MPS by an MPS of lower bond dimension in a way that is variationally optimal. This method is proven most useful if the MPS being approximated has some substructure, like being made up of an MPO times and MPS. In this case the method has lower complexity and requires less memory than stan-dard alternatives. We illustrate this with time evolution using an MPO that approximates the evolution operator, a power method for finding transfer matrix fixed points, and dynamical growing of bond dimension.
The generalization of this method to the (2+1)dimensional case can easily be envisioned, and would be interesting to investigate. An algorithm that variationally determines a PEPS approximation of some other PEPS-perhaps a projected entangled-pair operator (PEPO) times a PEPS-can readily be devised based on the algorithm in Ref. 24. The uses of such a method would be identical to the ones presented here: performing accurate and reliable time evolution, a power method for determining fixed points of non-hermitian PEPOs or PE-POs exhibiting spatial symmetry breaking, and growing of a PEPS bond dimension.