TL;DR
This paper introduces VUMPS, a variational algorithm combining DMRG and tangent space concepts, which improves convergence and precision for ground state calculations of 1D quantum systems, including long-range and 2D models.
Contribution
The paper presents VUMPS, a novel variational uniform MPS algorithm that outperforms existing methods like IDMRG and ITEBD in speed and accuracy for various quantum lattice models.
Findings
VUMPS converges faster than IDMRG and ITEBD.
VUMPS achieves higher precision in ground state energies.
VUMPS efficiently handles long-range interactions and 2D models.
Abstract
We combine the Density Matrix Renormalization Group (DMRG) with Matrix Product State tangent space concepts to construct a variational algorithm for finding ground states of one dimensional quantum lattices in the thermodynamic limit. A careful comparison of this variational uniform Matrix Product State algorithm (VUMPS) with infinite Density Matrix Renormalization Group (IDMRG) and with infinite Time Evolving Block Decimation (ITEBD) reveals substantial gains in convergence speed and precision. We also demonstrate that VUMPS works very efficiently for Hamiltonians with long range interactions and also for the simulation of two dimensional models on infinite cylinders. The new algorithm can be conveniently implemented as an extension of an already existing DMRG implementation.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16| W=6 | D=256 | D=512 |
|---|---|---|
| VUMPS | -0.672544677277 | -0.672724840927 |
| IDMRG | -0.672518318233 | -0.672724792693 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Variational optimization algorithms for uniform matrix product states
V. Zauner-Stauber
Vienna Center for Quantum Technology, University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria
L. Vanderstraeten
Ghent University, Faculty of Physics, Krijgslaan 281, 9000 Gent, Belgium
M.T. Fishman
Institute for Quantum Information and Matter, California Institute of Technology, Pasadena, California 91125, USA
F. Verstraete
Vienna Center for Quantum Technology, University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria
Ghent University, Faculty of Physics, Krijgslaan 281, 9000 Gent, Belgium
J. Haegeman
Ghent University, Faculty of Physics, Krijgslaan 281, 9000 Gent, Belgium
Abstract
We combine the Density Matrix Renormalization Group (DMRG) with Matrix Product State tangent space concepts to construct a variational algorithm for finding ground states of one dimensional quantum lattices in the thermodynamic limit. A careful comparison of this variational uniform Matrix Product State algorithm (VUMPS) with infinite Density Matrix Renormalization Group (IDMRG) and with infinite Time Evolving Block Decimation (ITEBD) reveals substantial gains in convergence speed and precision. We also demonstrate that VUMPS works very efficiently for Hamiltonians with long range interactions and also for the simulation of two dimensional models on infinite cylinders. The new algorithm can be conveniently implemented as an extension of an already existing DMRG implementation.
I Introduction
The strategy of renormalization group (RG) techniques to successively reduce a large number of microscopic degrees of freedom to a smaller set of effective degrees of freedom has led to powerful numerical and analytical methods to probe and understand the effective macroscopic behavior of both classical and quantum many body systems.Wilson (1975); Fisher (1974, 1998); Zinn-Justin (1996) However, it was not until the advent of White’s celebrated Density Matrix Renormalization Group (DMRG) White (1992, 1993) that variational RG methods reached unprecedented accuracy in numerically studying strongly correlated one dimensional quantum lattice systems at low temperature. The underlying variational ansatz of Matrix Product States (MPS) Fannes et al. (1992); Klümper et al. (1993); Östlund and Rommer (1995); Perez-Garcia et al. (2007); Verstraete et al. (2008); Schollwöck (2011); Haegeman et al. (2013a) belongs to a class of ansatzes known as Tensor Network States.Verstraete et al. (2008); Orús (2014); Bridgeman and Chubb (2016) These variational classes encode the many body wavefunction in terms of virtual entanglement degrees of freedom living on the boundary and thus satisfy an area law scaling of entanglement entropy per construction. As such, they provide a natural parameterization for the physical corner of Hilbert space, where low energy states of quantum many body systems ought to live in.Hastings (2007a); Eisert et al. (2010) MPS in particular are especially fit for studying ground states of strongly correlated one dimensional quantum systems with local interactions.Verstraete and Cirac (2006); Hastings (2007b); Haegeman et al. (2013b)
The variational parameters in MPS are contained within local tensors associated with the individual sites of the lattice system. For homogeneous systems, the global wave function can then be captured using just a single (or a small number of) such tensors, independent of the system size. They consequently offer very natural access to the thermodynamic limit, providing a clear advantage over other numerical approaches such as Exact Diagonalization or Quantum Monte Carlo.
On finite lattices, (one-site) DMRG implements the variational principle (energy minimization) by exploiting that the quantum state is a multilinear function of the local tensors. By fixing all but one tensors, the global eigenvalue problem is transformed into an effective eigenvalue problem for the local tensor.White (1992, 1993); Verstraete et al. (2004a); White (2005); McCulloch (2007); Schollwöck (2011) Using a translation invariant parameterization gives rise to an energy expectation value with a highly non-linear dependence on the tensor(s). Two different algorithms are widely used to obtain such an MPS in the thermodynamic limit. Infinite system DMRG (IDMRG)White (1992, 1993); McCulloch (2008) proceeds by performing regular DMRG on a successively growing lattice, inserting and optimizing over new tensors in the center of the lattice in each step only, effectively mimicking an infinite lattice by using a finite, albeit very large lattice. After convergence the most recently inserted tensors in the center are taken as a unit cell for an infinite MPS approximation of the ground state. An alternative approach is known as infinite time evolving block decimation (ITEBD).Vidal (2003, 2007) It works directly in the thermodynamic limit and is based on evolving an initial state in imaginary time by using a Trotter decomposition of the evolution operator.
We present a new variational algorithm, inspired by tangent space ideas,Haegeman et al. (2011, 2013a, 2016) that combines the advantages of IDMRG and ITEBD and addresses some of their shortcomings. As such it is directly formulated in the thermodynamic limit, but at the same time optimizes the state by solving effective eigenvalue problems, rather than employing imaginary time evolution. We find that it leads to a significant increase in efficiency in all of our test cases. The following section introduces MPS notations and definitions and presents our variational algorithm, heuristically motivated from the perspective of finite size DMRG. Sec. III illustrates the performance of our algorithm on various test cases, and compares to conventional IDMRG and ITEBD results. After the conclusion in Sec. IV, we provide further technical details in the appendices. Appendix A contains additional theoretical background: we derive the self-consistent conditions that characterize the variational minimum and provide additional motivation for our algorithm from the perspective of the MPS tangent space. Appendix B presents a suitable strategy to expand the bond dimension of translation invariant MPS. Appendix C explains how to construct effective Hamiltonians in the thermodynamic limit. These involve infinite geometric sums of the transfer matrix, which are further studied in Appendix D.
II A Variational Algorithm for Matrix Product States in the Thermodynamic Limit
In this section we introduce a variational algorithm for optimizing MPS directly in the thermodynamic limit. Because the algorithm strongly resembles conventional DMRG, we explain it by describing a single iteration step from the viewpoint of DMRG and show that only a few additional ingredients are needed to arrive at our variational algorithm. We only briefly motivate these extra ingredients for the sake of readability, and refer to Appendix A for additional explanations and more rigorous theoretical motivations. As such, the new algorithm can easily be implemented as an extension to an already existing (I)DMRG implementation.
We start by considering a setting familiar from conventional DMRG: a finite homogeneous one dimensional quantum lattice system, where every site corresponds to a level system. We label the sites by an integer and thus have a basis for the local Hilbert space on site . The total Hilbert space is spanned by the product basis . We assume the dynamics of the system to be governed by a translation invariant Hamiltonian .
We further consider a variational parameterization of a ground state approximation of the system, for now in terms of a finite size (site dependent) MPS, but we will ultimately be interested in the thermodynamic limit. DMRG proceeds to find the best variational ground state approximation by employing an alternating least squares minimization: It starts from some initial state and successively optimizes each of the individual MPS tensors site by site by solving effective (Hamiltonian) eigenvalue problems, in a sweeping process through the lattice until convergence, where each iteration depends on already optimized tensors from previous iterations (see e.g. Refs White, 1992, 1993; Verstraete et al., 2004a; McCulloch, 2007; Schollwöck, 2011).
We are now however interested in the thermodynamic limit (but will ignore the technical complications involving a rigorous definition of a Hilbert space in that limit). In that case the MPS ground state approximation will be given in terms of a translation invariant uniform MPS, described by a single MPS tensor (or a unit cell of tensors), repeated on all sites. Two immediate difficulties arise: Firstly, conventional DMRG updates the variational state site by site, thus breaking translation invariance. Secondly, the effective Hamiltonian for a single-site optimization has to be constructed from an infinite environment.
After briefly introducing the variational class of uniform MPS and introducing necessary notation and conventions (for further details see Sec. A.2), we describe how the new algorithm modifies DMRG accordingly to exactly account for these two issues in order to arrive at a variational ground state algorithm directly formulated in the thermodynamic limit.
II.1 Uniform MPS
A uniform MPS (uMPS) of bond dimension defined on an infinite translation invariant lattice is parameterized by a single collection of matrices for . The overall translation invariant variational state is then given by
[TABLE]
and can be represented diagrammatically as
[TABLE]
Exploiting the invariance of (1) under local gauge transformations , with invertible, the state can be cast into certain favorable representations, among them the left and right canonical representation
[TABLE]
or diagrammatically
[TABLE]
Here and correspond to the left and right reduced density matrices of a bipartition of the state respectively. We henceforth refer to () as a left (right) isometric tensor, or just a left (right) isometry.
Defining the left and right transfer matrices
[TABLE]
and using the notation and for vectorizations of a matrix in the dimensional “double layer” virtual space the transfer matrices act upon, the gauge conditions (2) are equivalent to
[TABLE]
i.e. and are the left and right dominant eigenvectors of , while and are the left and right dominant eigenvectors of .
As is standard practice in DMRG, we mix both of these representations and cast the state into the mixed canonical representation
[TABLE]
or diagrammatically,
[TABLE]
Here we have defined the center site tensor (known as the single-site wave function in DMRG)
[TABLE]
in terms of the bond matrix , which constitutes the (invertible) gauge transformation relating and via . The singular values of then encode the entanglement spectrum of the state. Indeed, using we can verify that the left and right reduced density matrices in (2) are given by and . Furthermore, ensures that (5a) and (5b) are translation invariant and that and can be shifted around arbitrarily. Normalization of the state, as well as of the reduced density matrices and , corresponds to the single condition .
For ease of notation we further introduce the following partial states
[TABLE]
with arbitrary, and use them to define the reduced basis states
[TABLE]
II.2 Effective Hamiltonian
The use of the mixed canonical representation (5a) in DMRG is of significant importance for the stability, as it reduces the minimization of the (global) energy expectation value with respect to into a standard (hermitian) eigenvalue problem, instead of a generalized one. The effective Hamiltonian for this eigenvalue problem is the system Hamiltonian projected onto the degrees of freedom of , and is known as the “reduced” or “superblock” Hamiltonian in DMRG.
We define the thermodynamic limit version of this reduced single-site Hamiltonian acting on as
[TABLE]
Additionally, we also define an effective Hamiltonian acting on the bond matrix as
[TABLE]
which does not appear directly in the context of DMRG, but will be needed later for a consistent update of the state without breaking translation invariance. It can be interpreted as a “zero site” effective Hamiltonian, which would feature in an optimization of the global energy expectation value with respect to the Schmidt values.
In an efficient implementation, these effective eigenvalue problems are typically solved using an iterative eigensolver, so that we only need to implement the action of and onto and .
While the energy expectation value is extensive and thus divergent in the thermodynamic limit, the effective Hamiltonians and are well defined and finite in the thermodynamic limit if one properly subtracts the current energy expectation value from the Hamiltonian . We demonstrate this procedure for the case of nearest neighbor interactions , where the two-site Hamiltonian acts on neighboring sites only. We refer to Appendix C for the case of long range interactions and for general Hamiltonians given in terms of Matrix Product Operators (MPOs).
In the case of nearest neighbor interactions, the action of onto splits up into four individual contributions, which follow from the decomposition (left block containing sites , center site , and right block containing sites ). The action of onto is given by
[TABLE]
where the first two terms correspond to the Hamiltonian terms and coupling the center site to the left and right block, respectively, and and sum up the contributions of all the Hamiltonian terms acting strictly to the left and to the right of the center site.
The environments and are usually constructed iteratively while sweeping through the (finite) lattice in conventional DMRG, or grown successively in every iteration in IDMRG. In the thermodynamic limit, these terms consist of a diverging number of individual local interaction contributions , and care needs to be taken in their construction.
Indeed, the contribution to comes from the Hamiltonian term and is given by . Likewise, is the contribution to stemming from . Here, we have used the definitions
[TABLE]
or diagrammatically
[TABLE]
Summing up all such local contributions gives rise to infinite geometric sums of the transfer matrices
[TABLE]
where can be presented diagrammatically as
[TABLE]
and likewise for .
The transfer matrix has a dominant eigenvalue of magnitude one, with corresponding left and right eigenvectors and . The projection is the energy density expectation value and is independent of . Subtracting the energy from the Hamiltonian as , we can write . The second term is exactly proportional to the left eigenvector of eigenvalue and therefore gives rise to a diverging contribution in the geometric sum, corresponding to the total energy of the left half infinite block. Since this contribution acts as the identity in the effective Hamiltonian [Eq. (11)], we can however safely discard this diverging contribution without changing the eigenvectors of . This corresponds to an overall energy shift of the left half infinite block such that . For the remaining part the geometric sum converges. With the same comments apply to the construction of .
We can evaluate and recursively as
[TABLE]
with initialization and . We can repeat these recursions until e.g. drops below some desired accuracy . This strategy is conceptually simple and closely resembles the successive construction of the environments in the context of (I)DMRG, but is not very efficient, as its performance is comparable to that of a power method.
A more efficient approach is to formally perform the geometric sums in (13) explicitly, and to iteratively solve the resulting two systems of equations
[TABLE]
for and to precision , as explained in Appendix D.
So far, we have discussed the action of . The action of onto follows simply from (11) by projecting onto or . Using the defining property of or , the result simplifies to
[TABLE]
The first two terms of (11) can be applied in operations111In many cases is sparse and the number of non-zero elements is usually of the order . The first two terms can then be applied in operations., and the last two terms in operations. For (16) the first term can be applied in operations, and the last two terms in operations. To generate the necessary terms for (11) and (16) we have to iteratively evaluate two infinite geometric sums, involving operations (when iteratively solving (15) the solutions from the previous iteration can be used as starting vectors to speed up convergence). A pseudocode summary for obtaining the necessary explicit terms of and and their applications onto a state is presented in Table 1.
II.3 Updating the state
In DMRG, we would update the state by replacing with the lowest eigenvector of and then shift the center site to the right by computing an orthogonal factorization , or to the left by computing . As such, the state gets updated by only replacing the current site with or , while leaving all other sites untouched. However, applying this scheme in our setting would immediately destroy translation invariance after a single step.
We want to construct an alternative scheme that applies global updates in order to preserve translation invariance at any time. Global updates can most easily be applied with an explicit uniform parameterization in terms of a single tensor . On the other hand, DMRG experience teaches us that the stability is greatly enhanced when applying updates at the level of and , which are isometrically related to the full state.
We therefore calculate the lowest eigenvector of like in DMRG, but additionally also the lowest eigenvector of . We then globally update the state by finding new and as the left and right isometric tensors that minimize and respectively. These minimization problems can be solved directly (not iteratively) and without inverting (see below). As shown in Appendix A, at the variational optimum the values of these objective functions go to zero, and current and are the lowest eigenvectors of and respectively.
For the remainder of this section we omit tildes and use the following matricizations of the 3-index tensors
[TABLE]
We thus want to extract updated and from updated and by solving
[TABLE]
In exact arithmetic, the solution of these minimization problems is known, namely will be the isometry in the polar decomposition of (and similar for , see Thm. IX.7.2 in Ref. Bhatia, 1997). Computing the singular value decompositions (SVD)
[TABLE]
we thus obtain
[TABLE]
Notice that close to (or at) an exact solution , the singular values contained in are the square of the singular values of , and might well fall below machine precision. Consequently, in finite precision arithmetic, corresponding singular vectors will not be accurately computed.
An alternative that has proven to be robust and still close to optimal is given by directly using the following left and right polar decompositions
[TABLE]
to obtain
[TABLE]
where matrices are hermitian and positive. Alternative isometric decompositions might be considered in Eq. (21), though it is important that they are unique (e.g. QR with positive diagonal in ) in order to have close to convergence.
II.4 The Algorithm: VUMPS
We are now ready to formulate our variational uniform MPS (VUMPS) algorithm. As shown in Appendix A, a variational minimum (vanishing energy gradient) in the manifold of uMPS is characterized by tensors , and satisfying the conditions
[TABLE]
Here, bold symbols denote vectorizations of the MPS tensors and matricizations of the effective Hamiltonians, and and are the lowest eigenvalues of the effective Hamiltonians.222These values can be different and depend on the subtraction scheme for the divergent energy expectation value. If is performed everywhere, we have . If, in the case of nearest neighbor interactions, we only substitute in the construction of and , but not in the local terms, we will have .
When iterating the steps outlined in the previous subsections, convergence is obtained when these conditions are satisfied. In particular, starting with a properly orthogonalized initial trial state of some bond dimension , we begin by solving the two eigenvalue problems for the effective Hamiltonians and . Since we are still far from the fixed point, the resulting lowest energy states and will in general not satisfy the gauge condition (23c) together with current .
Following the procedure of the previous section we can however find optimal approximations and for (23c) to arrive at an updated uMPS. Conversely, and will not be the correct lowest energy eigenstates of the new effective Hamiltonians and generated from . We then use the updated state and reiterate this process of alternately solving the effective eigenvalue problems, and finding optimal approximations for and to update the state. For a pseudocode summary of this algorithm, see Table 2.
We now elaborate on the various steps in the VUMPS algorithm. Firstly, extracting new from updated and can be done using the theoretically optimal (but numerically often inaccurate) Eq. (20) or the more robust Eq. (22), depending on the magnitude of the smallest singular value in . As a good uMPS approximation will always involve small singular values, Eq. (22) is preferable most of the time, except maybe during the first few iterations.
The maximum of the error quantities (18)
[TABLE]
provides an error measure for the fixed point condition in Eq. (23c) and is used as a global convergence criterion. It measures the precision of the current uMPS ground state approximation. Within every iteration, we use iterative methods (e.g. some variation of Lanczos) to find the eigenvectors and of the Hermitian operators and . As the goal is to drive the state towards the fixed point relations in Eqs. (23a) and (23b), it is not necessary to solve these eigenvalue problems to full machine precision. Rather, it is sufficient to use a tolerance chosen relative to . 333Further approximations to comparable accuracy can be made within the construction of the effective Hamiltonians, e.g. when determining and to precision . There, the approximations and for the true and needed for some of these operations are good enough, if is roughly of the same order of magnitude as . A value of of the order of has proven to work well in practice. It is also worthwhile to use tensors from the previous iteration as initial guess for the iterative solvers to speed up convergence.
As the main part of the algorithm works at fixed bond dimension (i.e. it is a single-site scheme in DMRG terminology), one might choose to increase the bond dimension before starting a new iteration. We have developed a subspace expansion technique that works directly in the thermodynamic limit and is explained in Appendix B.
While the true comparison of this algorithm with IDMRG White (1992); McCulloch (2008) and ITEBD Vidal (2007) will take place in Sec. III by gathering actual numerical simulation results, we can already compare the theoretical properties of these algorithms. Neither IDMRG or ITEBD is truly solving the variational problem in the sense of directly trying to satisfy the fixed point conditions Eqs. (23). IDMRG closely resembles regular DMRG on a successively growing lattice, as it inserts and optimizes over new tensors in the center of the lattice in each step. Tensors from previous steps are not updated, as this would render the cost prohibitive. When this approach converges, the resulting fixed point tensors in the center can be assumed to specify the unit cell of an infinite MPS. VUMPS has the immediate advantage that i) it directly works in the thermodynamic limit at all iterations and ii) it completely replaces the entire state after every iteration, thus moving faster through the variational manifold. In contrast, IDMRG keeps memory of earlier iterations and cannot guarantee a monotonically decreasing energy that converges to an optimum associated with a translation invariant MPS in which the effects of the boundary have completely disappeared. The advantages of VUMPS come with a greater computational cost per iteration, as two eigenvalue problems (for and for ) and – in the case of nearest neighbor interactions – two linear systems (for and ) have to be solved. IDMRG only solves a single eigenvalue problem and builds and step by step in every iteration. The latter approach is analogous to a power method for eigenvalue problems and, while very cheap, is expected to require many iteration steps to converge, especially for systems with large correlation lengths (e.g. close to criticality).
ITEBD Vidal (2007) is based on evolving an initial state in imaginary time by using a Trotter decomposition of the evolution operator. Like VUMPS, ITEBD works in the thermodynamic limit at any intermediate step, typically with a unit cell that depends on how the Hamiltonian was split into local terms in order to apply the Trotter decomposition. Furthermore, as every application of the evolution operator increases the virtual dimension of the MPS, truncation steps are required to restore the original (or any suitable) value of the bond dimension. While VUMPS can take big steps through the variational space, time steps in ITEBD have to be chosen sufficiently small (especially in the final steps of the algorithm) to eliminate the Trotter error, which negatively affects the rate of convergence (Ref. Phien et al., 2015 however proposes a scheme to effectively obtain a larger time step). Furthermore, the Trotter splitting essentially limits the applicability of ITEBD to short-range interactions and dictates the size of the unit cell of the resulting MPS, e.g. in the most common case of nearest neighbor interactions a two-site unit cell is obtained. (The approach of Ref. Pirvu et al., 2010 to obtain a translation invariant MPS is restricted to certain Hamiltonians, but see Ref. Zaletel et al., 2015 for an alternative proposal that can in fact also deal with long range interactions.)
Finally, we can also compare VUMPS to the more recent time dependent variational principle (TDVP),Haegeman et al. (2011) which was implemented as an alternative approach to simulate real and imaginary time evolution within the manifold of MPS by projecting the evolution direction onto the MPS tangent space. This approach can be applied to translation invariant MPS, independent of the type of Hamiltonian. When used to evolve in imaginary time, it can be identified as a covariant formulation of a gradient descent method, in that it evolves the state in the direction of the gradient of the energy functional, preconditioned with the metric of the manifold. As such, the energy decreases monotonically and at convergence, an exact (local) minimum is obtained, as characterized by the vanishing gradient. However, in its original formulation, TDVP was not formulated in a center site form and was therefore unstable and restricted to small time steps. For finite systems, a different formulation of the TDVP algorithm was provided in Ref. Haegeman et al., 2016, which allows for taking the limit of the imaginary time step to infinite, and then becomes provably equivalent to the single-site DMRG algorithm. VUMPS can be motivated from these developments, as explained in Appendix A.
We conclude this section by elaborating on how to incorporate symmetries in the algorithm. The construction of uMPS that is explicitly invariant under onsite unitary symmetries is equivalent to (I)DMRG McCulloch (2007); Schollwöck (2011) and (I)TEBD Hastings (2011); Singh et al. (2011), and it is immediately clear that the various steps in VUMPS have a corresponding covariant formulation. The same comments apply to time reversal symmetry, in which case everything can be implemented in real arithmetic, or to reflection symmetry, in which case and will be symmetric matrices and (which implies that and are also related). In all of these cases, the computational cost is reduced. However, explicitly imposing the symmetry in the MPS requires caution, as the physical system might have spontaneous symmetry breaking, or – more subtly – might be in a symmetry protected topological phase where the symmetries cannot be represented trivially on the MPS tensor.
In the case of spontaneous symmetry breaking, MPS algorithms tend to converge to maximally symmetry broken states for which the entanglement is minimal. This is also the case for VUMPS. One can control which state the algorithm converges to by suitably biasing the initial state or by adding small perturbation terms to the Hamiltonian which explicitly break the symmetry, and which are switched off after a few iterations.
Explicit conservation of translation symmetry was the very first requirement in the construction of VUMPS. In the case of spontaneous breaking of translation symmetry down to -site translation symmetry (as e.g. in the case of a state with antiferromagnetic order), enforcing one-site translation symmetry would result in a (non-injective) equal weight superposition of all symmetry broken uMPS ground state approximations. In order to reach an optimal accuracy with a given bond dimension, such a superposition of states is however undesirable, as the effective bond dimension is reduced to . In the case where this situation cannot be amended by a simple unitary transformation that restores one-site translation symmetry (such as e.g. flipping every second spin in the case of an antiferromagnet), it is preferable to choose an MPS ansatz with a -site unit cell, such that the state can spontaneously break translation symmetry. The generalization of the algorithm to multi-site unit cells is described in the next section.
II.5 Multi Site Unit Cell Implementations
We now generalize the VUMPS algorithm of the previous section for one-site translation invariant uMPS to the setting of translation invariance over sites. Such a uMPS ansatz is then parameterized by independent tensors , , which define the unit cell tensor
[TABLE]
where is a combined index. We can then write the variational state as
[TABLE]
and the left and right orthonormal forms are given by the relations
[TABLE]
where it is understood that and .
Defining the bond matrices as the gauge transformation that relates left and right canonical form via , we have and . We can then also cast in a mixed canonical form similar to (5a) with the center site tensor given by .
The variational minimum within this set of states is characterized by the following fixed point relations
[TABLE]
Notice that due to (27c), the relations for different are connected. There are several possible strategies for constructing algorithms which obtain states satisfying these conditions.
In the following we present two approaches which have shown good performance and stable convergence, which we shall term the “sequential” and “parallel” methods. But let us first elaborate on computing effective Hamiltonians for multi-site unit cells, which works similarly in both methods. We again restrict to the case of nearest neighbor interactions, such that the effective Hamiltonians are constructed similar as in Sec. II.2. To construct e.g. the left block Hamiltonian , we first collect all local contributions from a single unit cell in , before performing the geometric series of the transfer matrix, which now mediates a translation over an entire unit cell.
II.5.1 Sequential Algorithm
The sequential algorithm is inspired by finite size DMRG, in that we sweep through the unit cell, successively optimizing one tensor at a time while keeping tensors on other sites fixed. Notice that at site we however need two updated bond matrices and , in order to calculate updated from . We thus have to amend steps 7, 8 and 9 of the single-site algorithm in Table 2 by constructing and solving for two effective Hamiltonians and instead of a single one. The newly optimized tensors then get replaced in all unit cells of the infinite lattice, and contributions to the effective Hamiltonians have to be recalculated accordingly, before moving on to the next site. For a pseudocode summary see Algorithm 3 in Table 3.
One could now try to argue, that e.g. in a left to right sweep it is enough at site to calculate updated and only, and to use from the previous step at site as for calculating . This approach however fails, as the effective Hamiltonian used for calculating already contains updated , while the effective Hamiltonian used for calculating does not, and we cannot determine from and . Rather, has to be recalculated using an updated effective Hamiltonian, which exactly leads to the sequential Algorithm 3.
There is an additional subtlety that needs to be considered, in order for all tensors to fulfill the gauge constraints (27c) to current precision. Bond matrices are calculated as lowest energy eigenvectors of effective Hamiltonians and are therefore only determined up to a phase. Consider defined between sites and . At step it is updated as and used to calculate . In the next step however it is recalculated as (with an updated effective Hamiltonian) and used to determine . At the fixed point we should then have , but this is only true if there is no phase ambiguity, which would also consequently lead to a phase mismatch between and after step . This issue does not pose a problem for algorithm convergence (during calculations, matrices always appear as products of the form or and mismatching phases thus cancel out), but can be easily circumvented by employing a phase convention when calculating updated .
II.5.2 Parallel Algorithm
In the parallel approach, we choose to update an entire unit cell at once, using effective Hamiltonians generated from the same current state. To that end, we first generate all terms necessary for all and . For the case of nearest neighbor interactions, the contributions and to the left and right environment outside the unit cell can be shared, so that the corresponding geometric sum only needs to be computed once, and contributions inside the unit cell are obtained through successive applications of transfer matrices.
Next, we simultaneously and independently solve for the ground states and of all effective Hamiltonians at once. Once these are obtained we again simultaneously and independently determine all updated and , concluding one iteration for updating the entire unit cell. For a pseudocode summary see Algorithm 4 in Table 3.
II.5.3 Juxtaposition of Both Approaches
Several comments on the two presented algorithms are in order. First, the parallel algorithm requires substantially less computational effort, since the construction of the different effective Hamiltonians can recycle the calculation of the infinite geometric sum. Therefore, updating an entire unit cell only requires to evaluate two infinite geometric sums and effective eigenvalue problems. In the sequential algorithm, updating the environment after every tensor update requires to reevaluate the geometric sum, thus leading to infinite geometric sums and effective eigenvalue problems for updating the complete unit cell. Additionally, the parallel approach offers the possibility of parallelizing the solution of all eigenvalue problems in one iteration, while in the sequential approach only eigenvalue problems can be solved in parallel for each site. However, while sweeping through the unit cell in the sequential approach, initial guesses for solving the infinite geometric sums can be generated easily from the previous iterations, and are usually much better than the initial guesses in the parallel algorithm. Equivalently, updated obtained at site is a very good initial guess for its recalculation with updated environment on site . Overall, the computational cost for the parallel update is still much cheaper, albeit less than expected.
On the other hand, state convergence in terms of iterations is generally substantially faster in the sequential approach. This seems reasonable, as the optimization on a current site takes into account all previous optimization steps, whereas in the parallel approach, the optimizations on different sites within one iteration are independent of each other. This effect gets amplified with increasing unit cell size , and the performance of the parallel approach decreases, while the performance of the sequential approach seems more stable against increasing .
In conclusion, while updating the entire unit cell is computationally cheaper in the parallel approach, the sequential algorithm usually requires a substantially smaller number of iterations due to faster convergence. While there are instances where one approach clearly outperforms the other by far, such cases are rare and strongly depend on initial conditions, and generally both approaches show comparable performance. For comparison benchmark results see Sec. III.2.5.
III Test Cases and Comparison
In this section we test the performance of the new algorithm on several paradigmatic strongly correlated lattice models in the thermodynamic limit, with nearest neighbor as well as long range interactions. In Sec. III.1 we introduce and discuss the models under considerations. In Sec. III.2 we first test the convergence and stability of the single and multi-site implementations of the new algorithm. Lastly, we compare its performance against established conventional MPS methods for ground state search in Sec. III.3.
III.1 Models
As examples for spin chain models with nearest neighbor interactions we study the spin transverse field Ising (TFI) model
[TABLE]
and the XXZ model for general spin
[TABLE]
Here , and are spin representations of the generators of . The ground state energies are known exactly for the TFI model,Pfeuty (1970) and for also for the XXZ model.Takahashi (1999) For the XXZ model we focus on the isotropic antiferromagnetic case and take the result of Ref. Haegeman et al., 2011 for the ground state energy for as quasi-exact result.
As a further example for a system with nearest neighbor interactions we also study the Fermi Hubbard model
[TABLE]
where , are creation and annihilation operators of electrons of spin on site , and are the particle number operators. Again, the exact ground state energy is known.Lieb and Wu (1968); Essler et al. (2005)
As an example for an exactly solvable model with (algebraically decaying) long range interactions we consider the Haldane-Shastry modelHaldane (1988); Shastry (1988)
[TABLE]
where , and are again spin representations of the generators of . In order to efficiently compute the terms of the effective Hamiltonian (see Appendix C.1), we expand the distance function in a sum of exponentials, with maximum residual less than for a fit over sites.
Finally, as a state of the art problem of current interest, we also consider the two-dimensional antiferromagnetic Heisenberg model on a cylinder of infinite length, but finite circumference
[TABLE]
where denotes the location on the cylinder with and (assuming periodic boundary conditions, i.e. ). Following popular procedure (as e.g. also done in Ref. Osorio Iregui et al., 2017), we map the model onto a one-dimensional chain following a sawtooth path along the cylinder, such that the longest range spin-spin interactions are over sites. We take the energies obtained from Loop QMC for several cylinder circumferences presented in Table I in Ref. Osorio Iregui et al., 2017 as quasi-exact ground state energies.
III.2 Performance Benchmarks
We performed convergence benchmarks for several instances of the models introduced in the previous section, using simple implementations of VUMPS for single or multi-site unit cells presented in Algorithm 2, 3 and 4, without explicitly exploiting any symmetries. Hereto, we consider firstly the error in the variational energy density
[TABLE]
as a function of the number of iterations. Here is the exact analytic (or quasi-exact numerical) ground state energy density of the model under consideration. VUMPS as formulated in Table 2 has its internal convergence measure used to determine when to stop the iteration loop, as well as to set the tolerance in the iterative solvers used within every single outer iteration. However, as a more objective quantity that measures the distance to the variational minimum, we also compute the norm of the energy gradient; it is an absolute measure of convergence which is independent of any prior iterations, as opposed to relative changes in e.g. the energy or Schmidt spectrum between iterations. We denote this quantity as , the two-norm of a tensor , which can be worked out to be given by (see Appendix A.3)
[TABLE]
The efficient and accurate computation of the gradient norm is further discussed in Appendix A.4. To obtain the energy gradient of an -site unit cell, it is equivalent to determine the gradients for each site independently and to calculate the norm of the concatenation of all gradients.
A well known-property of the variational principle is that the energy expectation value itself converges quadratically faster than the state. When the state has converged to some accuracy , the energy density has already converged to precision which can therefore be well beyond machine precision. The convergence measure does however dictate the convergence of other observables which are not diagonal in the energy eigenbasis. Note, however, that we are here referring to convergence towards the value at the variational optimum, not towards the exact value. The error between the variational optimum and the exact ground state can be quantified using e.g. the energy variance, or – in the context of DMRG – the truncation error. Both quantities are also discussed in Appendix A.4. We show results for the truncation error further down in Sec. III.2.4, and when comparing VUMPS to IDMRG and ITEBD in Sec. III.3.
We show results for examples of 3 gapped and 5 critical systems in Fig. 1. Specifically, as examples for gapped systems we considered (a) the TFI model (28) in the symmetry broken ferromagnetic phase at , (b) the isotropic Heisenberg antiferromagnet, i.e. the XXZ model (29) at and (c) the XXZ model (29) in the symmetry broken antiferromagnetic phase at . As examples for gapless systems we considered (d) the isotropic Heisenberg antiferromagnet, i.e. the XXZ model (29) at , (e) the repulsive Fermi Hubbard model (30) at and half filling and (f) the Haldane-Shastry model (31). Finally, we also show results for the isotropic Heisenberg antiferromagnet on a cylinder (32) for different widths in (g) and (h).
Out of the gapped systems, only the antiferromagnetic ground state of (c) physically breaks translation invariance by spontaneously breaking the spin-flip symmetry; we therefore choose a two-site unit cell in this case. The critical systems physically show no spontaneous symmetry breaking. However, for uMPS ground state approximations in the gapless case it is often energetically beneficial to artificially break symmetries (which are restored in the limit of infinite bond dimension). In all three cases, the optimal uMPS ground state approximation artificially breaks a SU(2) symmetry and develops antiferromagnetic order, breaking translation invariance. We therefore choose a two-site unit cell in the case of the Hubbard model (e) and the Haldane-Shastry model (f). In the case of the Heisenberg antiferromagnet (d), translation invariance can be restored through a unitary transformation by rotating every second spin by around the -axis, transforming , and the artificially symmetry broken ground state becomes ferromagnetically ordered along the and directions. We can therefore choose a single-site unit cell for (d), and the staggered magnetization along is thus zero. A similar approach could be chosen to restore translation invariance also for the gapped antiferromagnet (c) and the Hubbard model (e), but we do not choose to do so for demonstrative reasons. For (g) and (h), we choose a similar strategy as for (d) and rotate every second spin by around the -axis, allowing the usage of a unit-cell size instead of , which would be necessary to accommodate an antiferromagnetically ordered ground state.
To summarize, we used the single-site Algorithm 2 for (a), (b) and (d), and the sequential Algorithm 3 with a two-site unit cell for (c), (e) and (f) and a -site unit cell for (g) and (h). For a comparison between the sequential and parallel approach see Sec. III.2.5.
III.2.1 General Convergence
Above all we observe that VUMPS shows unprecedented fast convergence, both in the energy density and the norm of its gradient , and excellent accuracy of the final ground state approximation. Observe in Fig. 1 that in all cases the energy is already well converged to machine precision after iterations, while the state is still quite some distance from the variational optimum, according to the gradient norm . Further optimizing the state, this quantity can also be converged to essentially machine precision (even in the presence of small Schmidt values), while the energy virtually doesn’t change anymore. The resulting final state then corresponds to the variationally optimal state for the given bond dimension. This is very useful in the case where the quantum state itself is required to be accurate to high precision, e.g. when used as a starting state for real time evolution or as a starting point to compute excited states and scattering thereof.Haegeman et al. (2012); Vanderstraeten et al. (2014, 2015)
The Schmidt spectrum of the ground state of the Heisenberg antiferromagnet at , converged to gradient norm , is depicted in Fig. 2. It can be seen that the degeneracies are reproduced perfectly to the same precision, without explicitly exploiting any symmetries in the implementation of the algorithm.
In cases where the final desired bond dimension is not known beforehand, one can successively enlarge a state of some small initial bond dimension every few iterations until the state fulfills the desired criteria, e.g. current bond dimension above some threshold, truncation error (see below) or smallest Schmidt value below some threshold, etc. This strategy is particularly useful when using an implementation exploiting physical symmetries of the system, such as e.g. conservation of magnetization or particle number, as the correct number and size of the required symmetry sectors in the MPS tensors is generally not known beforehand.Zauner-Stauber et al. (tion) On the other hand, if is known beforehand, it generally appears to be more efficient to immediately start from an initial state with . The gain in computational time due to the cheaper initial iterations with small bond dimension is usually outweighed by a considerable number of required additional iterations. On the other hand, for some hard problems (e.g. the Hubbard model) stability and convergence speed can profit from a strategy of sequentially increasing the bond dimension from some small initial value.
To conclude the discussion of general convergence, we plot the evolution of the Schmidt spectrum, as well as the gradient norm and various observables vs. iteration number during a ground state optimization for the TFI model (28) in the ferromagnetic phase at in Fig. 3. During the simulation we used a sequence of bond dimensions , where we started with an initial random state with and increased the bond dimension to the next value as soon as dropped below . We chose this set of bond dimension in order to not cut any degenerate multiplets of Schmidt values (see also Sec. III.2.3). It can be seen that the high lying Schmidt values converge quite quickly, while most of the computational time goes into converging the low lying Schmidt values. Moreover, there is quite some rearrangement of the small Schmidt values every time the gradient norm reaches a local maximum during a phase of non-monotonous evolution (see also next subsection). Lastly, from the evolution of the errors of the local observables and it is apparent that they require a substantially higher bond dimension of to reach the same accuracy as the energy, which is already correct to machine precision at .
III.2.2 Different Regimes of Gradient Norm Convergence
Depending on the complexity of the model, the gradient norm shows a period of irregular non-monotonous behavior before entering a regime of monotonous convergence. This can be understood as the (random) initial state having to adapt its initial structure (e.g. the Schmidt spectrum) to the requirements of the best variational ground state approximation. Monitoring the Schmidt values during this period nicely shows how groups of Schmidt values slightly rearrange to the correct structure (see e.g. Fig. 3) . This period is usually more dominant in critical systems – as can be seen in Fig. 1 (d) - (f), where it takes iterations – and of course strongly depends on the chosen initial state. One could argue, that the jumps in parameter space caused by the algorithm during this period are too big for the state to find the correct structure quickly, hindering a fast crossover to the regime of monotonous convergence. However, an approach of preconverging the state using smaller steps through parameter space – e.g. by means of imaginary time evolution with moderate time steps – has proven to be even slower in all cases tried. Thus, the best choice is still to use VUMPS during the entire optimization process. We want to emphasize here, that we have never observed a stagnation of the algorithm during this initial regime; the algorithm always reached the monotonous regime eventually in all cases, and instances where the algorithm remains in the irregular regime for an unusually long time are rare and only occur in the case of particularly hard problems.
As soon as the gradient norm reaches the monotonous regime, it always converges exponentially fast. Surprisingly, this is true even for critical systems, where one would in principle expect algebraic convergence. This can be qualitatively understood from the theory of finite entanglement scaling,Tagliacozzo et al. (2008); Pollmann et al. (2009); Stojevic et al. (2015) which states that the MPS approximation itself introduces a small perturbation away from criticality, and thus a finite gap. However, as VUMPS improves convergence speed over existing methods (see also Sec. III.3), it is ideally suited to study critical systems via the theory of finite entanglement scaling, which still requires that one finds the optimal MPS approximation in the first place.
III.2.3 Degenerate Schmidt Values
In the presence of multiplets of degenerate Schmidt values, the convergence rate is severely affected if the smallest few Schmidt values are part of an incomplete multiplet, i.e. if the last multiplet is “cut”. In that case the algorithm still shows stable convergence, albeit at a greatly reduced rate. For an example in the TFI model see Fig. 4. This issue can be easily circumvented by ensuring that the smallest few Schmidt values are part of a complete multiplet when dynamically increasing the bond dimension, or by choosing a viable (or reducing from some) fixed initial bond dimension.
III.2.4 Energy Convergence with Bond Dimension
In a careful MPS study, variational energies obtained for different bond dimension are compared in order to extrapolate to the exact limit. This can be done by plotting the energy as a function of bond dimension against the inverse of the bond dimension . The infinite limit is then obtained by fitting with a power law form and extrapolating to . In DMRG, another popular measure for the quality of an MPS approximation is given by the truncation error or discarded weight , defined in Eq. (62). The variational energy is found to scale linearly with White and Huse (1993); Legeza and Fáth (1996) and an extrapolation to is thus generally easier and more stable. For further details on assessing the quality of the ground state approximation we refer to Appendix A.4.
We show an example for both extrapolation schemes for the isotropic Heisenberg antiferromagnet in Fig. 5. The exact ground state energy is given by , or as numerical value by to 16 digits of precision. On the left we plot the energy vs. truncation error and obtain an estimate with 9 digits of precision from a linear fit . On the right we plot the energy vs. inverse bond dimension and obtain an estimate with 10 digits of precision from a power law fit . Comparing to we observe that has an error of , while has an error .
III.2.5 Multi Site Unit Cell Implementations
Lastly we discuss and compare the performance of the sequential and parallel algorithms for multi-site unit cells presented in Sec. II.5. As the two methods differ in their convergence with the number of iterations, as well as in the computational effort for each iteration of updating the entire unit cell, we compare the rate of convergence with absolute computing time in seconds. To that end we only time operations that are absolutely necessary for each algorithm, i.e. we do not time measurements, data storage etc. We further start from the same (random) initial state and keep the bond dimension fixed throughout the entire simulation for both methods to make the simulations as comparable as possible. All calculations are performed using a non-parallelized MATLAB implementation on a single core of a standard laptop CPU.
We find that in general for gapped systems, the sequential approach outperforms the parallel approach, while in critical systems no definite statement about better performance can be made. There are instances where one algorithm takes substantially longer than the other to reach the regime of monotonous convergence, but such cases are rare and strongly depend on the model and initial state. Once both algorithms are in the monotonous regime, convergence speed in terms of absolute computing time is similar, with the sequential approach generally taking longer for each iteration, but the parallel approach generally requiring more iterations to reach convergence. Overall both approaches thus generally show comparable performance, with the sequential approach appearing to be slightly more stable and reliable in the cases we considered. Fig. 6 shows examples for two gapped and two gapless systems.
III.3 Comparison with IDMRG and ITEBD
We further benchmark the performance of VUMPS against a standard two-site IDMRG implementation,White (1992, 1993); McCulloch (2008) and a standard two-site ITEBD implementation,Vidal (2007); Hastings (2011) and compare the rate of convergence of the energy error and the norm of the energy gradient between the three methods. For VUMPS, we solve the effective eigenvalue problems in each iteration to precision with the current precision according to (24). For IDMRG we solve the effective two-site eigenvalue problem in each iteration to precision , with the current orthogonality fidelity (see Sec. III.A in Ref. McCulloch, 2008). For ITEBD we employ a fourth order Suzuki Trotter decomposition and measure every 10 time steps. We use a sequence of time steps where we decrease the time step as soon as the change in Schmidt values per unit of imaginary time drops below a certain threshold. Naturally, the strategy of time step reduction should be optimized carefully for each model under consideration, however we choose the same strategy for all example cases to maintain comparability.
We also explicitly calculate all necessary quantities for obtaining a truly variational energy (e.g. by reorthogonalizing the unit cell) and for measuring the energy gradient, even if these quantities are not necessary for the respective algorithm itself.
As the three methods differ quite substantially in the number of iterations required for convergence, as well as the computational effort for each iteration, we again compare convergence against absolute computing time in seconds, where we only time operations that are absolutely necessary for each algorithm and we do not time measurements, data storage, reorthogonalizing, etc. We further start from the same (random) initial state and keep the bond dimension fixed throughout the entire simulation for all three methods to make the simulations as comparable as possible. Again, all calculations were performed using a non-parallelized MATLAB implementation on a single core of a standard laptop CPU.
We show example comparisons for two gapped and two critical models in Fig. 7, similar to the cases studied in the previous section. Specifically, we show results for (a) the gapped isotropic Heisenberg antiferromagnet, i.e. the XXZ model (29) at , (b) the XXZ model (29) in the gapped symmetry broken antiferromagnetic phase at , (c) the critical isotropic Heisenberg antiferromagnet, i.e. the XXZ model (29) at , and (d) the critical Fermi Hubbard model (30) at and half filling. We plot the energy error on the left and the gradient norm on the right, vs. absolute computing time in seconds. For VUMPS we used a single-site unit cell for (a) and (c), and a two-site unit cell for (b) and (d).
Above all we observe that VUMPS clearly outperforms both IDMRG and ITEBD by far, both in convergence speed and accuracy of the final state, especially for critical systems. In all shown cases, the final energy error of all three algorithms only differ by a few percent; VUMPS however always yields the best variational energy, often already after a few seconds, and thus converges in energy much faster than IDMRG or ITEBD – in the case of critical systems even by orders of magnitude. Observe that especially for the two critical systems (c) and (d) a large part of the computational time of IDMRG and ITEBD goes in converging the last few digits of the energy (see also insets in Fig. 7). For (d) in particular, the final energy error obtained by IDMRG is still almost higher than the value obtained by VUMPS.
In terms of convergence of the energy gradient , we observe that IDMRG and ITEBD perform quite poorly. Surprisingly, IDMRG usually stagnates at some value . ITEBD on the other hand would be in principle capable of converging essentially also to machine precision, albeit at prohibitively long simulation times, as the limiting factor appears to be the Trotter error, requiring very small time steps; we therefore also only reach values of with ITEBD within reasonable simulation times. 444The fact that ITEBD usually yields better variational states than IDMRG for the same bond dimension has already been observed in Ref. McCulloch, 2008. VUMPS on the other hand is always capable to converge essentially to machine precision, and does so – contrary to other methods – exponentially fast and with unprecedented speed. For instance, in the case of the Hubbard model in example (d), ITEBD only reached a gradient norm of after hours of absolute computing time, while VUMPS already reached this value after only seconds, and converged further to in seconds. IDMRG on the other hand stagnated at a quite high value of .
Finally, we also compare variational ground state energies for the Heisenberg antiferromagnet on a cylinder with results from Ref. Osorio Iregui et al., 2017; Osorio Iregui, tion. We have already seen in Fig. 1 that VUMPS shows excellent convergence speed also for this model. We show a comparison of obtained energies for cylinder circumferences and bond dimensions in Table 4. There it is shown that for all cases, VUMPS yields lower energies with differences roughly around the 5th-7th significant digit.555IDMRG simulations have been performed without a rotation by of every second spin and with a unit cell to accommodate an antiferromagnetically ordered ground state.Osorio Iregui (tion) We conclude that VUMPS is thus perfectly suited for a fast and accurate study of two dimensional models of current interest on infinite cylinders and improves over current state of the art methods.
III.3.1 Observables
We also measure and compare the regular and staggered (averaged) magnetizations and of the final state after convergence for the Hubbard model in example (d), as in this case all three methods use a two-site unit cell. The exact ground state is SU(2) symmetric and thus has zero magnetization; a finite ground state approximation however artificially breaks this symmetry. The final values for the regular magnetization are zero to machine precision for both VUMPS and IDMRG, but for ITEBD. The staggered magnetization is for VUMPS, for IDMRG, and for ITEBD. Both the regular and staggered magnetizations are thus smallest for the final state obtained from VUMPS. For IDMRG the staggered magnetization is highest, but the regular magnetization is zero, which in turn is finite for ITEBD. This result is not surprising, as VUMPS yields the best variational state out of the three methods.
III.3.2 Truncation Error
As a last figure of merit, popular in DMRG studies as a measure of the quality of the MPS ground state approximation and used for extrapolations to the exact infinite limit, we also calculate the truncation error or discarded weight of the final state, defined in Eq. (62). In the case of the Hubbard model in example (d), we obtain a truncation error from IDMRG, and a slightly lower from VUMPS.
IV Conclusion and Outlook
We have introduced a novel algorithm for calculating MPS ground state approximations of strongly correlated one dimensional quantum lattices models with nearest neighbor or long range interactions, in the thermodynamic limit. It combines ideas from conventional DMRG and tangent space methods by variationally optimizing a uniform MPS by successive solutions of effective eigenvalue problems. The algorithm can easily be implemented by extending an existing single-site (I)DMRG implementation with routines for i) calculating effective Hamiltonian contributions from infinite environments and ii) solving an effective “zero site” eigenvalue problem in addition to the usual single-site problem. The new algorithm is free of any ill-conditioned inverses and therefore does not suffer from small Schmidt values, contrary to other tangent space methods such as TDVP. Additionally, as it does not rely on imaginary time evolution, it is especially fit for studying systems with long range interactions.
We described and benchmarked implementations for uniform MPS with both single-site and multi-site unit cells. We observed that the new algorithm clearly outperforms existing methods such as IDMRG and ITEBD, both in convergence speed and accuracy of the final state at convergence. This is especially also the case for systems with long range interactions or simulations of two dimensional models on infinite cylinders. The energy converges with unprecedented speed after iterations, even in critical systems (where this is orders of magnitude faster than conventional methods). The algorithm further proceeds to converge the state to the variational optimum by minimizing the energy gradient essentially to machine precision; it does so exponentially fast, even for critical systems, contrary to other methods. The new algorithm is thus the perfect choice for studying critical systems. Additionally, a state converged to the variational optimum is particularly useful in cases where the quantum state itself is required to be accurate to high precision, e.g. when used as a starting state for real time evolution or for a variational calculation of elementary excitations.Haegeman et al. (2012); Vanderstraeten et al. (2014, 2015)
It is straightforward to include physical symmetries that come with good quantum numbers (such as e.g. conserved magnetization or particle number) after a proper definition of a symmetric uniform MPS unit cell, where absolute (diverging) values of these quantum numbers are replaced by densities. All steps of the algorithm then immediately also apply to MPS tensors with good quantum numbers. Symmetric ground states obtained this way are an excellent starting point for obtaining elementary excitations with well defined quantum numbers following Ref. Haegeman et al., 2012, which for instance enables to target elementary excitation that lie within a multi-particle continuum.Zauner-Stauber et al. (tion)
Within the same framework it is also very natural to recover real or imaginary time evolution by replacing the effective Hamiltonian ground state problems (23a) and (23b) by small finite time evolution steps, which yields the thermodynamic limit version of the time evolution algorithm presented in Ref. Haegeman et al., 2016. This enables e.g. to efficiently study real time evolution of quantum states on systems with long range interactions in the thermodynamic limit.Halimeh and Zauner-Stauber (2016); Halimeh et al. (2017)
We believe that the ideas presented in this paper should be relevant for other classes of tensor-network states as well. Specifically in the case of projected entangled pair states (PEPS) Verstraete and Cirac (2004), designed to capture ground states on two-dimensional quantum lattices, the further development of efficient variational algorithms – as an alternative to ITEBD inspired approaches – is still much desired.Corboz (2016); Vanderstraeten et al. (2016) In particular, it motivates the search for (approximate) canonical forms for PEPS, which would enable a translation of the VUMPS algorithm to the two-dimensional setting. Equivalently, an adaptation of the new algorithm for continuous Matrix Product State (cMPS) Verstraete and Cirac (2010) would allow for a much desired increase of efficiency for cMPS simulations, building on recent advances in that direction.Ganahl et al. (2017)
acknowledgments
The authors acknowledge inspiring and insightful discussions with M. Ganahl, M. Mariën and I.P. McCulloch. We thank B. Buyens for providing templates for drawing tensor network diagrams. We also thank J. Osorio Iregui for providing IDMRG data for infinite cylinder simulations. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 647905). The authors gratefully acknowledge support from the Austrian Science Fund (FWF): F4104 SFB ViCoM and F4014 SFB FoQuS (V.Z.-S. and F.V.), and GRW 1-N36 (M.F.). J.H. and L.V. are supporte d by the Research Foundation Flanders (FWO).
Appendix A Theoretical background
In this Appendix we reiterate definitions and concepts needed for the algorithm presented in Sec. II in more detail, and motivate the VUMPS algorithm from a variational perspective.
A.1 Variational principle on manifolds
The variational principle in quantum mechanics characterizes the ground state of a given Hamiltonian as the state which minimizes the normalized energy expectation value
[TABLE]
If (typically for computational reasons) we only have access to a subset of Hilbert space, the variational principle still gives a way to find an approximation to the true ground state, namely by solving the minimization problem within the restricted set. If this subset is a linear subspace spanned by a number of basis vectors , we obtain a generalized eigenvalue problem
[TABLE]
for the expansion coefficients in . This is known as the Rayleigh-Ritz method, and by orthonormalizing the basis it clearly amounts to projecting the full time-independent Schrödinger equation into the variational subspace.
If, more generally, we have a variational ansatz which depends analytically on a number of complex parameters, as encoded in the complex vector , a variational minimum is characterized by a vanishing gradient of the energy expectation value, i.e.
[TABLE]
with the (formally independent) complex conjugate of , and the complex derivatives with respect to the ’th component of and and
[TABLE]
Eq. (35) can be interpreted as a Galerkin condition: It forces the residual of the full Schrödinger eigenvalue equation – which does not have an exact solution in the variational subset – to be orthogonal to the space spanned by the states . If the variational subset is a manifold, these states can be interpreted as a basis for the tangent space of the manifold at the point of the variational optimum. Hence, geometrically, the residual has to be orthogonal to the manifold (and thus to its tangent space) at the point of the variational optimum. Interpreting Eq. (35) as a Galerkin condition on the ground state eigenvalue problem is useful because it can be generalized to other eigenvalue problems which do not necessarily have a variational characterization (and thus no gradient), as e.g. when the operator is non-Hermitian. Indeed, a similar approach as is developed here was described for finding fixed points of transfer matrices, encoded as matrix product operators, in Ref. Haegeman and Verstraete, 2016.
However, before discussing Eq. (35) in the context of MPS, let us conclude this section by relating it to the time-dependent variational principle (TDVP).Haegeman et al. (2011) Geometrically, the TDVP also amounts to an orthogonal projection of the equation of motion (the time-dependent Schrödinger equation) onto the tangent space of the variational manifold. In the case of imaginary time evolution, it can be written as
[TABLE]
where
[TABLE]
is the Gram matrix of the tangent vectors and thus the metric of the manifold. The right hand side of Eq. (36) is again the gradient of the objective function, and the TDVP will thus converge when it reaches a variational optimum where Eq. (35) is satisfied. However, the metric in the left hand side shows that the TDVP equation is not a normal gradient flow, but rather a proper covariant gradient flow that takes the geometry of the manifold and its embedding into the Hilbert space into account. We can thus also associate a quantum state with the gradient, which is given by
[TABLE]
where we have omitted the arguments and , is the inverse of the metric and
[TABLE]
is the projector onto the tangent space at the point in the variational manifold . This latter expression is only valid when the variational parameters are proper coordinates for the manifold (i.e. a bijective mapping). While this is not the case for MPS because of gauge freedom (see below), the geometrical interpretation for the Galerkin condition
[TABLE]
remains valid; the correct expression for the MPS tangent-space projector will be discussed in more detail in the following section. Independent of whether a variational algorithm is based on a gradient flow, the Hilbert space norm of the gradient provides an objective measure for the convergence of the state towards the variational optimum. Note that this is different from the standard Euclidean norm of the naive gradient vector with components .
A.2 Manifold of Uniform MPS and its tangent space
We have already introduced the set of uniform MPS and some of its properties in Sec II.1. For completeness, we restate the definition as
[TABLE]
Such states can be rigorously defined in the thermodynamic limit and correspond to the so-called purely generated finitely correlated states of Ref. Fannes et al., 1992. However, for all practical purposes, we can interpret this state as a large but finite MPS, which happens to be uniformly parameterized (and therefore translation invariant) in the bulk. We have now introduced boundary vectors living at . They represent just one way to close the matrix product near the boundaries, but any other behavior is equally fine and won’t affect the bulk properties, provided that the MPS tensor has the property of injectivity.Perez-Garcia et al. (2007) This condition is generically fulfilled and is in one-to-one correspondence with the MPS transfer matrix
[TABLE]
having a unique eigenvalue of largest magnitude, with corresponding eigenvectors that can be reinterpreted as full rank positive matrices. A proper normalization of the state is obtained by rescaling the tensor such that this largest magnitude eigenvalue is .
We can then use gauge invariance to transform the tensor into the left and right canonical representations (2) discussed in the main text. These representations are themselves related via the bond matrix as , which allows to write the state in the mixed canonical representation (5) familiar from DMRG. We can also make contact with the representation commonly used in TEBDVidal (2003, 2007) by additional unitary gauge transforms. We make diagonal by considering the SVD and gauge transform , . The singular values are then the Schmidt values of a bipartition of the state and we can obtain , or equivalently
[TABLE]
and .
It was proven in Ref. Haegeman et al., 2014 that the set of injective uMPS constitute a (complex) manifold. We now construct the tangent space projector for a uMPS in the thermodynamic limit.
Let us first discuss generic tangent vectors in this tangent space. To define them via the partial derivatives of , we require that the latter is represented using a uniform parameterization or gauge choice of the tensor throughout. By applying the chain rule, the derivative will give rise to a uniform superposition of states where a single tensor is replaced by a new tensor .Haegeman et al. (2013a, 2014) But it is clear that we can afterwards change gauges again and absorb the gauge factors in the tensor that parameterizes the tangent vector, so as to obtain the most general tangent vector representation
[TABLE]
The multiplicative gauge freedom of the MPS translates into an additive gauge freedom in the tangent space, i.e. a transformation with leaves invariant, as can readily be verified by explicit substitution. We can exploit these gauge degrees of freedom to impose e.g. the left tangent space gauge
[TABLE]
or the right tangent space gauge
[TABLE]
Strictly speaking, these conditions can only be imposed for tangent vectors . This is no restriction as we can always evaluate the contribution in the direction of separately. Under either of these two gauge constraints, we indeed have , but more importantly, the overlap between two tangent vectors simplifies to
[TABLE]
This corresponds to an Euclidean inner product for the tensors and thus to an orthonormal basis for the tangent space. The diverging factor arises because a tangent vector contains a sum over all lattice sites, i.e. its norm is extensive. Fortunately, this diverging factor will drop out in all computations.
Given Eq. (39) we need to derive the explicit form of the projector onto the tangent space, a derivation that was written down in Ref. Haegeman et al., 2016 for the case of finite MPS. The tangent vector resulting from the orthogonal projection of a general translation invariant state onto the tangent space can be readily found by solving the minimization problem
[TABLE]
or, equivalently,
[TABLE]
In order to use Eq. (46) for the first term, we however need to impose the constraint in Eq. (44) or (45). In the former case, this will add a term to the objective function, with and corresponding Lagrange multipliers. The solution is readily obtained by demanding , where Eq. (46) simply results in . The overlap between a tangent vector and is given by
[TABLE]
so that its derivative is easily obtained by omitting the tensor from the diagram and interpreting the open legs as defining the indices of a new tensor. Without the Lagrange multiplier, we would simply obtain
[TABLE]
With the additional constraint, the solution is still straightforward. It can be easily verified that the correct value of the Lagrange multiplier is such that the additional term acts as a projection
[TABLE]
or similarly
[TABLE]
if we would have chosen the right gauge of Eq. (45).
By inserting the solution for back into Eq. (43), we can read off the tangent space projector. While the value of depends on the gauge condition, the resulting projector is of course gauge independent and given by
[TABLE]
We can represent the tangent space projector as
[TABLE]
by defining the partial projectors
[TABLE]
We can verify that by using
[TABLE]
A.3 Gradient and Effective Hamiltonians
As discussed in the beginning of this section, a variational optimum can be characterized geometrically as
[TABLE]
where (unit normalization is assumed). Since the Galerkin condition is automatically ensured in the direction of the MPS itself, the only non-trivial information of Eq. (52) is thus contained in the part of the tangent space orthogonal to . This is convenient, as was actually constructed as the projector onto the part of tangent space orthogonal to in the first place. While this implies that the subtraction does not contribute, it is convenient to keep it around, as it ensures that the individual terms in the final expression are finite in the thermodynamic limit.
Applying the tangent space projection as in the previous section to the state gives rise to a tangent vector of the form Eq. (43), with
[TABLE]
where originates from applying and from applying . By writing itself in a compatible gauge for every individual term, we obtain the diagrammatic expressions
[TABLE]
and
[TABLE]
We can thus obtain and by acting with the effective Hamiltonians and introduced in (9) and (10) in the main text onto and
[TABLE]
Even without subtracting the energy, the two choices of (which are related by the additive gauge transform with ) will be finite in the thermodynamic limit. However, the individual tensors and will have a divergent contribution proportional to and , respectively. Indeed, as discussed in the main text for the case of nearest neighbor interactions, the effective Hamiltonians and have a divergent contribution corresponding to the total energy times the identity operator. It is thus by subtracting (or for the local terms) that these divergences are canceled. Appendix C provides a detailed description of the construction of the effective Hamiltonian for other types of interactions and illustrates explicitly that the diverging contributions cancel exactly.
A variational extremum is characterized by , which leads to for either gauge choice, as these choices completely fix the gauge freedom. This gives rise to the following simultaneous conditions:
[TABLE]
However, because the gauge transformation that relates and is unique up to a factor (for injective MPS), and have to be proportional, and we have actually obtained the eigenvalue equations
[TABLE]
As we are looking for a variational minimum, the eigenvalues and should be the lowest eigenvalues of the effective Hamiltonians and . Depending on how we have regularized the divergent contributions, these eigenvalues might be different. If we have completely subtracted the energy expectation value from every term, we then have .
A.4 Convergence and error measures
While neither VUMPS nor IDMRG or ITEBD directly use the gradient itself, the Hilbert space norm of can be used as an objective convergence measure to indicate how far the current state is from the variational optimum. For either choice of , we obtain , with the diverging number of sites and the 2-norm of the tensor . Its square is given by
[TABLE]
where the equalities follow from . Note that none of these expressions are well suited for numerically evaluating the norm close to convergence, as they involve subtracting quantities that are almost equal, especially when the state is close to convergence.
An alternative strategy for evaluating is by using the matrix notation for tensors (17) to write . Since is an isometry, we can extend it to a unitary matrix , where contains an orthonormal basis for the -dimensional null space of , i.e. . As the 2-norm is unitarily invariant, we can write
[TABLE]
The second equality follows from and the third from the null space property of . We then obtain as the Frobenius norm of a single matrix, which can be calculated accurately as a sum of strictly positive numbers. For further reference, we reshape into a tensor and similarly introduce a tensor via the defining relations
[TABLE]
While the norm of the gradient provides a measure for the quality of approaching the variational minimum, it does not provide any information about the quality of the (u)MPS approximation to the true ground state itself. In the context of (two-site) DMRG schemes, a popular measure is the truncation error, as it is naturally accessible throughout the algorithm (see e.g. Ref. White and Huse, 1993; Legeza and Fáth, 1996; Schollwöck, 2011). But also within VUMPS we can compute this quantity by first writing the state in a mixed canonical form with a two-site center block
[TABLE]
The two-cite center tensor (known as the two-site wave function in standard DMRG) has an associated effective Hamiltonian . We can compute its lowest eigenvector and compute its singular value decomposition (by first reshaping it to a matrix) . The truncation error then corresponds to the discarded weight
[TABLE]
when truncating the inner bond dimension of this two-site tensor to its original value .
A more generic measure for the error in the variational approximation is given by the energy variance . This quantity is also used in the context of e.g. variational Monte Carlo and various other methods. We can systematically decompose into various parts: the projection onto is automatically zero by the definition of . The projection onto the tangent space is zero when we are at the variational minimum. Next, we can project onto the space of all 2-site variations, which is given by states of the form
[TABLE]
In the case of nearest neighbor Hamiltonians, this space captures completely, namely by choosing with the two-site center tensor defined in the previous paragraph and the local terms of the Hamiltonian, with the current expectation value subtracted, i.e. . However, there is again an additive representation redundancy (gauge freedom) which enables us to choose representations satisfying e.g. a left gauge condition (). The advantage of this representation is again that it facilitates the calculation of the norm, as with the diverging number of sites. The projection of onto this space can be worked out similarly as for the tangent space, and leads to the general result (for any Hamiltonian)
[TABLE]
with a single application of the two-site effective Hamiltonian. Using , we can rewrite the first form of as
[TABLE]
We now also apply to the right hand side and recognize . But since at the variational minimum, we obtain at the variational minimum
[TABLE]
and in particular
[TABLE]
We can also relate to the truncation error defined in the previous paragraph. For the truncation error arising in the context of two-site DMRG schemes, the lowest eigenvector of the two-site effective Hamiltonian is used. When the DMRG algorithm has converged, the rank approximation of should again be the original two-site center tensor . But this means that we can construct and exactly from the singular vectors corresponding to the singular values that were truncated away, and thus that the truncation error is given by . This definition is close to , except that the latter uses the tensor arising from applying the two-site effective Hamiltonian once. As and are anyway close, we can think of as providing the leading order correction from to in the sense of a Krylov scheme. Indeed, in the first iteration of the Lanczos method, the eigenvector would be approximated in the form . Since the first term drops out when projecting onto and , the DMRG truncation error and will be of the same order of magnitude.
Note, however, that only captures the full energy variance (per site) for nearest neighbor Hamiltonians, whose action on is completely contained within the space of two-site variations as noticed above. In that case, we can see that the only term that survives in after projection onto and is the local term, where acts on the two-site center tensor. We can thus also write
[TABLE]
We can also relate this to the truncation step in (I)TEBD, where we would apply to every two-site block of the state. The resulting truncation would lead to a discarded weight of the order .
The considerations regarding the projection of onto the space of two-site variations can also be used to devise a scheme for expanding the bond dimension of the uMPS. This approach is presented in the next section.
Appendix B Dynamic Control of the Bond Dimension
A characteristic feature of two-site implementations of conventional MPS methods – such as e.g. (I)TEBD or (I)DMRG – is that the bond dimension of the MPS is automatically increased in every iteration and has to be truncated in order to remain at a finite maximum bond dimension. This truncation step lies at the basis of why such schemes for finding ground states will never truly converge to the variational minimum up to machine precision, as observed in the results. Indeed, even in finite size simulations, two-site DMRG is used to initialize the state and one-site DMRG to obtain final convergence. However, the truncation step in two-site methods has the advantage that the bond dimension can be dynamically increased (or decreased) according to some quality constraint, such as the magnitude of the smallest Schmidt-value or the discarded weight. Especially in the presence of symmetry, this is important to automatically obtain the correct symmetry sectors within the virtual MPS space.
The VUMPS algorithm presented in the main text is variational from the start and therefore works at fixed bond dimension , i.e. it is a one-site scheme in DMRG terminology. Alternative subspace expansion strategies for dynamically increasing the bond dimension in such one-site schemes have been proposed.White (2005); Hubig et al. (2015) These methods use information from acting with the global Hamiltonian onto the current state to either add a tiny perturbation to the current MPS or to generate a larger basis in which the effective eigenvalue problem is solved.
We have developed a similar subspace expansion technique that works for a uMPS in the thermodynamic limit. It is based on projecting the full action of the Hamiltonian onto the space of two-site variations, as developed in the previous section. There we have found the representation
[TABLE]
Even when we have not yet reached the variational minimum, the first term (on line 1) or the first two terms (on line 2) are captured in the tangent space, and only the last term (on either line) contains a new search direction. To capture it completely, we would need to expand the bond dimension from value to . If we want to expand to a new dimension , we can use a singular value decomposition to compute the rank approximation of . By keeping only the largest singular values, and are left and right isometries of size and respectively. As remarked in the previous section, in the case of nearest neighbor interactions, the projection of onto and does not require the full two-site effective Hamiltonian but reduces to the local term.
We do not directly update the current MPS, but rather write it in an en expanded basis in a mixed canonical form with matrices
[TABLE]
With these initial tensors, we can now start a new iteration of VUMPS. Note that we can straightforwardly update the environments used to construct the effective Hamiltonians into this expanded basis, which is necessary if we want to use them as initial guess.
Appendix C Explicit Construction of Effective Hamiltonians
In this section we describe how to efficiently apply the effective Hamiltonians and onto the center site tensor and bond matrix and how the necessary individual terms are explicitly constructed. Such a procedure is needed for solving the effective eigenvalue problems (23a) and (23b) by means of an iterative eigensolver. The case of systems with nearest neighbor interaction has already been discussed in Sec. II.2. In the following we consider the cases of Hamiltonians with long range interactions in Sec. C.1 and general Hamiltonians given in terms of Matrix Product Operators (MPOs) in Sec. C.2.
C.1 Long Range Interactions
Consider Hamiltonians with long range interactions of the form , where is itself an infinite sum
[TABLE]
and operators act on a single-site and commute when acting on different sites . Without loss of generality, we restrict to a single pair of (bounded) operators , which commute when acting on different sites 666For fermionic Hamiltonians with long range interactions, we can employ a Jordan-Wigner transformation to spin operators, introducing a string operator counting the number of fermions between and ; geometric sums of transfer matrices in (70) then turn into geometric sums of (string) operator transfer matrices.. The generalization to Hamiltonians containing several terms of that form is straight forward. Furthermore, we assume distance functions that are bounded in the sense of , such that , and that can be well approximated by a sum of exponentials, i.e.
[TABLE]
with and . In practice, for an infinite system we fit with a suitable number of exponentials over a distance large enough, such that and the largest residuals are below some desired threshold.
Examples of Hamiltonians that fall in this class are the transverse field Ising (TFI) model or XXZ model with power-law interactions,Ruelle (1968); Dyson (1969); Cardy (1981) as well as the famous Haldane-Shastry model,Haldane (1988); Shastry (1988) for which the ground state is exactly known.
Similar to the case of nearest neighbor interactions in Sec. II.2, the effective Hamiltonians factorize into a number of terms which can all be applied efficiently. For these are five terms, out of which four are already familiar from the case of nearest neighbor interactions. Two of these are the left and right block Hamiltonians and with infinitely many local contributions from acting on sites strictly left or right of the current center site, and the other two are the terms containing interactions between the center site and the left and right block respectively, i.e. where partially acts on . For long range interactions we have one additional term, containing infinitely many interaction terms between the left and the right block only without involving the center site, i.e. where acts to the left of the current center site, and acts to the right.
To construct all these terms we start by defining the operator transfer matrices
[TABLE]
The current energy density expectation value can thus be written as
[TABLE]
or using (68)
[TABLE]
Since the geometric series converge and we can perform them explicitly. We proceed by defining
[TABLE]
These terms can again either be calculated recursively by explicitly evaluating the geometric sums term by term until convergence, or more efficiently by iteratively solving the following systems of linear equations
[TABLE]
using iterative methods.
We represent these terms by the diagrams
[TABLE]
and collect all such terms into single left and right environment contributions
[TABLE]
and further
[TABLE]
We can then write for the energy density
[TABLE]
Comparing with (70) we have thus defined
[TABLE]
With these definitions at hand we can write the left and right block Hamiltonians as
[TABLE]
These equations are exactly the same as Eq. (13) for the case of nearest neighbor interactions, but with different and . We can thus evaluate the geometric sums recursively or by solving a linear system iteratively, as explained in Sec. II.2. Note that we again start by applying an energy shift and similar for , such that .
We are now ready to formulate the action of onto as
[TABLE]
The additional factor of in the sum in the last term arises due to adding an additional site between the left and right operators . Similarly, the action of onto becomes
[TABLE]
In (79) the first two terms can be applied in , the second two in and the last term in operations, and in (80) the first two terms in and the last term in operations. In general we have to perform iterative inversions involving operations and collect terms to arrive at the necessary terms for (79) and (80), where the solutions from the previous iteration can be used as starting vectors to speed up convergence.
If there are additional simple single or nearest neighbor two-site terms present in the Hamiltonian, appropriate terms as described in Sec. II.2 can be added. For a pseudocode summary for obtaining the necessary explicit terms of and for Hamiltonians with long range interactions, and their applications onto a state, required for solving the effective eigenvalue problems using an iterative eigensolver, see Table 5.
C.2 General Hamiltonians given in terms of MPOs
Consider the Hamiltonian given in terms of an infinite Matrix Product Operator (MPO) Verstraete et al. (2004b); McCulloch (2007); Pirvu et al. (2010); Crosswhite and Bacon (2008); Crosswhite et al. (2008); Fröwis et al. (2010); Michel and McCulloch (2010) with 4-index MPO elements with and and we call the MPO bond dimension. In terms of the operator valued matrices the Hamiltonian can then be written as
[TABLE]
where contains operators acting on site only and and are operator valued boundary vectors.
An example for such an MPO decomposition for the Transverse Field Ising (TFI) Hamiltonian with exponentially decaying long range interaction
[TABLE]
with is given by
[TABLE]
where and are Pauli matrices. For the TFI Hamiltonian we thus have and the limit corresponds to the nearest neighbor interaction case.
In order to efficiently apply the effective Hamiltonians and , it is necessary to determine the left and right (quasi) fixed points and of the MPO transfer matrices
[TABLE]
where – similar to MPS tensors – and are collections of matrices of dimension , with . These two objects are in fact the thermodynamic limit versions of the objects defined in Eq. (190) and (191) in Ref. Schollwöck, 2011.
Typically, MPO representations of (quasi)local Hamiltonians (such as e.g. Eq. (81)) are of Schur form,Michel and McCulloch (2010) such that the MPO transfer matrix contains Jordan blocks and that the dominant eigenvalue is one and of twofold algebraic degeneracy. Such MPO transfer matrices therefore technically do not have well defined fixed points. We can however find quasi fixed points and , that are fixed points up to a term contributing to the energy density expectation value in one of the elements of and . An application of onto both quasi fixed points will therefore accumulate an additional term contributing to the extensive global energy expectation value. Similar to the terms in Eq. (13) or Eq. (78) in the previous cases involving infinite geometric sums, we can however safely discard these diverging contributions, which is equivalent to setting the energy expectation values of the semi-infinite left and right half of the system to zero (see below).
In the following we briefly reiterate the procedure of Ref. Michel and McCulloch, 2010 to systematically determine and from given and . The obtained solutions will of course contain the results of Sec. II.2 and Sec. C.1 as special cases.
Without loss of generality we assume to be of lower triangular form, i.e. . Furthermore, we assume the typical case of any nonzero diagonal elements being proportional to the identity, i.e. , where and usually , as is e.g. the case in (81). 777The generalization to nontrivial operators on the diagonal, such as e.g. Pauli- operators in the case of long range fermion hopping, is straight forward. By defining the result of the action of the MPO transfer matrix as
[TABLE]
the system of fixed point equations can be written as
[TABLE]
Notice that due to the lower triangular structure of , the terms and only contain contributions from and and we can solve (85) and (86) recursively, starting with for and with for , which initially amounts to and . Terms with are particularly simple and simply reduce to the identification and .
Terms with where , now result in solutions of the form
[TABLE]
equivalent to terms such as (72) stemming from infinite geometric sums of (weighted) MPS transfer matrices.
Equivalently, terms with then result in relations of the form
[TABLE]
which in general do not have a formal solution, since the left hand sides of these equations live in the subspace orthogonal to the dominant eigenspaces of , while the right hand sides generally do have contributions in the dominant eigenspace. We can however discard these contributions by projecting onto the complementary subspace, and then obtain and by solving the systems of equations (see also Appendix D)
[TABLE]
We have encountered exactly the same type of equations in (15) when evaluating infinite geometric sums of transfer matrices, after a constant shift in energy to remove diverging terms. The MPO formalism thus automatically yields these contributions in a form where the sums have already been explicitly performed.
Such a situation typically occurs only for the final terms in the recursive solution of the fixed point equations, i.e. for and . A concrete evaluation (see below) of the discarded terms in these cases shows that they correspond to contributions to the energy density expectation value, i.e. discarding these terms is equivalent to a constant shift in energy, such that the energy density is zero and we have . After applying once onto the quasi fixed points we thus have for the first element of and the last element of
[TABLE]
i.e. the fixed point relations only hold up to an additive diagonal correction for these elements. These corrections correspond to the energy density expectation value
[TABLE]
and they can in fact be used for its evaluation.
As a concrete example, for the long range TFI Hamiltonian given by MPO (81) we obtain
[TABLE]
and
[TABLE]
Having determined the left and right quasi fixed points of the MPO transfer matrices, it is now particularly easy to calculate the action of the effective Hamiltonians onto as
[TABLE]
and equivalently the action of onto as
[TABLE]
which can be performed in , respective operations. In total we also have to perform an iterative inversion for each diagonal element of .
This framework is very flexible, general and powerful, once a routine for determining quasi fixed points of general MPO transfer matrices has been implemented. The effective Hamiltonians of Sec. II.2 and Sec. C.1 are contained within as special cases. A pseudocode summary for obtaining the necessary explicit terms of and for Hamiltonians given in terms of an MPO, and their applications onto a state, required for solving the effective eigenvalue problems using an iterative eigensolver, is presented in Table 6.
Appendix D Geometric Sums of Transfer Matrices
We wish to evaluate terms involving infinite geometric sums of the form
[TABLE]
Such expressions typically arise in situations where one sums up contributions of successive applications of onto some fixed virtual boundary vector , with the initial contribution being the boundary vector itself. This is reflected in the above expression by summing from and using the definition .
We assume a spectral decomposition of the transfer matrix given by
[TABLE]
where the left and right eigenvectors are mutually orthonormal, i.e. . Note that is in general not hermitian and thus .
For a generic injective normalized state, has a unique eigenvalue of largest magnitude given by , whereas all other eigenvalues are contained in the unit circle ().
We divide into dominant and complementary subspaces and get for powers of
[TABLE]
We can safely perform the geometric sum for all eigenvalues , while contributes a formally diverging term
[TABLE]
The interpretation of this diverging contribution depends on the situation. By using the projectors
[TABLE]
onto the dominant and complementary subspaces we define the projected transfer matrix
[TABLE]
We realize that the spectral decomposition of has a component of
[TABLE]
and therefore identify the second term in (100) as
[TABLE]
For the geometric sum we then obtain
[TABLE]
with a diverging contribution from . Plugging into (96) we finally get
[TABLE]
Usually it is not necessary to calculate the full matrix expression of , but to just act with it onto some or . The diverging contributions can typically be safely discarded, as they correspond to a constant (albeit infinite) offset of some extensive observable (e.g. the Hamiltonian). The action of the finite remaining part can be calculated efficiently by iteratively solving the linear system of equations of the type or
[TABLE]
with inhomogeneities and . One can then efficiently compute and by employing an iterative Krylov subspace method such as bicgstabvan der Vorst (1992) or gmres.Saad and Schultz (1986) For such methods only the implementation of a (left or right) action of onto a vector is necessary, which can be done efficiently with operations. If the transfer matrix is in left or right canonical form, we recover the linear systems in Eqs. (15) and (91)
[TABLE]
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Wilson (1975) K. G. Wilson, Rev. Mod. Phys. 47 , 773 (1975) .
- 2Fisher (1974) M. E. Fisher, Rev. Mod. Phys. 46 , 597 (1974) .
- 3Fisher (1998) M. E. Fisher, Rev. Mod. Phys. 70 , 653 (1998) .
- 4Zinn-Justin (1996) J. Zinn-Justin, Quantum Field Theory and Critical Phenomena (Oxford University Press, 1996).
- 5White (1992) S. R. White, Phys. Rev. Lett. 69 , 2863 (1992) .
- 6White (1993) S. R. White, Phys. Rev. B 48 , 10345 (1993) .
- 7Fannes et al. (1992) M. Fannes, B. Nachtergaele, and R. Werner, Comm. Math. Phys. 144 , 443 (1992) .
- 8Klümper et al. (1993) A. Klümper, A. Schadschneider, and J. Zittartz, Europhys. Lett. 24 , 293 (1993) .
