This page is a collection and summary of some generic improvements
without any claim to completeness. These are essentially tricks which
are relatively independent of the actual time-evolution method and can
mostly be implemented in all of the methods we just reviewed. They are
not strictly necessary to implement in conjunction with any method and
as such will not be benchmarked in detail later, but they are useful
to keep in mind in case of particularly hard or challenging problems.
Combining Heisenberg and Schrödinger picture time evolutionn
In the context of MPS methods, combining Schrödinger and Heisenberg
picture[1][2] time evolution was first proposed in
Ref. 3. Considering a time-dependent observable
⟨ϕ∣O^(t)∣ψ⟩
between two arbitrary states, in the Schrödinger picture we would
evaluate
⟨ϕ∣O^(t)∣ψ⟩=(⟨ϕ∣eitH^)O^(e−itH^∣ψ⟩).
That is, we apply time-evolution operators to the states
∣ϕ⟩ and ∣ψ⟩ to obtain time-evolved states
∣ϕ(t)⟩ and ∣ψ(t)⟩ and then evaluate the
time-independent observable between them. The maximum time t
obtainable is then typically limited by the entanglement growth in the
states, resulting in larger and larger bond dimensions or errors.
In comparison, the Heisenberg picture would see us time evolve the
operator O^ as
⟨ϕ∣O^(t)∣ψ⟩=⟨ϕ∣(eitH^O^e−itH^)∣ψ⟩
while keeping the states ∣ϕ⟩, ∣ψ⟩ static. Again,
the maximal obtainable time t is limited by the entanglement growth in
the operator O^ and the maximal bond dimension we can use to
represent it. (Note that it is difficult to compare errors
between MPS and MPO truncations, as the error of the MPO truncation is
given by the operator norm whereas during an MPO compression, we only
control the 2-norm of the operator[4].) If we now
combine the two evolutions as
we can obtain times t1+t2 while only requiring MPO and MPS bond
dimensions typical of times t1 and t2 respectively. Note that in
this case, the computationally limiting operation is no longer the
time evolution itself but the evaluation of observables given as the
tensor network of a large-w MPO between two large-m
MPS[3].
Complex time steps
Let us assume that we have a time-evolution operator U^′(δ)=1^−iδH^ which
is exact to first order. Applying this operator to a state will result
in an error O(δ2) compared to the exact evolution with the
operator U^(δ)=e−iδH^. Repeating
the process T/δ times to obtain a state at final time T, we
incur an error O(δ2)δT=O(δ).
However, if we allow complex intermediate steps δ1 and
δ2, we can solve
Choosing these values for δ1,2 then results in a third-order
error per time step and a second-order error overall. This choice of
time steps is suggested in particular in combination with the MPO
\WII method to obtain a better error per time step. The cost of the
method only grows linearly with the number of evolution operators and,
e.g., four operators U^′(δ1,2,3,4) are required
for a third-order error [5]. Four steps for the third-order approximation have an algebraic form
The drawback is the loss of unitarity at each individual time step
which may be disadvantageous. Furthermore, if the time evolution is
purely imaginary (e.g., for finite-temperature calculations) and the
Hamiltonian does not contain complex coefficients, one may avoid
complex arithmetic entirely and only use real floating-point scalars
for 50% less memory usage and an approximately four-fold speed-up on
matrix multiplications. Unfortunately, it is then impossible to use
this trick to reduce the time-step error.
Green’s functions 1: Removal of low-lying states
This trick was first proposed in Ref. 7
and is relatively straightforward to implement. Assume a ground state
∣0⟩ as an MPS obtained via DMRG and let us shift the
Hamiltonian such that this state has energy E0=0,
H^∣0⟩=0. We are then interested in the time-dependent
observable
x(t)=⟨0∣A^(t)B^∣0⟩
where A^ and B^ are typically local operators such as
creators or annihilators. The evolution of B^∣0⟩ is
generically non-trivial and if we want to capture all frequency
contributions in x(t), we need to evolve until at least times
t′=E11 where E1 is the energy of the lowest
eigenstate with non-zero energy ∣1⟩ contained in B^∣0⟩. In contrast, to capture contributions of higher-energy
states ∣n⟩ with energies En>E1, we only need to evolve
to shorter times t′′=En1<t′.
However, a few low-lying eigenstates can often be calculated also with
DMRG by orthogonalizing against previously-found eigenstates. Hence if
we run DMRG multiple times, we can obtain not just the ground state
∣0⟩ but also further eigenstates ∣1⟩, ∣2⟩
etc. If we use quantum numbers and B^ changes the quantum number
of the state, these additional eigenstates should be calculated in the
quantum number sector of B^∣0⟩. If we then orthogonalize
B^∣0⟩ against ∣1⟩, ∣2⟩ etc., we remove
the contributions which rotate (in real-time) or decay (in
imaginary-time evolutions) the slowest and hence require the longest
time evolutions. The evolution of the removed states can then be done
exactly as we know both their energy and initial weight in B^∣0⟩. Even missing one of the eigenstates due to convergence
problems with DMRG does not introduce an error but merely decreases
the effectivity of the method.
Green’s functions 2: Linear prediction
Calculating dynamical structure factors or more generally spectral
functions from time-dependent data requires two Fourier
transformations: first, one needs to transform from real-space to
momentum-space and second from real-time to frequency. The former
transformation is typically unproblematic, but the latter
transformation suffers either from overdamping or strong spectral
leakage if the available maximal time t is insufficient. Linear
prediction[8][9][10] assumes that the real-time
momentum-space Green’s function G(k,t) is composed of multiple
distinct exponentially decaying and oscillating contributions arising
from a distinct pole structure of G(k,ω). If this is true for a
time series x1,x2,x3,…, an additional data point
x~n can be approximated well by the form
x~n=−i=1∑paixn−i−1
with suitably chosen coefficients ai independent of n. We hence
first compute a finite time series which is still as long as we can
manage with time-dependent MPS methods. Subsequently, we need to find
coefficients ai such that the above holds for the data we computed
exactly. Using those ai, we can then extend the time series to
arbitrarily long times to generate a sufficiently long time series
that a subsequent Fourier transform only requires minimal damping and
hence provides for clear features.
It is useful to divide the available calculated data into three
segments: first, one should discard an interval
[x0,…,xD−1] at the beginning which captures short-time
physics irrelevant and untypical of the longer-time behavior. Second,
a fitting interval [xD,…,xD+N−1] of N points is
selected over which the coefficients ai are minimized. Third, trust
in the prediction is increased if it coincides with additional
calculated data [xD+N,…,xmax] outside the
fitting interval.
To select the ai, we want to minimize the error
ϵ=k=D∑D+N−1∣x~k−xk∣2.
Note that to evaluate x~k, D must be larger than the
number of coefficients p. The coefficient vector a is obtained as
respectively. Once the ai are obtained, data ideally can be
generated initially for the interval
[xD+N,…,xmax] and, once verified to coincide
with the calculated data, extended to arbitrary times.
Several numerical pitfalls need to be considered here: First, the
matrix R may be singular. Two possible remedies include
addition of a small shift ε or reduction of the number of
parameters p. Ideally the latter should be considered, but may lead
to problems finding the optimal non-singular p. Second, if we
construct the vector
xn=[xn−1,…,xn−p]T
we can move it forward one step as
x~n+1=Axn
where the matrix A is of the form
A=⎝⎛−a110⋮0−a201⋱⋯−a300⋱0⋯⋯⋯⋱1−ap00⋮0⎠⎞
and its eigenvalues αi contain the frequencies and dampings of
the aforementioned oscillations and exponential decays. As such,
αi>1 are unphysical and need to be dealt with, it appears
[10] that setting those contributions to zero works
best.
Purification insertion point (PIP)
(Above) Graphical representation of the expectation value of an out-of-time-ordered correlator within the purification framework of finite-temperature MPS calculations.
Calculating out-of-time-ordered correlators (OTOC) allows us to
measure the scrambling of quantum information and finds many
interesting and current applications. In general an OTOC of operators
W^,V^ is given as an ensemble average
wherein we have suppressed the time-ordered terms and define the OTOC
as the out-of-time ordered part FβV^,W^(t). At
finite temperature, we have to use a purification to evaluate this
quantity. If we would calculate the time evolutions in FβV^,W^(t) naively by direct evolution only in the physical
degrees of freedom we would require O(N2) time steps to
obtain the OTOC at time t=Nδ.
Clearly, the growing numerical expenses forbid to reach both large
system sizes and long time scales t. Graphically representing this
process (see above) immediately suggests to transform the operators in the OTOC in
some way as to evenly distribute the required time evolutions leading
to only linear scaling of effort in time t. In the following, we
will explain how to transform these operators in the purification
picture and alter the purification insertion point (PIP). For related
work in the framework of matrix-product operators,
cf. Ref. 11.
Consider the ensemble average FX^,Y^,Z^,β≡Tr{ρ^(β)Z^Y^X^} for some
global operators X^,Y^,Z^ at inverse temperature
β. Using the cyclic property of the trace the ensemble average
can now be written as expectation value in the enlarged Hilbert space
where we have introduced the purified finite temperature state
∣2β⟩≡ρ^(2β)∣0⟩ based on the infinite
temperature state ∣0⟩. A graphical representation of recasting
the trace into an expectation value is given by the two networks (a)
and (b) in the figure below with out-going indices representing row
vectors and in-going indices column vectors.
(Above) Different choices of operator purifications for ensemble average FX^,Y^,Z^,β=Tr{ρ^(β)Z^Y^X^}.
From the pictographical representation we motivate the infinite temperature state ∣0⟩ to be represented by a rank (2,0) tensor ∣0⟩≡∑a,bˉDa,bˉ∣a⟩∣bˉ⟩ and correspondingly ⟨0∣ by a rank (0,2) tensor ⟨0∣≡∑a,bˉDa,bˉ⟨a∣⟨bˉ∣, where we have placed a bar over those indices labeling ancilla degrees of freedom.
These tensors have to fulfill the orthogonality conditions
bˉ∑Da,bˉDc,bˉ=δca,a∑Da,bˉDa,cˉ=δcˉbˉ
so that the tensor elements can be choosen to be Da,bˉ≡Da,bˉδbˉaˉ and Da,bˉ≡Da,bˉδaˉbˉ.
When contracted over physical degrees of freedom, the action of these tensors is to convert row vectors into column vectors and vice versa
DO^D†=a,c∑bˉ,dˉ∑Da,bˉO^acDc,dˉ∣bˉ⟩⟨dˉ∣=bˉ,dˉ∑O^cˉaˉ∣aˉ⟩⟨cˉ∣=O^t.
If we now interpret indices carrying a bar as maps between ancilla degrees of freedom we can reformulate the purification in terms of the D tensors
FX^,Y^,Z^,β=a,c,…,g,bˉ∑Da,bˉρ^ca(2β)Z^dcY^edX^feρ^gf(2β)Dg,bˉ.
Inserting identities on the physical Hilbert space between ρ^ and X^ as well as X^ and Y^ and making explicit use of the representation of D^ we obtain
FX^,Y^,Z^,β=a,c,d,g,bˉ,eˉ,fˉ,el,er,fl,fr∑Da,bˉρca(2β)Z^dcY^eldδerelDel,eˉDer,eˉX^flerδfrflDfl,fˉDfr,fˉρgfr(2β)Dg,bˉ=a,c,d,bˉ,eˉ,fˉ,el∑Da,bˉρca(2β)Z^dcY^eldDel,eˉact on HAX^eˉfˉρfˉbˉ(2β)=⟨0∣(ρ^t(2β)X^t)A⊗(ρ^(2β)Z^Y^)P∣0⟩
so that now ∑fˉX^eˉfˉρ^fˉbˉ(2β)≡ρ^t(2β)X^t are acting on the ancilla space HA, i.e., we have shifted the purification insertion point.
Again these manipulations can be represented efficiently in a graphical notation and are given as (c) in the figure above.
Using this procedure, we can rewrite the OTOC FβV^,W^(t) as
FβV^,W^(t)=Re[Tr{ρ^(2β)U^†(t)V^†U^(t)W^†U^†(t)V^U^(t)W^ρ^(2β)}]=Re[⟨0∣(U^†(t)V^U^(t)W^ρ^(2β))At⊗(ρ^(2β)U^†(t)V^†U^(t)W^†)P∣0⟩].
Defining the initial states
∣W⟩∣Wβ⟩≡W^P†⊗1^A∣0⟩,≡ρ^P(2β)⊗(W^∗ρ^A(2β))∣0⟩,
and their purified time evolutions
∣W(t)⟩∣Wβ(t)⟩≡U^P∗(t)⊗U^A∗(t)∣W⟩,≡U^P∗(t)⊗U^A∗(t)∣Wβ⟩
the OTOC can be obtained by calculating the expectation value
FβV^,W^(t)=Re[⟨Wβ(t)∣V^P†⊗V^At∣W(t)⟩].
We hence only need N steps to evaluate all expectation values for times t=Nδ.
From a more general point of view shifting the purification insertion
point in the OTOCs reformulates the multiple Schrödinger time
evolutions of the physical system in the canonical choice of the PIP
into a Heisenberg time evolution on both the physical and ancilla
system of a generalized observable
V^P†⊗V^A†t.
Local basis optimization
(Above) Local basis optimization matrices Uj are inserted on the
physical legs of the MPS to transform a large physical basis
(indicated by doubled lines) into a smaller effective basis of the MPS
tensors Mj.
While the dimension of local Hilbert spaces is typically very limited
in spin and electronic systems, bosonic systems potentially require a
large local dimension σ=O(100). As this local
dimension typically enters at least quadratically in some operations
on matrix-product states, some way to dynamically select the most
relevant subspace of the local Hilbert space is potentially extremely
helpful. The local basis
transformation[12][13][14] method
provides for just this: by inserting an additional matrix Uj on
each physical leg of the MPS, the large physical dimension is
transformed into a smaller effective basis. The rank-3 MPS tensor then
only has to work with the smaller effective basis. The method was
adapted for TEBD time evolution in Ref. 14 but it
is also straightforward to use in the other time-evolution methods
presented here. For the MPO \wiii and global Krylov methods, only the
MPO-MPS product has to be adapted to generate additionally a new
optimal local basis after each step. The DMRG, TDVP and local Krylov
method translate[15][16][17] directly in much the same way as
DMRG.
Matrix product simulations of non-equilibrium steady states of quantum spin chains, Tomaz Prosen, Marko Znidaric, Journal of Statistical Mechanics: Theory and Experiment2009, P02035 (2009)