Table Of Content

Variational optimization algorithms for uniform matrix product states V. Zauner-Stauber,1 L. Vanderstraeten,2 M.T. Fishman,3 F. Verstraete,1,2 and J. Haegeman2 1Vienna Center for Quantum Technology, University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria 2Ghent University, Faculty of Physics, Krijgslaan 281, 9000 Gent, Belgium 3Institute for Quantum Information and Matter, California Institute of Technology, Pasadena, California 91125, USA 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 di- mensional 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 inconvergencespeedand precision. We alsodemonstrate thatVUMPSworksveryefficiently forHamiltonianswithlongrangeinteractions. Thenewalgorithmcanbeconvenientlyimplemented as an extension of an already existing DMRG implementation. I. INTRODUCTION Using a translation invariant parameterization gives rise 7 to an energy expectation value with a highly non-linear 1 dependence on the tensor(s). Two different algorithms Thestrategyofrenormalizationgroup(RG)techniques 0 are widely used to obtain such an MPS in the thermo- 2 to successively reduce a large number of microscopic de- dynamic limit. Infinite system DMRG (IDMRG)5,6,24 grees of freedom to a smaller set of effective degrees of n proceeds by performing regular DMRG ona successively freedom has led to powerful numerical and analytical a growing lattice, inserting and optimizing over new ten- J methods to probe and understand the effective macro- sors in the center of the lattice in each step only, effec- scopic behavior of both classical and quantum many 4 tivelymimickinganinfinitelatticebyusingafinite,albeit 2 body systems.1–4 However, it was not until the advent very large lattice. After convergence the most recently of White’s celebrated Density Matrix Renormalization inserted tensors in the center are taken as a unit cell for ] Group (DMRG)5,6 that variational RG methods reached h an infinite MPS approximation of the ground state. An unprecedentedaccuracyinnumericallystudyingstrongly p alternative approach is known as infinite time evolving correlated one dimensional quantum lattice systems at - block decimation (ITEBD).25,26 It works directly in the t low temperature. The underlying variational ansatz of n thermodynamic limit and is based on evolving an initial a Matrix Product States (MPS)7–13 belongs to a class of stateinimaginarytimebyusingaTrotterdecomposition u ansatzes known as Tensor Network States.11,14,15 These of the evolution operator. q variational classes encode the many body wavefunction [ We present a new variational algorithm, inspired by in terms of virtual entanglement degrees of freedom liv- 1 ing on the boundary and thus satisfy an area law scaling tangentspaceideas,13,27,28 thatcombinestheadvantages ofIDMRGandITEBDandaddressessomeoftheirshort- v of entanglement entropy per construction. As such, they comings. Assuchitisdirectlyformulatedinthethermo- 5 provide a natural parameterization for the physical cor- 3 nerofHilbertspace,wherelowenergystatesofquantum dynamic limit, but at the same time optimizes the state 0 manybodysystemsoughttolivein.16,17 MPSinparticu- bysolvingeffectiveeigenvalueproblems,ratherthanem- 7 ployingimaginarytimeevolution. Wefindthatitleadsto larareespeciallyfitforstudyinggroundstatesofstrongly 0 a significant increase in efficiency in all of our test cases. correlated one dimensional quantum systems with local 1. interactions.18–20 ThefollowingsectionintroducesMPSnotationsanddef- 0 initions and presents our variational algorithm, heuristi- 7 The variational parameters in MPS are contained callymotivatedfromtheperspectiveoffinitesizeDMRG. 1 withinlocaltensorsassociatedwiththeindividualsitesof Sec. III illustrates the performance of our algorithm on : the lattice system. For homogeneous systems, the global v varioustestcases,andcomparestoconventionalIDMRG wave function can then be captured using just a single i and ITEBD results. After the conclusion in Sec. IV, we X (or a small number of) such tensors, independent of the provide further technical details in the appendices. Ap- r system size. They consequently offer very natural access pendixAcontainsadditionaltheoreticalbackground: we a tothethermodynamiclimit, providingaclearadvantage derivetheself-consistentconditionsthatcharacterizethe overothernumericalapproachessuchasExactDiagonal- variational minimum and provide additional motivation ization or Quantum Monte Carlo. for our algorithm from the perspective of the MPS tan- On finite lattices, (one-site) DMRG implements the gent space. Appendix B presents a suitable strategy to variational principle (energy minimization) by exploit- expandthebonddimensionoftranslationinvariantMPS. ing that the quantum state is a multilinear function of AppendixCexplainshowtoconstructeffectiveHamilto- the local tensors. By fixing all but one tensors, the nians in the thermodynamic limit. These involve infinite global eigenvalue problem is transformed into an effec- geometric sums of the transfer matrix, which are further tive eigenvalue problem for the local tensor.5,6,12,21–23 studied in Appendix D. 2 II. A VARIATIONAL ALGORITHM FOR A. Uniform MPS MATRIX PRODUCT STATES IN THE THERMODYNAMIC LIMIT A uniform MPS (uMPS) of bond dimension D defined on an infinite translation invariant lattice is parameter- In this section we introduce a variational algorithm ized by a single collection of d matrices As ∈ CD×D for foroptimizingMPSdirectlyinthethermodynamiclimit. s=1,...,d. Theoveralltranslationinvariantvariational Because the algorithm strongly resembles conventional state is then given by DMRG,weexplainitbydescribingasingleiterationstep (cid:88)(cid:16) (cid:17) from the viewpoint of DMRG and show that only a few |Ψ(A)(cid:105)= ...Asn−1AsnAsn+1... |s(cid:105) (1) additional ingredients are needed to arrive at our vari- s ational algorithm. We only briefly motivate these extra and can be represented diagrammatically as ingredients for the sake of readability, and refer to Ap- pendix A for additional explanations and more rigorous |Ψ(A)(cid:105)=... A A A A A ... theoretical motivations. As such, the new algorithm can easily be implemented as an extension to an already ex- isting (I)DMRG implementation. Exploiting the invariance of (1) under local gauge transformationsAs →XAsX−1,withX ∈CD×D invert- We start by considering a setting familiar from con- ible,thestatecanbecastintocertainfavorablerepresen- ventional DMRG: a finite homogeneous one dimensional tations, among them the left and right canonical repre- quantum lattice system, where every site corresponds to sentation a d level system. We label the sites by an integer n and thushaveabasis{|s(cid:105)n,s=1,...,d}forthelocalHilbert (cid:88)AsL†AsL =11 (cid:88)AsLRAsL† =R (2a) spaceonsiten. ThetotalHilbertspaceisspannedbythe (cid:78) s s product basis |s(cid:105) = |s(cid:105) . We assume the dynamics (cid:88) (cid:88) n n AsAs† =11 As†LAs =L, (2b) of the system to be governed by a translation invariant R R R R Hamiltonian H. s s We further consider a variational parameterization of or diagrammatically a ground state approximation of the system, for now in terms of a finite size (site dependent) MPS, but we AL AL willultimatelybeinterestedinthethermodynamiclimit. DMRGproceedstofindthebestvariationalgroundstate = R = R approximationbyemployinganalternatingleastsquares minimization: It starts from some initial state and suc- A¯L A¯L cessively optimizes each of the individual MPS tensors site by site by solving effective (Hamiltonian) eigen- AR AR value problems, in a sweeping process through the lat- = L = L tice until convergence, where each iteration depends on already optimized tensors from previous iterations (see e.g. Refs 5, 6, 12, 21, and 23). A¯R A¯R We are now however interested in the thermodynamic Here L and R correspond to the left and right reduced limit n ∈ Z (but will ignore the technical complications densitymatricesofabipartitionofthestaterespectively. involving a rigorous definition of a Hilbert space in that WehenceforthrefertoA (A )asaleft(right)isometric L R limit). In that case the MPS ground state approxima- tensor, or just a left (right) isometry. tion will be given in terms of a translation invariant uni- Defining the left and right transfer matrices form MPS, described by a single MPS tensor (or a unit cell of N tensors), repeated on all sites. Two immediate T =(cid:88)A¯s ⊗As (3) L/R L/R L/R difficulties arise: Firstly, conventional DMRG updates s the variational state site by site, thus breaking transla- and using the notation (x| and |x) for vectorizations of tion invariance. Secondly, the effective Hamiltonian for a D×D matrix x in the D2 dimensional “double layer” a single-site optimization has to be constructed from an virtual space the transfer matrices act upon, the gauge infinite environment. conditions (2) are equivalent to After briefly introducing the variational class of uni- form MPS and introducing necessary notation and con- (11|T =(11| T |R)=|R) (4a) L L ventions (for further details see Sec. A2), we describe T |11)=|11) (L|T =(L|, (4b) R R how the new algorithm modifies DMRG accordingly to exactly account for these two issues in order to arrive at i.e. 11 and R are the left and right dominant eigenvectors a variational ground state algorithm directly formulated of T , while L and 11 are the left and right dominant L in the thermodynamic limit. eigenvectors of T . R 3 AsisstandardpracticeinDMRG,wemixbothofthese states representations and cast the state into the mixed canon- |Ψ(α,s,β)(cid:105)=|Ψα(cid:105)|s(cid:105)|Ψβ(cid:105) (8a) ical representation AC L R |Ψ(α,β)(cid:105)=|Ψα(cid:105)|Ψβ(cid:105). (8b) C L R |Ψ(A)(cid:105)=(cid:88)(...Asn−1AsnAsn+1...)|s(cid:105) (5a) L C R s B. Effective Hamiltonian =(cid:88)(...Asn−1AsnCAsn+1Asn+2...)|s(cid:105), (5b) L L R R The use of the mixed canonical representation (5a) in s DMRG is of significant importance for the stability, as it reduces the minimization of the (global) energy ex- or diagrammatically, pectationvalue(cid:104)Ψ|H|Ψ(cid:105)/(cid:104)Ψ|Ψ(cid:105)withrespecttoA into C a standard (hermitian) eigenvalue problem, instead of a generalizedone. TheeffectiveHamiltonianforthiseigen- value problem is the system Hamiltonian H projected |Ψ(A)(cid:105)=... AL AL AC AR AR ... onto the degrees of freedom of AC, and is known as the “reduced” or “superblock” Hamiltonian in DMRG. We define the thermodynamic limit version of this re- =... AL AL AL C AR AR ... duced single-site Hamiltonian acting on AC as H (α(cid:48),s(cid:48),β(cid:48)) =(cid:104)Ψ(α(cid:48),s(cid:48),β(cid:48))|H|Ψ(α,s,β)(cid:105) (9) AC(α,s,β) AC AC Here we have defined the center site tensor As (known AL AL AR AR C as the single-site wave function Ψs in DMRG) =··· H ··· AsC =AsLC =CAsR A¯L A¯L A¯R A¯R (6) AC = AL C = C AR Additionally, we also define an effective Hamiltonian acting on the bond matrix C as H (α(cid:48),β(cid:48)) =(cid:104)Ψ(α(cid:48),β(cid:48))|H|Ψ(α,β)(cid:105) (10) C(α,β) C C intermsofthebondmatrixC,whichconstitutesthe(in- vertible) gauge transformation relating A and A via AL AL AR AR L R As = CAsC−1. The singular values of C then encode thLe entangRlement spectrum of the state. Indeed, using =··· H ··· AsC = CAs we can verify that the left and right re- duLced densitRy matrices in (2) are given by L=C†C and A¯L A¯L A¯R A¯R R = CC†. Furthermore, AsC = CAs ensures that (5a) L R which does not appear directly in the context of DMRG, and (5b) are translation invariant and that A and C C but will be needed later for a consistent update of the can be shifted around arbitrarily. Normalization of the state without breaking translation invariance. It can be state,aswellasofthereduceddensitymatricesLandR, interpreted as a “zero site” effective Hamiltonian, which correspondstothesinglecondition(cid:107)C(cid:107)2 =Tr(CC†)=1. 2 would feature in an optimization of the global energy Foreaseofnotationwefurtherintroducethefollowing expectation value with respect to the Schmidt values. partial states In an efficient implementation, these effective eigen- value problems are typically solved using an iterative eigensolver,sothatweonlyneedtoimplementtheaction |ΨαL(cid:105)=(cid:88)(...AsLn−1AsLn)α|...sn−1sn(cid:105) (7a) of HAC and HC onto AC and C. While the energy expectation value is extensive and s thus divergent in the thermodynamic limit, the effective =... AL AL AL α Hamiltonians H and H are well defined and finite in AC C the thermodynamic limit if one properly subtracts the |Ψα(cid:105)=(cid:88)(AsnAsn+1...) |s s ...(cid:105) (7b) current energy expectation value from the Hamiltonian R R R α n n+1 s H. Wedemonstratethisprocedureforthecaseofnearest (cid:80) = α AR AR AR ... neighborinteractionsH = nhn,n+1,wherethetwo-site Hamiltonian h acts on neighboring sites n,n + 1 n,n+1 only. We refer to Appendix C for the case of long range interactions and for general Hamiltonians given in terms withnarbitrary,andusethemtodefinethereducedbasis of Matrix Product Operators (MPOs). 4 Inthecaseofnearestneighborinteractions, theactionofH ontoA splitsupintofourindividualcontributions, AC C which follow from the decomposition |Ψ(cid:105)=(cid:80) As |Ψα(cid:105)|s(cid:105)|Ψβ(cid:105) (left block containing sites n<0, center site α,β,s C,(α,β) L R n=0, and right block containing sites n>0). The action of H onto A is given by AC C A(cid:48)s =(cid:88)htsAt†AkA(cid:96) +hstAkA(cid:96) At †+H As +AsH C k(cid:96) L L C k(cid:96) C R R L C C R tk(cid:96) AL AC AC AR AC AC (11) A(cid:48)C = h + h + HL + HR A¯L A¯R where the first two terms correspond to the Hamiltonian terms h and h coupling the center site to the left and −1,0 0,1 rightblock,respectively,andH andH sumupthecontributionsofalltheHamiltoniantermsh actingstrictly L R n,n+1 to the left and to the right of the center site. TheenvironmentsH andH areusuallyconstructed where (H | can be presented diagrammatically as L R L iteratively while sweeping through the (finite) lattice in   conventional DMRG, or grown successively in every it- AL AL AL eration in IDMRG. In the thermodynamic limit, these     terms consist of a diverging number of individual local HL = hL 11+ + + ... interaction contributions h , and care needs to be   n,n+1 taken in their construction. A¯L A¯L A¯L Indeed, the kth contribution to (H | comes from the L and likewise for |H ). Hamiltonian term h and is given by (h |[T ]k−1. R −k−1,−k L L The transfer matrix T has a dominant eigenvalue of Likewise, [T ]k−1|h ) is the kth contribution to |H ) L R R R magnitude one, with corresponding left and right eigen- stemming from h . Here, we have used the defini- k,k+1 vectors (11| and |R). The projection (h |[T ]k|R) = tions L L (h |R) is the energy density expectation value e = L hL =(cid:88)hskt(cid:96)AtL†AsL†AkLA(cid:96)L (cid:104)tΨhe|he−nke−rg1,y−kh˜|Ψ=(cid:105) ahn−d eis11infrdoempetnhdeenHtaomfilkt.oniSaunb,twraectcinang stk(cid:96) (12) write (h | = (h˜ | + e(11|. The second term is exactly h =(cid:88)hstAkA(cid:96) At †As†, proportioLnal to tLhe left eigenvector of eigenvalue 1 and R k(cid:96) R R R R stk(cid:96) therefore gives rise to a diverging contribution in the ge- ometric sum, corresponding to the total energy of the or diagrammatically left half infinite block. Since this contribution acts as the identity in the effective Hamiltonian H [Eq. (11)], AC wecanhoweversafelydiscardthisdivergingcontribution AL AL without changing the eigenvectors of H . This corre- AC sponds to an overall energy shift of the left half infinite hL = h blocksuchthat(H |R)=0. Fortheremainingpart(h˜ | L L A¯L A¯L the geometric sum converges. With |h˜R) = |hR)−e|11) the same comments apply to the construction of |H ). R AR AR We can evaluate HL and HR recursively as (H[n+1]|=(H[n]|T +(h˜ | hR = h L L L L (14) |H[n+1])=T |H[n])+|h˜ ) R R R R A¯R A¯R with initialization (H[0]| = (h˜ | and |H[0]) = |h˜ ). We L L R R Summing up all such local contributions gives rise to in- can repeat these recursions until e.g. (cid:107)H[n+1] −H[n] (cid:107) L/R L/R finite geometric sums of the transfer matrices TL/R drops below some desired accuracy (cid:15)S. This strategy is conceptually simple and closely resembles the succes- ∞ ∞ sive construction of the environments in the context of (cid:88) (cid:88) (HL|=(hL| [TL]k |HR)= [TR]k|hR), (13) (I)DMRG, but is not very efficient, as its performance is k=0 k=0 comparable to that of a power method. 5 Algorithm 1 Explicit terms of effective Hamiltonians with nearest neighbor interactions and their application Input: two-site Hamiltonian h, current uMPS tensors A , A in left and right gauge, left dominant eigenvector (L| of T , L R R right dominant eigenvector |R) of T , desired precision (cid:15) for terms involving infinite geometric sums L S Output: Explicit terms of effective Hamiltonians H and H , updated A(cid:48) and C(cid:48) AC C C 1: function HeffTerms(H =h,AL,AR,L,R,(cid:15)S) (cid:46) Calculates explicit terms of effective Hamiltonians 2: Calculate hL and hR from (12) 3: Calculate HL and HR by iteratively solving (14) or (preferably) (15), to precision (cid:15)S 4: HAC ←{h,AL,AR,HL,HR} 5: HC ←{h,AL,AR,HL,HR} 6: return HAC,HC 7: end function 8: function ApplyHAC(AC,HAC) (cid:46) Terms of HAC from HeffTerms(H,AL,AR,L,R,(cid:15)S) 9: Calculate updated A(cid:48) from (11) C 10: return A(cid:48) C 11: end function 12: function ApplyHC(C,HC) (cid:46) Terms of HC from HeffTerms(H,AL,AR,L,R,(cid:15)S) 13: Calculate updated C(cid:48) from (16) 14: return C(cid:48) 15: end function Table I. Pseudocode for obtaining the explicit terms of the effective Hamiltonians H and H for systems with nearest AC C neighbor interactions and their applications onto a state. A more efficient approach is to formally perform the for (H | and |H ) to precision (cid:15) , as explained in Ap- L R S geometricsumsin(13)explicitly, andtoiterativelysolve pendix D. the resulting two systems of equations So far, we have discussed the action of H . The ac- AC tionofH ontoC followssimplyfrom(11)byprojecting C (HL|[11−TL+|R)(11|]=(hL|−(hL|R)(11| ontoALorAR. UsingthedefiningpropertyofHLorHR, (15) the result simplifies to [11−T +|11)(L|]|H )=|h )−|11)(L|h ) R R R R C(cid:48) =(cid:88)hstAs†AkCA(cid:96) At †+H C+CH . k(cid:96) L L R R L R stk(cid:96) AL C AR C C (16) C(cid:48) = h + HL + HR A¯L A¯R The first two terms of (11) can be applied in O(d4D3) C. Updating the state operations29, and the last two terms in O(dD3) opera- tions. For (16) the first term can be applied in O(d4D3) InDMRG,wewouldupdatethestatebyreplacingA operations, and the last two terms in O(D3) operations. C with the lowest eigenvector A˜ of H and then shift Togeneratethenecessarytermsfor(11)and(16)wehave C AC the center site to the right by computing an orthogonal to iteratively evaluate two infinite geometric sums, in- factorization A˜s = A˜sC˜ , or to the left by computing volving O(D3) operations (when iteratively solving (15) C L R A˜s = C˜ A˜s. As such, the state gets updated by only the solutions from the previous iteration can be used as C L R startingvectorstospeedupconvergence). Apseudocode replacing the current site with A˜sL or A˜sR, while leaving summary for obtaining the necessary explicit terms of allothersitesuntouched. However,applyingthisscheme H and H and their applications onto a state is pre- in our setting would immediately destroy translation in- seAntCed in TaCble I. variance after a single step. We want to construct an alternative scheme that ap- plies global updates in order to preserve translation in- variance at any time. Global updates can most easily 6 be applied with an explicit uniform parameterization in where matrices P are hermitian and positive. Alter- termsofasingletensorA. Ontheotherhand,DMRGex- native isometric decompositions might be considered in perienceteachesusthatthestabilityisgreatlyenhanced Eq. (21), though it is important that they are unique when applying updates at the level of A and C, which (e.g. QR with positive diagonal in R) in order to have C are isometrically related to the full state. P[(cid:96)/r] ≈P[(cid:96)/r] close to convergence. We therefore calculate the lowest eigenvector A˜ of AC C C H like in DMRG, but additionally also the lowest AC eigenvector C˜ of HC. We then globally update the D. The Algorithm: VUMPS state by finding new A˜ and A˜ as the left and right L R isometric tensors that minimize (cid:80) (cid:107)A˜sC˜ − A˜s(cid:107)2 and s L C Wearenowreadytoformulateourvariationaluniform (cid:80)s(cid:107)C˜A˜sR−A˜sC(cid:107)2 respectively. Theseminimizationprob- MPS (VUMPS) algorithm. As shown in Appendix A, a lems can be solved directly (not iteratively) and without variational minimum (vanishing energy gradient) in the invertingC˜ (seebelow). AsshowninAppendixA,atthe manifoldofuMPSischaracterizedbytensorsA , C and L variational optimum the values of these objective func- A satisfying the conditions R tions go to zero, and current A and C are the lowest C eigenvectors of H and H respectively. H A =E A (23a) AC C AC C AC C For the remainder of this section we omit tildes and H C =E C (23b) C C use the following matricizations of the 3-index tensors As =AsC =CAs. (23c) C L R A =As L,(sα,β) L,(α,β) HereboldsymbolsdenotevectorizationsoftheMPSten- AR,(α,sβ) =AsR,(α,β) (17) sorsandmatricizationsoftheeffectiveHamiltonians,and A[C(cid:96)],(sα,β) =AC[r,](α,sβ) =AsC,(α,β). EHAamCilatnodniaEnCs.3a1re the lowest eigenvalues of the effective We thus want to extract updated AL and AR from up- When iterating the steps outlined in the previous sub- dated AC and C by solving sections, convergence is obtained when these conditions are satisfied. In particular, starting with a properly or- (cid:15) = min (cid:107)A[(cid:96)]−A C(cid:107) (18a) L C L 2 thogonalized initial trial state |Ψ(A)(cid:105) of some bond di- A†LAL=11 mensionD,webeginbysolvingthetwoeigenvalueprob- (cid:15) = min (cid:107)A[r]−CA (cid:107) . (18b) lems for the effective Hamiltonians H and H . Since R ARA†R=11 C R 2 we are still far from the fixed point, thAeCresultinCg lowest energy states A˜ and C˜ will in general not satisfy the In exact arithmetic, the solution of these minimization C gauge condition (23c) together with current A . problems is known, namely A will be the isometry in L/R L the polar decomposition of A[l]C† (and similar for A , Following the procedure of the previous section we C R can however find optimal approximations A˜s and A˜s see Thm. IX.7.2 in Ref. 30). Computing the singular L R for (23c) to arrive at an updated uMPS. Conversely, A˜ value decompositions (SVD) C and C˜ will not be the correct lowest energy eigenstates A[C(cid:96)]C† =U[(cid:96)]Σ[(cid:96)]V[(cid:96)]† C†AC[r] =U[r]Σ[r]V[r]† (19) of the new effective Hamiltonians HA˜C and HC˜ gener- ated from A˜ . We then use the updated state and we thus obtain L/R reiterate this process of alternately solving the effective A =U[(cid:96)]V[(cid:96)]† A =U[r]V[r]†. (20) eigenvalueproblems,andfindingoptimalapproximations L R for A and A to update the state. For a pseudocode L R Notice that close to (or at) an exact solution As = summary of this algorithm, see Table II. C AsC =CAs, the singular values contained in Σ[(cid:96)/r] are We now elaborate on the various steps in the VUMPS L R the square of the singular values of C, and might well algorithm. Firstly, extracting new A˜ from updated L/R fallbelowmachineprecision. Consequently,infinitepre- A˜ and C˜ can be done using the theoretically opti- C cisionarithmetic,correspondingsingularvectorswillnot mal (but numerically often inaccurate) Eq. (20) or the be accurately computed. morerobustEq.(22),dependingonthemagnitudeofthe An alternative that has proven to be robust and still smallest singular value in C˜. As a good uMPS approxi- close to optimal is given by directly using the following mationwillalwaysinvolvesmallsingularvalues,Eq.(22) left and right polar decompositions is preferable most of the time, except maybe during the A[(cid:96)] =U[(cid:96)] P[(cid:96)] C =U[(cid:96)]P[(cid:96)] (21a) first few iterations. C AC AC C C The maximum of the error quantities (18) A[r] =P[r]U[r] C =P[r]U[r] (21b) C AC AC C C (cid:15) =max((cid:15) ,(cid:15) ) (24) prec L R to obtain A =U[(cid:96)] U[(cid:96)]† A =U[r]†U[r], (22) providesanerrormeasureforthefixedpointconditionin L AC C R C AC Eq. (23c) and is used as a global convergence criterion. 7 Algorithm 2 variational uMPS algorithm for single-site unit cells Input: Hamiltonian H, initial uMPS A , A , C, convergence threshold (cid:15) L R Output: uMPS approximation A , A , C of ground state of H, fulfilling fixed point relations (23a), (23b) and (23c) up to L R precision (cid:15) 1: procedure VUMPS(H,AL,AR,C,(cid:15)) 2: initialize current precision (cid:15)prec >(cid:15) 3: while (cid:15)prec >(cid:15) do 4: (optional) Dynamically adjust bond dimension following Appendix B 5: Calculate explicit terms of effective Hamiltonians HAC,HC ←HeffTerms(H,AL,AR,L,R,(cid:15)S ≤(cid:15)prec) from Algorithm 1, 5 or 6 6: Calculate ground state A˜C of effective Hamiltonian HAC to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHAC(A ,H ) from Algorithm 1, 5 or 6 C AC 7: Calculate ground state C˜ of effective Hamiltonian HC to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 C 8: Calculate new A˜L and A˜R from A˜C and C˜ using (20) or (22), depending on singular values of C˜ 9: Evaluate new (cid:15)L and (cid:15)R from (18) 10: (optional) Calculate current expectation values 11: Set (cid:15)prec ←max((cid:15)L,(cid:15)R) and replace AL ←A˜L, AR ←A˜R and C ←C˜ 12: end while 13: return AL, AR, C 14: end procedure Table II. Pseudocode of the VUMPS algorithm described in Sec. IID. Terms within step 5 involving the evaluation of infinite geometric sums usually require the left dominant eigenvector L of T and the right dominant eigenvector R of T , for which R L L = C†C and R = CC† with current C are a good enough approximation to current precision (cid:15) (see main text). Notice prec that this algorithm is free of any possibly ill-conditioned inverses and therefore has no convergence issues in the presence of small Schmidt values. It also does not require expensive reorthogonalizations of the state at intermediate iterations. It measures the precision of the current uMPS ground prohibitive. Whenthisapproachconverges,theresulting state approximation. Within every iteration, we use it- fixedpointtensorsinthecentercanbeassumedtospecify erative methods (e.g. some variation of Lanczos) to find theunitcellofaninfiniteMPS.VUMPShastheimmedi- the eigenvectors A˜ and C˜ of the Hermitian operators ate advantage that i) it directly works in the thermody- C H and H . As the goal is to drive the state towards namiclimitatalliterationsandii)itcompletelyreplaces AC C the fixed point relations in Eqs. (23a) and (23b), it is the entire state after every iteration, thus moving faster not necessary to solve these eigenvalue problems to full through the variational manifold. In contrast, IDMRG machine precision. Rather, it is sufficient to use a toler- keeps memory of earlier iterations and cannot guarantee ance (cid:15) chosen relative to (cid:15) .32 A value of (cid:15) of the a monotonically decreasing energy that converges to an H prec H order of (cid:15) /100 has proven to work well in practice. It optimum associated with a translation invariant MPS in prec is also worthwhile to use tensors from the previous iter- whichtheeffectsoftheboundaryhavecompletelydisap- ation as initial guess for the iterative solvers to speed up peared. The advantages of VUMPS come with a greater convergence. computationalcostperiteration,astwoeigenvalueprob- Asthemainpartofthealgorithmworksatfixedbond lems(forAC andforC)and–inthecaseofnearestneigh- dimension(i.e.itisasingle-siteschemeinDMRGtermi- bor interactions – two linear systems (for HL and HR) nology),onemightchoosetoincreasethebonddimension havetobesolved. IDMRGonlysolvesasingleeigenvalue D before starting a new iteration. We have developed a problem and builds HL and HR step by step in every it- subspace expansion technique that works directly in the eration. The latter approach is analogous to a power thermodynamic limit and is explained in Appendix B. method for eigenvalue problems and, while very cheap, While the true comparison of this algorithm with is expected to require many iteration steps to converge, IDMRG5,24 and ITEBD26 will take place in Sec. III by especially for systems with large correlation lengths (e.g. gathering actual numerical simulation results, we can al- close to criticality). ready compare the theoretical properties of these algo- ITEBD26 isbasedonevolvinganinitialstateinimagi- rithms. Neither IDMRG or ITEBD is truly solving the narytimebyusingaTrotterdecompositionoftheevolu- variationalprobleminthesenseofdirectlytryingtosat- tion operator. Like VUMPS, ITEBD works in the ther- isfy the fixed point conditions Eqs. (23). IDMRG closely modynamiclimitatanyintermediatestep,typicallywith resembles regular DMRG on a successively growing lat- a unit cell that depends on how the Hamiltonian was tice, as it inserts and optimizes over new tensors in the split into local terms in order to apply the Trotter de- center of the lattice in each step. Tensors from previ- composition. Furthermore, as every application of the ous steps are not updated, as this would render the cost evolution operator increases the virtual dimension of the 8 MPS, truncation steps are required to restore the origi- algorithmstendtoconvergetomaximallysymmetrybro- nal(oranysuitable)valueofthebonddimension. While ken states for which the entanglement is minimal. This VUMPScantakebigstepsthroughthevariationalspace, isalsothecaseforVUMPS.Onecancontrolwhichstate timestepsinITEBDhavetobechosensufficientlysmall the algorithm converges to by suitably biasing the ini- (especially in the final steps of the algorithm) to elimi- tial state or by adding small perturbation terms to the natetheTrottererror,whichnegativelyaffectstherateof Hamiltonian which explicitly break the symmetry, and convergence. (Ref.33howeverproposesaschemetoeffec- which are switched off after a few iterations. tively obtain a larger time step). Furthermore, the Trot- Explicit conservation of translation symmetry was the tersplittingessentiallylimitstheapplicabilityofITEBD very first requirement in the construction of VUMPS. to short-range interactions and dictates the size of the In the case of spontaneous breaking of translation sym- unit cell of the resulting MPS, e.g. in the most common metry down to N-site translation symmetry (as e.g. in case of nearest neighbor interactions a two-site unit cell the case of a state with antiferromagnetic order), en- isobtained. (TheapproachofRef.34toobtainatransla- forcing one-site translation symmetry would result in a tioninvariantMPSisrestrictedtocertainHamiltonians, (non-injective) equal weight superposition of all symme- but see Ref. 35 for an alternative proposal that can in try broken uMPS ground state approximations. In order fact also deal with long range interactions.) to reach an optimal accuracy with a given bond dimen- Finally, we can also compare VUMPS to the more sion,suchasuperpositionofN statesishoweverundesir- recent time dependent variational principle (TDVP),27 able, astheeffectivebonddimensionisreducedtoD/N. which was implemented as an alternative approach to In the case where this situation cannot be amended by simulate real and imaginary time evolution within the a simple unitary transformation that restores one-site manifold of MPS by projecting the evolution direction translation symmetry (such as e.g. flipping every second onto the MPS tangent space. This approach can be ap- spin in the case of an antiferromagnet), it is preferable plied to translation invariant MPS, independent of the to choose an MPS ansatz with a N-site unit cell, such type of Hamiltonian. When used to evolve in imaginary that the state can spontaneously break translation sym- time, it can be identified as a covariant formulation of metry. The generalization of the algorithm to multi-site a gradient descent method, in that it evolves the state unit cells is described in the next section. in the direction of the gradient of the energy functional, preconditionedwiththemetricofthemanifold. Assuch, the energy decreases monotonically and at convergence, E. Multi Site Unit Cell Implementations an exact (local) minimum is obtained, as characterized by the vanishing gradient. However, in its original for- We now generalize the VUMPS algorithm of the pre- mulation,TDVPwasnotformulatedinacentersiteform vious section for one-site translation invariant uMPS to and was therefore unstable and restricted to small time the setting of translation invariance over N sites. Such steps. For finite systems, a different formulation of the a uMPS ansatz is then parameterized by N independent TDVP algorithm was provided in Ref. 28, which allows tensorsA(k)s ∈CD×d×D,k =1,...,N,whichdefinethe fortakingthelimitoftheimaginarytimesteptoinfinite, unit cell tensor and then becomes provably equivalent to the single-site DMRGalgorithm. VUMPScanbemotivatedfromthese As =A(1)s1...A(N)sN, (25) developments, as explained in Appendix A. We conclude this section by elaborating on how to in- wheres=(s1,...,sN)isacombinedindex. Wecanthen corporate symmetries in the algorithm. The construc- write the variational state as tuinointaoryf usyMmPmSettrhiaets iiss eeqxupilvicailtelnytintova(rIi)aDntMuRnGd1e2r,2o3nasnitde |Ψ(A)(cid:105)=(cid:88)(...Asn−1AsnAsn+1...)|s(cid:105) (I)TEBD36,37, and it is immediately clear that the var- s ious steps in VUMPS have a corresponding covariant andtheleftandrightorthonormalformsaregivenbythe formulation. The same comments apply to time rever- relations sal symmetry, in which case everything can be imple- mented in real arithmetic, or to reflection symmetry, in (cid:88)A(k)s†A(k)s =11 which case C and As will be symmetric matrices and L L C s AsR = AsLT (which implies that HL and HR are also re- (cid:88)A(k)s R(k)A(k)s† =R(k−1) (26a) lated). In all of these cases, the computational cost is L L s reduced. However, explicitly imposing the symmetry in the MPS requires caution, as the physical system might and have spontaneous symmetry breaking, or – more subtly – might be in a symmetry protected topological phase (cid:88)A(k)s A(k)s† =11 R R where the symmetries cannot be represented trivially on s (26b) the MPS tensor. (cid:88)A(k)s†L(k−1)A(k)s =L(k), R R In the case of spontaneous symmetry breaking, MPS s 9 where it is understood that N +1≡1 and 0≡N. determine A˜(k) from A˜(k) and C˜(k −1) . Rather, R C R Defining the bond matrices C(k) as the gauge trans- C˜(k) has to be recalculated using an updated effective L formation that relates left and right canonical form via Hamiltonian, whichexactly leadsto thesequentialAlgo- C(k−1)A(k)s =A(k)sC(k),wehaveR(k)=C(k)C(k)† rithm 3. R L and L(k)=C(k)†C(k). We can then also cast |Ψ(A)(cid:105) in Thereisanadditionalsubtletythatneedstobeconsid- amixedcanonicalformsimilarto(5a)withthecentersite ered,inorderforalltensorstofulfillthegaugeconstraints tensor given by A(k)s =A(k)sC(k)=C(k−1)A(k)s. (27c) to current precision. Bond matrices C˜(k) are cal- C L R The variational minimum within this set of states is culated as lowest energy eigenvectors of effective Hamil- characterized by the following 3N fixed point relations tonians H and are therefore only determined up to a C(k) phase. Consider C(k) defined between sites k and k+1. HA(k)CA(k)C =EA(k)CA(k)C (27a) At step k it is updated as C˜(k)R and used to calculate HC(k)C(k)=EC(k)C(k) (27b) A˜(k)s. Inthenextstepk+1howeveritisrecalculatedas L A(k)s =A(k)sC(k)=C(k−1)A(k)s. (27c) C˜(k+1) (with an updated effective Hamiltonian) and C L R L usedtodetermineA˜(k+1)s. Atthefixedpointweshould R Notice that due to (27c), the relations for different k are then have C˜(k) = C˜(k+1) = C(k), but this is only connected. There are several possible strategies for con- R L trueifthereisnophaseambiguity,whichwouldalsocon- structing algorithms which obtain states satisfying these sequently lead to a phase mismatch between A˜(k) and conditions. L C˜(k)afterstepk+1. Thisissuedoesnotposeaproblem Inthefollowingwepresenttwoapproacheswhichhave for algorithm convergence (during calculations, matrices shown good performance and stable convergence, which C(k) always appear as products of the form C(k)†C(k) we shall term the “sequential” and “parallel” methods. or C(k)C(k)† and mismatching phases thus cancel out), But let us first elaborate on computing effective Hamil- but can be easily circumvented by employing a phase tonians for multi-site unit cells, which works similarly in convention when calculating updated C˜(k). both methods. We again restrict to the case of nearest neighbor interactions, such that the effective Hamiltoni- ans are constructed similar as in Sec. IIB. To construct 2. Parallel Algorithm e.g. the left block Hamiltonian H , we first collect all L local contributions from a single unit cell in h , before L performing the geometric series of the transfer matrix, Intheparallelapproach,wechoosetoupdateanentire whichnowmediatesatranslationoveranentireunitcell. unit cell at once, using effective Hamiltonians generated from the same current state. To that end, we first gen- erate all terms necessary for all H and H . For A(k)C C(k) 1. Sequential Algorithm the case of nearest neighbor interactions, the contribu- tions H and H to the left and right environment out- L R The sequential algorithm is inspired by finite size side the unit cell can be shared, so that the correspond- DMRG, in that we sweep through the unit cell, succes- ing geometric sum only needs to be computed once, and sively optimizing one tensor at a time while keeping ten- contributions inside the unit cell are obtained through sorsonothersitesfixed. Noticethatatsitekwehowever successive applications of transfer matrices. need two updated bond matrices C˜(k) =C˜(k−1) and Next, we simultaneously and independently solve for C˜(k) = C˜(k), in order to calculate Lupdated A˜(k)s the ground states A˜(k)C and C˜(k) of all 2N effective R L/R Hamiltonians at once. Once these are obtained we again from A˜(k)s ≈ A˜(k)s C˜(k) ≈ C˜(k) A˜(k)s. We thus C L R L R simultaneouslyandindependentlydetermineallupdated havetoamendsteps5,6and7ofthesingle-sitealgorithm A˜(k) and A˜(k) , concluding one iteration for updat- L R in Table II by constructing and solving for two effective ing the entire unit cell. For a pseudocode summary see HamiltoniansH andH insteadofasingleone. C(k−1) C(k) Algorithm 4 in Table III. Thenewlyoptimizedtensorsthengetreplacedinall unit cellsoftheinfinitelattice,andcontributionstotheeffec- tive Hamiltonians have to be recalculated accordingly, 3. Juxtaposition of Both Approaches before moving on to the next site. For a pseudocode summary see Algorithm 3 in Table III. Severalcommentsonthetwopresentedalgorithmsare One could now try to argue, that e.g. in a left to right sweep it is enough at site k to calculate updated A˜(k) in order. First, the parallel algorithm requires substan- C and C˜(k) = C˜(k) only, and to use C˜(k − 1) from tially less computational effort, since the construction of R R the different effective Hamiltonians H can recycle the previous step at site k − 1 as C˜(k)L for calculat- the calculation of the infinite geometricA(sku)mC. Therefore, ing A˜(k) . This approach however fails, as the effective R updatinganentireunitcellonlyrequirestoevaluatetwo HamiltonianusedforcalculatingA˜(k) alreadycontains C infinitegeometricsumsand2N effectiveeigenvalueprob- updated A˜(k − 1)L/R, while the effective Hamiltonian lems. In the sequential algorithm, updating the envi- used for calculating C˜(k−1) does not, and we cannot ronment after every tensor update requires to reevaluate R 10 Algorithm 3 sequential variational uMPS algorithm for multi-site unit cells Input: Hamiltonian H, initial uMPS {A },{A },{C} of an N-site unit cell, convergence threshold (cid:15) L R Output: uMPS approximation {A },{A },{C} of ground state of H, fulfilling fixed point relations (27a), (27b) and (27c) L R up to precision (cid:15). 1: procedure VUMPSMultiSequential(H,{AL},{AR},{C},(cid:15)) 2: initialize current precision (cid:15)prec >(cid:15) 3: while (cid:15)prec >(cid:15) do 4: for n=1,...,N do 5: (optional) Dynamically adjust bond dimension following Appendix B 6: Calculate explicit terms of effective Hamiltonians from a multi-site version H ,H ,H ←HeffTermsMulti(H,{A },{A },{L},{R},(cid:15) ≤(cid:15) ) of Algorithm 1, 5 or 6 A(n)C C(n−1) C(n) L R S prec 7: eCiaglecnusloaltveerg,rocaulnlidngstAatpepAl˜yCHoAfCeff(Cec,tHive Ham)iflrtoomniaAnlgHorAit(nh)mC 1to, 5proerci6sion (cid:15)H <(cid:15)prec using an iterative A(n)C 8: Calculate ground state C˜L of effective Hamiltonian HC(n−1) to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 (cid:46) To ensure gauge consistency, C(n−1) employ a phase convention for C˜ L 9: Calculate ground state C˜R of effective Hamiltonian HC(n) to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 (cid:46) To ensure gauge consistency, C(n) employ a phase convention for C˜ R 10: Calculate new A˜L from A˜C and C˜R using (20) or (22), depending on singular values of C˜R 11: Calculate new A˜R from A˜C and C˜L using (20) or (22), depending on singular values ofC˜L 12: Evaluate new (cid:15)L(n) and (cid:15)R(n) from (18a) and (18b) 13: Replace A(n)L ←A˜L, A(n)R ←A˜R, C(n−1)←C˜L and C(n)←C˜R 14: end for 15: Set (cid:15)prec ←max({(cid:15)L},{(cid:15)R}) 16: (optional) Calculate current expectation values 17: end while 18: return {AL},{AR},{C} 19: end procedure Algorithm 4 parallel variational uMPS algorithm for multi-site unit cells Input: Hamiltonian H, initial uMPS {A },{A },{C} of an N-site unit cell, convergence threshold (cid:15) L R Output: uMPS approximation {A },{A },{C} of ground state of H, fulfilling fixed point relations (27a), (27b) and (27c) L R up to precision (cid:15). 1: procedure VUMPSMultiParallel(H,{AL},{AR},{C},(cid:15)) 2: initialize current precision (cid:15)prec >(cid:15) 3: while (cid:15)prec >(cid:15) do 4: (optional) Dynamically adjust bond dimension following Appendix B 5: for n=1,...,N do 6: Calculate explicit terms of effective Hamiltonians from a multi-site version H ,H ←HeffTermsMulti(H,{A },{A },{L},{R},(cid:15) ≤(cid:15) ) of Algorithm 1, 5 or 6 A(n)C C(n) L R S prec 7: eCiaglecnusloaltveerg,rocaulnlidngstAatpepAl˜y(nH)ACCo(fCe,ffHective H) farmomiltoAnligaonriHthAm(n1)C, 5toorp6recision (cid:15)H <(cid:15)prec using an iterative A(n)C 8: Calculate ground state C˜(n) of effective Hamiltonian HC(n) to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 C(n−1) 9: end for 10: for n=1,...,N do 11: Calculate new A˜(n)L from A˜(n)C and C˜(n) using (20) or (22), depending on singular values of C˜(n) 12: Calculate new A˜(n)R from A˜(n)C and C˜(n−1) using (20) or (22), depending on singular values of C˜(n−1) 13: Evaluate new (cid:15)L(n) and (cid:15)R(n) from (18a) and (18b) 14: end for 15: Replace {AL}←{A˜L}, {AR}←{A˜R} and {C}←{C˜} 16: (optional) Calculate current expectation values 17: Set (cid:15)prec ←max({(cid:15)L},{(cid:15)R}) 18: end while 19: return {AL},{AR},{C} 20: end procedure Table III. Pseudocode for the two approaches for a multi-site unit cell implementation described in Sec. IIE. Algorithm 3 sweeps through the unit cell and sequentially updates tensors site by site, replacing updated tensors in all unit cells before moving on to the next site. Algorithm 4 updates the entire unit cell at once by independently updating tensors on each site.

