Table Of ContentVariational 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.