ebook img

Variational optimization algorithms for uniform matrix product states PDF

1.6 MB·
by  V. Zauner-Stauber
#journals #arxiv
Save to my drive
Quick download
Download
Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Variational optimization algorithms for uniform matrix product states

Variational optimization algorithms for uniform matrix product states V. Zauner-Stauber,1 L. Vanderstraeten,2 M.T. Fishman,3 F. Verstraete,1,2 and J. Haegeman2 1Vienna Center for Quantum Technology, University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria 2Ghent University, Faculty of Physics, Krijgslaan 281, 9000 Gent, Belgium 3Institute for Quantum Information and Matter, California Institute of Technology, Pasadena, California 91125, USA 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 di- mensional 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 inconvergencespeedand precision. We alsodemonstrate thatVUMPSworksveryefficiently forHamiltonianswithlongrangeinteractions. Thenewalgorithmcanbeconvenientlyimplemented as an extension of an already existing DMRG implementation. I. INTRODUCTION Using a translation invariant parameterization gives rise 7 to an energy expectation value with a highly non-linear 1 dependence on the tensor(s). Two different algorithms Thestrategyofrenormalizationgroup(RG)techniques 0 are widely used to obtain such an MPS in the thermo- 2 to successively reduce a large number of microscopic de- dynamic limit. Infinite system DMRG (IDMRG)5,6,24 grees of freedom to a smaller set of effective degrees of n proceeds by performing regular DMRG ona successively freedom has led to powerful numerical and analytical a growing lattice, inserting and optimizing over new ten- J methods to probe and understand the effective macro- sors in the center of the lattice in each step only, effec- scopic behavior of both classical and quantum many 4 tivelymimickinganinfinitelatticebyusingafinite,albeit 2 body systems.1–4 However, it was not until the advent very large lattice. After convergence the most recently of White’s celebrated Density Matrix Renormalization inserted tensors in the center are taken as a unit cell for ] Group (DMRG)5,6 that variational RG methods reached h an infinite MPS approximation of the ground state. An unprecedentedaccuracyinnumericallystudyingstrongly p alternative approach is known as infinite time evolving correlated one dimensional quantum lattice systems at - block decimation (ITEBD).25,26 It works directly in the t low temperature. The underlying variational ansatz of n thermodynamic limit and is based on evolving an initial a Matrix Product States (MPS)7–13 belongs to a class of stateinimaginarytimebyusingaTrotterdecomposition u ansatzes known as Tensor Network States.11,14,15 These of the evolution operator. q variational classes encode the many body wavefunction [ We present a new variational algorithm, inspired by in terms of virtual entanglement degrees of freedom liv- 1 ing on the boundary and thus satisfy an area law scaling tangentspaceideas,13,27,28 thatcombinestheadvantages ofIDMRGandITEBDandaddressessomeoftheirshort- v of entanglement entropy per construction. As such, they comings. Assuchitisdirectlyformulatedinthethermo- 5 provide a natural parameterization for the physical cor- 3 nerofHilbertspace,wherelowenergystatesofquantum dynamic limit, but at the same time optimizes the state 0 manybodysystemsoughttolivein.16,17 MPSinparticu- bysolvingeffectiveeigenvalueproblems,ratherthanem- 7 ployingimaginarytimeevolution. Wefindthatitleadsto larareespeciallyfitforstudyinggroundstatesofstrongly 0 a significant increase in efficiency in all of our test cases. correlated one dimensional quantum systems with local 1. interactions.18–20 ThefollowingsectionintroducesMPSnotationsanddef- 0 initions and presents our variational algorithm, heuristi- 7 The variational parameters in MPS are contained callymotivatedfromtheperspectiveoffinitesizeDMRG. 1 withinlocaltensorsassociatedwiththeindividualsitesof Sec. III illustrates the performance of our algorithm on : the lattice system. For homogeneous systems, the global v varioustestcases,andcomparestoconventionalIDMRG wave function can then be captured using just a single i and ITEBD results. After the conclusion in Sec. IV, we X (or a small number of) such tensors, independent of the provide further technical details in the appendices. Ap- r system size. They consequently offer very natural access pendixAcontainsadditionaltheoreticalbackground: we a tothethermodynamiclimit, providingaclearadvantage derivetheself-consistentconditionsthatcharacterizethe overothernumericalapproachessuchasExactDiagonal- variational minimum and provide additional motivation ization or Quantum Monte Carlo. for our algorithm from the perspective of the MPS tan- On finite lattices, (one-site) DMRG implements the gent space. Appendix B presents a suitable strategy to variational principle (energy minimization) by exploit- expandthebonddimensionoftranslationinvariantMPS. ing that the quantum state is a multilinear function of AppendixCexplainshowtoconstructeffectiveHamilto- the local tensors. By fixing all but one tensors, the nians in the thermodynamic limit. These involve infinite global eigenvalue problem is transformed into an effec- geometric sums of the transfer matrix, which are further tive eigenvalue problem for the local tensor.5,6,12,21–23 studied in Appendix D. 2 II. A VARIATIONAL ALGORITHM FOR A. Uniform MPS MATRIX PRODUCT STATES IN THE THERMODYNAMIC LIMIT A uniform MPS (uMPS) of bond dimension D defined on an infinite translation invariant lattice is parameter- In this section we introduce a variational algorithm ized by a single collection of d matrices As ∈ CD×D for foroptimizingMPSdirectlyinthethermodynamiclimit. s=1,...,d. Theoveralltranslationinvariantvariational Because the algorithm strongly resembles conventional state is then given by DMRG,weexplainitbydescribingasingleiterationstep (cid:88)(cid:16) (cid:17) from the viewpoint of DMRG and show that only a few |Ψ(A)(cid:105)= ...Asn−1AsnAsn+1... |s(cid:105) (1) additional ingredients are needed to arrive at our vari- s ational algorithm. We only briefly motivate these extra and can be represented diagrammatically as ingredients for the sake of readability, and refer to Ap- pendix A for additional explanations and more rigorous |Ψ(A)(cid:105)=... A A A A A ... theoretical motivations. As such, the new algorithm can easily be implemented as an extension to an already ex- isting (I)DMRG implementation. Exploiting the invariance of (1) under local gauge transformationsAs →XAsX−1,withX ∈CD×D invert- We start by considering a setting familiar from con- ible,thestatecanbecastintocertainfavorablerepresen- ventional DMRG: a finite homogeneous one dimensional tations, among them the left and right canonical repre- quantum lattice system, where every site corresponds to sentation a d level system. We label the sites by an integer n and thushaveabasis{|s(cid:105)n,s=1,...,d}forthelocalHilbert (cid:88)AsL†AsL =11 (cid:88)AsLRAsL† =R (2a) spaceonsiten. ThetotalHilbertspaceisspannedbythe (cid:78) s s product basis |s(cid:105) = |s(cid:105) . We assume the dynamics (cid:88) (cid:88) n n AsAs† =11 As†LAs =L, (2b) of the system to be governed by a translation invariant R R R R Hamiltonian H. s s We further consider a variational parameterization of or diagrammatically a ground state approximation of the system, for now in terms of a finite size (site dependent) MPS, but we AL AL willultimatelybeinterestedinthethermodynamiclimit. DMRGproceedstofindthebestvariationalgroundstate = R = R approximationbyemployinganalternatingleastsquares minimization: It starts from some initial state and suc- A¯L A¯L cessively optimizes each of the individual MPS tensors site by site by solving effective (Hamiltonian) eigen- AR AR value problems, in a sweeping process through the lat- = L = L tice until convergence, where each iteration depends on already optimized tensors from previous iterations (see e.g. Refs 5, 6, 12, 21, and 23). A¯R A¯R We are now however interested in the thermodynamic Here L and R correspond to the left and right reduced limit n ∈ Z (but will ignore the technical complications densitymatricesofabipartitionofthestaterespectively. involving a rigorous definition of a Hilbert space in that WehenceforthrefertoA (A )asaleft(right)isometric L R limit). In that case the MPS ground state approxima- tensor, or just a left (right) isometry. tion will be given in terms of a translation invariant uni- Defining the left and right transfer matrices form MPS, described by a single MPS tensor (or a unit cell of N tensors), repeated on all sites. Two immediate T =(cid:88)A¯s ⊗As (3) L/R L/R L/R difficulties arise: Firstly, conventional DMRG updates s the variational state site by site, thus breaking transla- and using the notation (x| and |x) for vectorizations of tion invariance. Secondly, the effective Hamiltonian for a D×D matrix x in the D2 dimensional “double layer” a single-site optimization has to be constructed from an virtual space the transfer matrices act upon, the gauge infinite environment. conditions (2) are equivalent to After briefly introducing the variational class of uni- form MPS and introducing necessary notation and con- (11|T =(11| T |R)=|R) (4a) L L ventions (for further details see Sec. A2), we describe T |11)=|11) (L|T =(L|, (4b) R R how the new algorithm modifies DMRG accordingly to exactly account for these two issues in order to arrive at i.e. 11 and R are the left and right dominant eigenvectors a variational ground state algorithm directly formulated of T , while L and 11 are the left and right dominant L in the thermodynamic limit. eigenvectors of T . R 3 AsisstandardpracticeinDMRG,wemixbothofthese states representations and cast the state into the mixed canon- |Ψ(α,s,β)(cid:105)=|Ψα(cid:105)|s(cid:105)|Ψβ(cid:105) (8a) ical representation AC L R |Ψ(α,β)(cid:105)=|Ψα(cid:105)|Ψβ(cid:105). (8b) C L R |Ψ(A)(cid:105)=(cid:88)(...Asn−1AsnAsn+1...)|s(cid:105) (5a) L C R s B. Effective Hamiltonian =(cid:88)(...Asn−1AsnCAsn+1Asn+2...)|s(cid:105), (5b) L L R R The use of the mixed canonical representation (5a) in s DMRG is of significant importance for the stability, as it reduces the minimization of the (global) energy ex- or diagrammatically, pectationvalue(cid:104)Ψ|H|Ψ(cid:105)/(cid:104)Ψ|Ψ(cid:105)withrespecttoA into C a standard (hermitian) eigenvalue problem, instead of a generalizedone. TheeffectiveHamiltonianforthiseigen- value problem is the system Hamiltonian H projected |Ψ(A)(cid:105)=... AL AL AC AR AR ... onto the degrees of freedom of AC, and is known as the “reduced” or “superblock” Hamiltonian in DMRG. We define the thermodynamic limit version of this re- =... AL AL AL C AR AR ... duced single-site Hamiltonian acting on AC as H (α(cid:48),s(cid:48),β(cid:48)) =(cid:104)Ψ(α(cid:48),s(cid:48),β(cid:48))|H|Ψ(α,s,β)(cid:105) (9) AC(α,s,β) AC AC Here we have defined the center site tensor As (known AL AL AR AR C as the single-site wave function Ψs in DMRG) =··· H ··· AsC =AsLC =CAsR A¯L A¯L A¯R A¯R (6) AC = AL C = C AR Additionally, we also define an effective Hamiltonian acting on the bond matrix C as H (α(cid:48),β(cid:48)) =(cid:104)Ψ(α(cid:48),β(cid:48))|H|Ψ(α,β)(cid:105) (10) C(α,β) C C intermsofthebondmatrixC,whichconstitutesthe(in- vertible) gauge transformation relating A and A via AL AL AR AR L R As = CAsC−1. The singular values of C then encode thLe entangRlement spectrum of the state. Indeed, using =··· H ··· AsC = CAs we can verify that the left and right re- duLced densitRy matrices in (2) are given by L=C†C and A¯L A¯L A¯R A¯R R = CC†. Furthermore, AsC = CAs ensures that (5a) L R which does not appear directly in the context of DMRG, and (5b) are translation invariant and that A and C C but will be needed later for a consistent update of the can be shifted around arbitrarily. Normalization of the state without breaking translation invariance. It can be state,aswellasofthereduceddensitymatricesLandR, interpreted as a “zero site” effective Hamiltonian, which correspondstothesinglecondition(cid:107)C(cid:107)2 =Tr(CC†)=1. 2 would feature in an optimization of the global energy Foreaseofnotationwefurtherintroducethefollowing expectation value with respect to the Schmidt values. partial states In an efficient implementation, these effective eigen- value problems are typically solved using an iterative eigensolver,sothatweonlyneedtoimplementtheaction |ΨαL(cid:105)=(cid:88)(...AsLn−1AsLn)α|...sn−1sn(cid:105) (7a) of HAC and HC onto AC and C. While the energy expectation value is extensive and s thus divergent in the thermodynamic limit, the effective =... AL AL AL α Hamiltonians H and H are well defined and finite in AC C the thermodynamic limit if one properly subtracts the |Ψα(cid:105)=(cid:88)(AsnAsn+1...) |s s ...(cid:105) (7b) current energy expectation value from the Hamiltonian R R R α n n+1 s H. Wedemonstratethisprocedureforthecaseofnearest (cid:80) = α AR AR AR ... neighborinteractionsH = nhn,n+1,wherethetwo-site Hamiltonian h acts on neighboring sites n,n + 1 n,n+1 only. We refer to Appendix C for the case of long range interactions and for general Hamiltonians given in terms withnarbitrary,andusethemtodefinethereducedbasis of Matrix Product Operators (MPOs). 4 Inthecaseofnearestneighborinteractions, theactionofH ontoA splitsupintofourindividualcontributions, AC C which follow from the decomposition |Ψ(cid:105)=(cid:80) As |Ψα(cid:105)|s(cid:105)|Ψβ(cid:105) (left block containing sites n<0, center site α,β,s C,(α,β) L R n=0, and right block containing sites n>0). The action of H onto A is given by AC C A(cid:48)s =(cid:88)htsAt†AkA(cid:96) +hstAkA(cid:96) At †+H As +AsH C k(cid:96) L L C k(cid:96) C R R L C C R tk(cid:96) AL AC AC AR AC AC (11) A(cid:48)C = h + h + HL + HR A¯L A¯R where the first two terms correspond to the Hamiltonian terms h and h coupling the center site to the left and −1,0 0,1 rightblock,respectively,andH andH sumupthecontributionsofalltheHamiltoniantermsh actingstrictly L R n,n+1 to the left and to the right of the center site. TheenvironmentsH andH areusuallyconstructed where (H | can be presented diagrammatically as L R L iteratively while sweeping through the (finite) lattice in   conventional DMRG, or grown successively in every it- AL AL AL eration in IDMRG. In the thermodynamic limit, these     terms consist of a diverging number of individual local HL = hL 11+ + + ... interaction contributions h , and care needs to be   n,n+1 taken in their construction. A¯L A¯L A¯L Indeed, the kth contribution to (H | comes from the L and likewise for |H ). Hamiltonian term h and is given by (h |[T ]k−1. R −k−1,−k L L The transfer matrix T has a dominant eigenvalue of Likewise, [T ]k−1|h ) is the kth contribution to |H ) L R R R magnitude one, with corresponding left and right eigen- stemming from h . Here, we have used the defini- k,k+1 vectors (11| and |R). The projection (h |[T ]k|R) = tions L L (h |R) is the energy density expectation value e = L hL =(cid:88)hskt(cid:96)AtL†AsL†AkLA(cid:96)L (cid:104)tΨhe|he−nke−rg1,y−kh˜|Ψ=(cid:105) ahn−d eis11infrdoempetnhdeenHtaomfilkt.oniSaunb,twraectcinang stk(cid:96) (12) write (h | = (h˜ | + e(11|. The second term is exactly h =(cid:88)hstAkA(cid:96) At †As†, proportioLnal to tLhe left eigenvector of eigenvalue 1 and R k(cid:96) R R R R stk(cid:96) therefore gives rise to a diverging contribution in the ge- ometric sum, corresponding to the total energy of the or diagrammatically left half infinite block. Since this contribution acts as the identity in the effective Hamiltonian H [Eq. (11)], AC wecanhoweversafelydiscardthisdivergingcontribution AL AL without changing the eigenvectors of H . This corre- AC sponds to an overall energy shift of the left half infinite hL = h blocksuchthat(H |R)=0. Fortheremainingpart(h˜ | L L A¯L A¯L the geometric sum converges. With |h˜R) = |hR)−e|11) the same comments apply to the construction of |H ). R AR AR We can evaluate HL and HR recursively as (H[n+1]|=(H[n]|T +(h˜ | hR = h L L L L (14) |H[n+1])=T |H[n])+|h˜ ) R R R R A¯R A¯R with initialization (H[0]| = (h˜ | and |H[0]) = |h˜ ). We L L R R Summing up all such local contributions gives rise to in- can repeat these recursions until e.g. (cid:107)H[n+1] −H[n] (cid:107) L/R L/R finite geometric sums of the transfer matrices TL/R drops below some desired accuracy (cid:15)S. This strategy is conceptually simple and closely resembles the succes- ∞ ∞ sive construction of the environments in the context of (cid:88) (cid:88) (HL|=(hL| [TL]k |HR)= [TR]k|hR), (13) (I)DMRG, but is not very efficient, as its performance is k=0 k=0 comparable to that of a power method. 5 Algorithm 1 Explicit terms of effective Hamiltonians with nearest neighbor interactions and their application Input: two-site Hamiltonian h, current uMPS tensors A , A in left and right gauge, left dominant eigenvector (L| of T , L R R right dominant eigenvector |R) of T , desired precision (cid:15) for terms involving infinite geometric sums L S Output: Explicit terms of effective Hamiltonians H and H , updated A(cid:48) and C(cid:48) AC C C 1: function HeffTerms(H =h,AL,AR,L,R,(cid:15)S) (cid:46) Calculates explicit terms of effective Hamiltonians 2: Calculate hL and hR from (12) 3: Calculate HL and HR by iteratively solving (14) or (preferably) (15), to precision (cid:15)S 4: HAC ←{h,AL,AR,HL,HR} 5: HC ←{h,AL,AR,HL,HR} 6: return HAC,HC 7: end function 8: function ApplyHAC(AC,HAC) (cid:46) Terms of HAC from HeffTerms(H,AL,AR,L,R,(cid:15)S) 9: Calculate updated A(cid:48) from (11) C 10: return A(cid:48) C 11: end function 12: function ApplyHC(C,HC) (cid:46) Terms of HC from HeffTerms(H,AL,AR,L,R,(cid:15)S) 13: Calculate updated C(cid:48) from (16) 14: return C(cid:48) 15: end function Table I. Pseudocode for obtaining the explicit terms of the effective Hamiltonians H and H for systems with nearest AC C neighbor interactions and their applications onto a state. A more efficient approach is to formally perform the for (H | and |H ) to precision (cid:15) , as explained in Ap- L R S geometricsumsin(13)explicitly, andtoiterativelysolve pendix D. the resulting two systems of equations So far, we have discussed the action of H . The ac- AC tionofH ontoC followssimplyfrom(11)byprojecting C (HL|[11−TL+|R)(11|]=(hL|−(hL|R)(11| ontoALorAR. UsingthedefiningpropertyofHLorHR, (15) the result simplifies to [11−T +|11)(L|]|H )=|h )−|11)(L|h ) R R R R C(cid:48) =(cid:88)hstAs†AkCA(cid:96) At †+H C+CH . k(cid:96) L L R R L R stk(cid:96) AL C AR C C (16) C(cid:48) = h + HL + HR A¯L A¯R The first two terms of (11) can be applied in O(d4D3) C. Updating the state operations29, and the last two terms in O(dD3) opera- tions. For (16) the first term can be applied in O(d4D3) InDMRG,wewouldupdatethestatebyreplacingA operations, and the last two terms in O(D3) operations. C with the lowest eigenvector A˜ of H and then shift Togeneratethenecessarytermsfor(11)and(16)wehave C AC the center site to the right by computing an orthogonal to iteratively evaluate two infinite geometric sums, in- factorization A˜s = A˜sC˜ , or to the left by computing volving O(D3) operations (when iteratively solving (15) C L R A˜s = C˜ A˜s. As such, the state gets updated by only the solutions from the previous iteration can be used as C L R startingvectorstospeedupconvergence). Apseudocode replacing the current site with A˜sL or A˜sR, while leaving summary for obtaining the necessary explicit terms of allothersitesuntouched. However,applyingthisscheme H and H and their applications onto a state is pre- in our setting would immediately destroy translation in- seAntCed in TaCble I. variance after a single step. We want to construct an alternative scheme that ap- plies global updates in order to preserve translation in- variance at any time. Global updates can most easily 6 be applied with an explicit uniform parameterization in where matrices P are hermitian and positive. Alter- termsofasingletensorA. Ontheotherhand,DMRGex- native isometric decompositions might be considered in perienceteachesusthatthestabilityisgreatlyenhanced Eq. (21), though it is important that they are unique when applying updates at the level of A and C, which (e.g. QR with positive diagonal in R) in order to have C are isometrically related to the full state. P[(cid:96)/r] ≈P[(cid:96)/r] close to convergence. We therefore calculate the lowest eigenvector A˜ of AC C C H like in DMRG, but additionally also the lowest AC eigenvector C˜ of HC. We then globally update the D. The Algorithm: VUMPS state by finding new A˜ and A˜ as the left and right L R isometric tensors that minimize (cid:80) (cid:107)A˜sC˜ − A˜s(cid:107)2 and s L C Wearenowreadytoformulateourvariationaluniform (cid:80)s(cid:107)C˜A˜sR−A˜sC(cid:107)2 respectively. Theseminimizationprob- MPS (VUMPS) algorithm. As shown in Appendix A, a lems can be solved directly (not iteratively) and without variational minimum (vanishing energy gradient) in the invertingC˜ (seebelow). AsshowninAppendixA,atthe manifoldofuMPSischaracterizedbytensorsA , C and L variational optimum the values of these objective func- A satisfying the conditions R tions go to zero, and current A and C are the lowest C eigenvectors of H and H respectively. H A =E A (23a) AC C AC C AC C For the remainder of this section we omit tildes and H C =E C (23b) C C use the following matricizations of the 3-index tensors As =AsC =CAs. (23c) C L R A =As L,(sα,β) L,(α,β) HereboldsymbolsdenotevectorizationsoftheMPSten- AR,(α,sβ) =AsR,(α,β) (17) sorsandmatricizationsoftheeffectiveHamiltonians,and A[C(cid:96)],(sα,β) =AC[r,](α,sβ) =AsC,(α,β). EHAamCilatnodniaEnCs.3a1re the lowest eigenvalues of the effective We thus want to extract updated AL and AR from up- When iterating the steps outlined in the previous sub- dated AC and C by solving sections, convergence is obtained when these conditions are satisfied. In particular, starting with a properly or- (cid:15) = min (cid:107)A[(cid:96)]−A C(cid:107) (18a) L C L 2 thogonalized initial trial state |Ψ(A)(cid:105) of some bond di- A†LAL=11 mensionD,webeginbysolvingthetwoeigenvalueprob- (cid:15) = min (cid:107)A[r]−CA (cid:107) . (18b) lems for the effective Hamiltonians H and H . Since R ARA†R=11 C R 2 we are still far from the fixed point, thAeCresultinCg lowest energy states A˜ and C˜ will in general not satisfy the In exact arithmetic, the solution of these minimization C gauge condition (23c) together with current A . problems is known, namely A will be the isometry in L/R L the polar decomposition of A[l]C† (and similar for A , Following the procedure of the previous section we C R can however find optimal approximations A˜s and A˜s see Thm. IX.7.2 in Ref. 30). Computing the singular L R for (23c) to arrive at an updated uMPS. Conversely, A˜ value decompositions (SVD) C and C˜ will not be the correct lowest energy eigenstates A[C(cid:96)]C† =U[(cid:96)]Σ[(cid:96)]V[(cid:96)]† C†AC[r] =U[r]Σ[r]V[r]† (19) of the new effective Hamiltonians HA˜C and HC˜ gener- ated from A˜ . We then use the updated state and we thus obtain L/R reiterate this process of alternately solving the effective A =U[(cid:96)]V[(cid:96)]† A =U[r]V[r]†. (20) eigenvalueproblems,andfindingoptimalapproximations L R for A and A to update the state. For a pseudocode L R Notice that close to (or at) an exact solution As = summary of this algorithm, see Table II. C AsC =CAs, the singular values contained in Σ[(cid:96)/r] are We now elaborate on the various steps in the VUMPS L R the square of the singular values of C, and might well algorithm. Firstly, extracting new A˜ from updated L/R fallbelowmachineprecision. Consequently,infinitepre- A˜ and C˜ can be done using the theoretically opti- C cisionarithmetic,correspondingsingularvectorswillnot mal (but numerically often inaccurate) Eq. (20) or the be accurately computed. morerobustEq.(22),dependingonthemagnitudeofthe An alternative that has proven to be robust and still smallest singular value in C˜. As a good uMPS approxi- close to optimal is given by directly using the following mationwillalwaysinvolvesmallsingularvalues,Eq.(22) left and right polar decompositions is preferable most of the time, except maybe during the A[(cid:96)] =U[(cid:96)] P[(cid:96)] C =U[(cid:96)]P[(cid:96)] (21a) first few iterations. C AC AC C C The maximum of the error quantities (18) A[r] =P[r]U[r] C =P[r]U[r] (21b) C AC AC C C (cid:15) =max((cid:15) ,(cid:15) ) (24) prec L R to obtain A =U[(cid:96)] U[(cid:96)]† A =U[r]†U[r], (22) providesanerrormeasureforthefixedpointconditionin L AC C R C AC Eq. (23c) and is used as a global convergence criterion. 7 Algorithm 2 variational uMPS algorithm for single-site unit cells Input: Hamiltonian H, initial uMPS A , A , C, convergence threshold (cid:15) L R Output: uMPS approximation A , A , C of ground state of H, fulfilling fixed point relations (23a), (23b) and (23c) up to L R precision (cid:15) 1: procedure VUMPS(H,AL,AR,C,(cid:15)) 2: initialize current precision (cid:15)prec >(cid:15) 3: while (cid:15)prec >(cid:15) do 4: (optional) Dynamically adjust bond dimension following Appendix B 5: Calculate explicit terms of effective Hamiltonians HAC,HC ←HeffTerms(H,AL,AR,L,R,(cid:15)S ≤(cid:15)prec) from Algorithm 1, 5 or 6 6: Calculate ground state A˜C of effective Hamiltonian HAC to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHAC(A ,H ) from Algorithm 1, 5 or 6 C AC 7: Calculate ground state C˜ of effective Hamiltonian HC to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 C 8: Calculate new A˜L and A˜R from A˜C and C˜ using (20) or (22), depending on singular values of C˜ 9: Evaluate new (cid:15)L and (cid:15)R from (18) 10: (optional) Calculate current expectation values 11: Set (cid:15)prec ←max((cid:15)L,(cid:15)R) and replace AL ←A˜L, AR ←A˜R and C ←C˜ 12: end while 13: return AL, AR, C 14: end procedure Table II. Pseudocode of the VUMPS algorithm described in Sec. IID. Terms within step 5 involving the evaluation of infinite geometric sums usually require the left dominant eigenvector L of T and the right dominant eigenvector R of T , for which R L L = C†C and R = CC† with current C are a good enough approximation to current precision (cid:15) (see main text). Notice prec that this algorithm is free of any possibly ill-conditioned inverses and therefore has no convergence issues in the presence of small Schmidt values. It also does not require expensive reorthogonalizations of the state at intermediate iterations. It measures the precision of the current uMPS ground prohibitive. Whenthisapproachconverges,theresulting state approximation. Within every iteration, we use it- fixedpointtensorsinthecentercanbeassumedtospecify erative methods (e.g. some variation of Lanczos) to find theunitcellofaninfiniteMPS.VUMPShastheimmedi- the eigenvectors A˜ and C˜ of the Hermitian operators ate advantage that i) it directly works in the thermody- C H and H . As the goal is to drive the state towards namiclimitatalliterationsandii)itcompletelyreplaces AC C the fixed point relations in Eqs. (23a) and (23b), it is the entire state after every iteration, thus moving faster not necessary to solve these eigenvalue problems to full through the variational manifold. In contrast, IDMRG machine precision. Rather, it is sufficient to use a toler- keeps memory of earlier iterations and cannot guarantee ance (cid:15) chosen relative to (cid:15) .32 A value of (cid:15) of the a monotonically decreasing energy that converges to an H prec H order of (cid:15) /100 has proven to work well in practice. It optimum associated with a translation invariant MPS in prec is also worthwhile to use tensors from the previous iter- whichtheeffectsoftheboundaryhavecompletelydisap- ation as initial guess for the iterative solvers to speed up peared. The advantages of VUMPS come with a greater convergence. computationalcostperiteration,astwoeigenvalueprob- Asthemainpartofthealgorithmworksatfixedbond lems(forAC andforC)and–inthecaseofnearestneigh- dimension(i.e.itisasingle-siteschemeinDMRGtermi- bor interactions – two linear systems (for HL and HR) nology),onemightchoosetoincreasethebonddimension havetobesolved. IDMRGonlysolvesasingleeigenvalue D before starting a new iteration. We have developed a problem and builds HL and HR step by step in every it- subspace expansion technique that works directly in the eration. The latter approach is analogous to a power thermodynamic limit and is explained in Appendix B. method for eigenvalue problems and, while very cheap, While the true comparison of this algorithm with is expected to require many iteration steps to converge, IDMRG5,24 and ITEBD26 will take place in Sec. III by especially for systems with large correlation lengths (e.g. gathering actual numerical simulation results, we can al- close to criticality). ready compare the theoretical properties of these algo- ITEBD26 isbasedonevolvinganinitialstateinimagi- rithms. Neither IDMRG or ITEBD is truly solving the narytimebyusingaTrotterdecompositionoftheevolu- variationalprobleminthesenseofdirectlytryingtosat- tion operator. Like VUMPS, ITEBD works in the ther- isfy the fixed point conditions Eqs. (23). IDMRG closely modynamiclimitatanyintermediatestep,typicallywith resembles regular DMRG on a successively growing lat- a unit cell that depends on how the Hamiltonian was tice, as it inserts and optimizes over new tensors in the split into local terms in order to apply the Trotter de- center of the lattice in each step. Tensors from previ- composition. Furthermore, as every application of the ous steps are not updated, as this would render the cost evolution operator increases the virtual dimension of the 8 MPS, truncation steps are required to restore the origi- algorithmstendtoconvergetomaximallysymmetrybro- nal(oranysuitable)valueofthebonddimension. While ken states for which the entanglement is minimal. This VUMPScantakebigstepsthroughthevariationalspace, isalsothecaseforVUMPS.Onecancontrolwhichstate timestepsinITEBDhavetobechosensufficientlysmall the algorithm converges to by suitably biasing the ini- (especially in the final steps of the algorithm) to elimi- tial state or by adding small perturbation terms to the natetheTrottererror,whichnegativelyaffectstherateof Hamiltonian which explicitly break the symmetry, and convergence. (Ref.33howeverproposesaschemetoeffec- which are switched off after a few iterations. tively obtain a larger time step). Furthermore, the Trot- Explicit conservation of translation symmetry was the tersplittingessentiallylimitstheapplicabilityofITEBD very first requirement in the construction of VUMPS. to short-range interactions and dictates the size of the In the case of spontaneous breaking of translation sym- unit cell of the resulting MPS, e.g. in the most common metry down to N-site translation symmetry (as e.g. in case of nearest neighbor interactions a two-site unit cell the case of a state with antiferromagnetic order), en- isobtained. (TheapproachofRef.34toobtainatransla- forcing one-site translation symmetry would result in a tioninvariantMPSisrestrictedtocertainHamiltonians, (non-injective) equal weight superposition of all symme- but see Ref. 35 for an alternative proposal that can in try broken uMPS ground state approximations. In order fact also deal with long range interactions.) to reach an optimal accuracy with a given bond dimen- Finally, we can also compare VUMPS to the more sion,suchasuperpositionofN statesishoweverundesir- recent time dependent variational principle (TDVP),27 able, astheeffectivebonddimensionisreducedtoD/N. which was implemented as an alternative approach to In the case where this situation cannot be amended by simulate real and imaginary time evolution within the a simple unitary transformation that restores one-site manifold of MPS by projecting the evolution direction translation symmetry (such as e.g. flipping every second onto the MPS tangent space. This approach can be ap- spin in the case of an antiferromagnet), it is preferable plied to translation invariant MPS, independent of the to choose an MPS ansatz with a N-site unit cell, such type of Hamiltonian. When used to evolve in imaginary that the state can spontaneously break translation sym- time, it can be identified as a covariant formulation of metry. The generalization of the algorithm to multi-site a gradient descent method, in that it evolves the state unit cells is described in the next section. in the direction of the gradient of the energy functional, preconditionedwiththemetricofthemanifold. Assuch, the energy decreases monotonically and at convergence, E. Multi Site Unit Cell Implementations an exact (local) minimum is obtained, as characterized by the vanishing gradient. However, in its original for- We now generalize the VUMPS algorithm of the pre- mulation,TDVPwasnotformulatedinacentersiteform vious section for one-site translation invariant uMPS to and was therefore unstable and restricted to small time the setting of translation invariance over N sites. Such steps. For finite systems, a different formulation of the a uMPS ansatz is then parameterized by N independent TDVP algorithm was provided in Ref. 28, which allows tensorsA(k)s ∈CD×d×D,k =1,...,N,whichdefinethe fortakingthelimitoftheimaginarytimesteptoinfinite, unit cell tensor and then becomes provably equivalent to the single-site DMRGalgorithm. VUMPScanbemotivatedfromthese As =A(1)s1...A(N)sN, (25) developments, as explained in Appendix A. We conclude this section by elaborating on how to in- wheres=(s1,...,sN)isacombinedindex. Wecanthen corporate symmetries in the algorithm. The construc- write the variational state as tuinointaoryf usyMmPmSettrhiaets iiss eeqxupilvicailtelnytintova(rIi)aDntMuRnGd1e2r,2o3nasnitde |Ψ(A)(cid:105)=(cid:88)(...Asn−1AsnAsn+1...)|s(cid:105) (I)TEBD36,37, and it is immediately clear that the var- s ious steps in VUMPS have a corresponding covariant andtheleftandrightorthonormalformsaregivenbythe formulation. The same comments apply to time rever- relations sal symmetry, in which case everything can be imple- mented in real arithmetic, or to reflection symmetry, in (cid:88)A(k)s†A(k)s =11 which case C and As will be symmetric matrices and L L C s AsR = AsLT (which implies that HL and HR are also re- (cid:88)A(k)s R(k)A(k)s† =R(k−1) (26a) lated). In all of these cases, the computational cost is L L s reduced. However, explicitly imposing the symmetry in the MPS requires caution, as the physical system might and have spontaneous symmetry breaking, or – more subtly – might be in a symmetry protected topological phase (cid:88)A(k)s A(k)s† =11 R R where the symmetries cannot be represented trivially on s (26b) the MPS tensor. (cid:88)A(k)s†L(k−1)A(k)s =L(k), R R In the case of spontaneous symmetry breaking, MPS s 9 where it is understood that N +1≡1 and 0≡N. determine A˜(k) from A˜(k) and C˜(k −1) . Rather, R C R Defining the bond matrices C(k) as the gauge trans- C˜(k) has to be recalculated using an updated effective L formation that relates left and right canonical form via Hamiltonian, whichexactly leadsto thesequentialAlgo- C(k−1)A(k)s =A(k)sC(k),wehaveR(k)=C(k)C(k)† rithm 3. R L and L(k)=C(k)†C(k). We can then also cast |Ψ(A)(cid:105) in Thereisanadditionalsubtletythatneedstobeconsid- amixedcanonicalformsimilarto(5a)withthecentersite ered,inorderforalltensorstofulfillthegaugeconstraints tensor given by A(k)s =A(k)sC(k)=C(k−1)A(k)s. (27c) to current precision. Bond matrices C˜(k) are cal- C L R The variational minimum within this set of states is culated as lowest energy eigenvectors of effective Hamil- characterized by the following 3N fixed point relations tonians H and are therefore only determined up to a C(k) phase. Consider C(k) defined between sites k and k+1. HA(k)CA(k)C =EA(k)CA(k)C (27a) At step k it is updated as C˜(k)R and used to calculate HC(k)C(k)=EC(k)C(k) (27b) A˜(k)s. Inthenextstepk+1howeveritisrecalculatedas L A(k)s =A(k)sC(k)=C(k−1)A(k)s. (27c) C˜(k+1) (with an updated effective Hamiltonian) and C L R L usedtodetermineA˜(k+1)s. Atthefixedpointweshould R Notice that due to (27c), the relations for different k are then have C˜(k) = C˜(k+1) = C(k), but this is only connected. There are several possible strategies for con- R L trueifthereisnophaseambiguity,whichwouldalsocon- structing algorithms which obtain states satisfying these sequently lead to a phase mismatch between A˜(k) and conditions. L C˜(k)afterstepk+1. Thisissuedoesnotposeaproblem Inthefollowingwepresenttwoapproacheswhichhave for algorithm convergence (during calculations, matrices shown good performance and stable convergence, which C(k) always appear as products of the form C(k)†C(k) we shall term the “sequential” and “parallel” methods. or C(k)C(k)† and mismatching phases thus cancel out), But let us first elaborate on computing effective Hamil- but can be easily circumvented by employing a phase tonians for multi-site unit cells, which works similarly in convention when calculating updated C˜(k). both methods. We again restrict to the case of nearest neighbor interactions, such that the effective Hamiltoni- ans are constructed similar as in Sec. IIB. To construct 2. Parallel Algorithm e.g. the left block Hamiltonian H , we first collect all L local contributions from a single unit cell in h , before L performing the geometric series of the transfer matrix, Intheparallelapproach,wechoosetoupdateanentire whichnowmediatesatranslationoveranentireunitcell. unit cell at once, using effective Hamiltonians generated from the same current state. To that end, we first gen- erate all terms necessary for all H and H . For A(k)C C(k) 1. Sequential Algorithm the case of nearest neighbor interactions, the contribu- tions H and H to the left and right environment out- L R The sequential algorithm is inspired by finite size side the unit cell can be shared, so that the correspond- DMRG, in that we sweep through the unit cell, succes- ing geometric sum only needs to be computed once, and sively optimizing one tensor at a time while keeping ten- contributions inside the unit cell are obtained through sorsonothersitesfixed. Noticethatatsitekwehowever successive applications of transfer matrices. need two updated bond matrices C˜(k) =C˜(k−1) and Next, we simultaneously and independently solve for C˜(k) = C˜(k), in order to calculate Lupdated A˜(k)s the ground states A˜(k)C and C˜(k) of all 2N effective R L/R Hamiltonians at once. Once these are obtained we again from A˜(k)s ≈ A˜(k)s C˜(k) ≈ C˜(k) A˜(k)s. We thus C L R L R simultaneouslyandindependentlydetermineallupdated havetoamendsteps5,6and7ofthesingle-sitealgorithm A˜(k) and A˜(k) , concluding one iteration for updat- L R in Table II by constructing and solving for two effective ing the entire unit cell. For a pseudocode summary see HamiltoniansH andH insteadofasingleone. C(k−1) C(k) Algorithm 4 in Table III. Thenewlyoptimizedtensorsthengetreplacedinall unit cellsoftheinfinitelattice,andcontributionstotheeffec- tive Hamiltonians have to be recalculated accordingly, 3. Juxtaposition of Both Approaches before moving on to the next site. For a pseudocode summary see Algorithm 3 in Table III. Severalcommentsonthetwopresentedalgorithmsare One could now try to argue, that e.g. in a left to right sweep it is enough at site k to calculate updated A˜(k) in order. First, the parallel algorithm requires substan- C and C˜(k) = C˜(k) only, and to use C˜(k − 1) from tially less computational effort, since the construction of R R the different effective Hamiltonians H can recycle the previous step at site k − 1 as C˜(k)L for calculat- the calculation of the infinite geometricA(sku)mC. Therefore, ing A˜(k) . This approach however fails, as the effective R updatinganentireunitcellonlyrequirestoevaluatetwo HamiltonianusedforcalculatingA˜(k) alreadycontains C infinitegeometricsumsand2N effectiveeigenvalueprob- updated A˜(k − 1)L/R, while the effective Hamiltonian lems. In the sequential algorithm, updating the envi- used for calculating C˜(k−1) does not, and we cannot ronment after every tensor update requires to reevaluate R 10 Algorithm 3 sequential variational uMPS algorithm for multi-site unit cells Input: Hamiltonian H, initial uMPS {A },{A },{C} of an N-site unit cell, convergence threshold (cid:15) L R Output: uMPS approximation {A },{A },{C} of ground state of H, fulfilling fixed point relations (27a), (27b) and (27c) L R up to precision (cid:15). 1: procedure VUMPSMultiSequential(H,{AL},{AR},{C},(cid:15)) 2: initialize current precision (cid:15)prec >(cid:15) 3: while (cid:15)prec >(cid:15) do 4: for n=1,...,N do 5: (optional) Dynamically adjust bond dimension following Appendix B 6: Calculate explicit terms of effective Hamiltonians from a multi-site version H ,H ,H ←HeffTermsMulti(H,{A },{A },{L},{R},(cid:15) ≤(cid:15) ) of Algorithm 1, 5 or 6 A(n)C C(n−1) C(n) L R S prec 7: eCiaglecnusloaltveerg,rocaulnlidngstAatpepAl˜yCHoAfCeff(Cec,tHive Ham)iflrtoomniaAnlgHorAit(nh)mC 1to, 5proerci6sion (cid:15)H <(cid:15)prec using an iterative A(n)C 8: Calculate ground state C˜L of effective Hamiltonian HC(n−1) to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 (cid:46) To ensure gauge consistency, C(n−1) employ a phase convention for C˜ L 9: Calculate ground state C˜R of effective Hamiltonian HC(n) to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 (cid:46) To ensure gauge consistency, C(n) employ a phase convention for C˜ R 10: Calculate new A˜L from A˜C and C˜R using (20) or (22), depending on singular values of C˜R 11: Calculate new A˜R from A˜C and C˜L using (20) or (22), depending on singular values ofC˜L 12: Evaluate new (cid:15)L(n) and (cid:15)R(n) from (18a) and (18b) 13: Replace A(n)L ←A˜L, A(n)R ←A˜R, C(n−1)←C˜L and C(n)←C˜R 14: end for 15: Set (cid:15)prec ←max({(cid:15)L},{(cid:15)R}) 16: (optional) Calculate current expectation values 17: end while 18: return {AL},{AR},{C} 19: end procedure Algorithm 4 parallel variational uMPS algorithm for multi-site unit cells Input: Hamiltonian H, initial uMPS {A },{A },{C} of an N-site unit cell, convergence threshold (cid:15) L R Output: uMPS approximation {A },{A },{C} of ground state of H, fulfilling fixed point relations (27a), (27b) and (27c) L R up to precision (cid:15). 1: procedure VUMPSMultiParallel(H,{AL},{AR},{C},(cid:15)) 2: initialize current precision (cid:15)prec >(cid:15) 3: while (cid:15)prec >(cid:15) do 4: (optional) Dynamically adjust bond dimension following Appendix B 5: for n=1,...,N do 6: Calculate explicit terms of effective Hamiltonians from a multi-site version H ,H ←HeffTermsMulti(H,{A },{A },{L},{R},(cid:15) ≤(cid:15) ) of Algorithm 1, 5 or 6 A(n)C C(n) L R S prec 7: eCiaglecnusloaltveerg,rocaulnlidngstAatpepAl˜y(nH)ACCo(fCe,ffHective H) farmomiltoAnligaonriHthAm(n1)C, 5toorp6recision (cid:15)H <(cid:15)prec using an iterative A(n)C 8: Calculate ground state C˜(n) of effective Hamiltonian HC(n) to precision (cid:15)H <(cid:15)prec using an iterative eigensolver, calling ApplyHC(C,H ) from Algorithm 1, 5 or 6 C(n−1) 9: end for 10: for n=1,...,N do 11: Calculate new A˜(n)L from A˜(n)C and C˜(n) using (20) or (22), depending on singular values of C˜(n) 12: Calculate new A˜(n)R from A˜(n)C and C˜(n−1) using (20) or (22), depending on singular values of C˜(n−1) 13: Evaluate new (cid:15)L(n) and (cid:15)R(n) from (18a) and (18b) 14: end for 15: Replace {AL}←{A˜L}, {AR}←{A˜R} and {C}←{C˜} 16: (optional) Calculate current expectation values 17: Set (cid:15)prec ←max({(cid:15)L},{(cid:15)R}) 18: end while 19: return {AL},{AR},{C} 20: end procedure Table III. Pseudocode for the two approaches for a multi-site unit cell implementation described in Sec. IIE. Algorithm 3 sweeps through the unit cell and sequentially updates tensors site by site, replacing updated tensors in all unit cells before moving on to the next site. Algorithm 4 updates the entire unit cell at once by independently updating tensors on each site.

See more

The list of books you might like

Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.
DMCA & Copyright: Dear all, most of the website is community built, users are uploading hundred of books everyday, which makes really hard for us to identify copyrighted material, please contact us if you want any material removed.
Copyright @ PDFdrive.to, 2022 | We love our users