Krylov Approximation

other, previous to each time step, with the help of so called swap gates^[40]. However, this is not only cumbersome to implement, but also increases the numerical effort.

by searching for the vector|ui of the Krylov subspace that minimises the distance to the exact time evolved state ˆU(δt)|ψi

U(δt)ˆ |ψ(t)i ≈ arg min

|ui∈K( ˆH,|ψ(t)i)

|||ui −Uˆ(δt)|ψ(t)i||² =:|ψ⁰(t+δt)i. (4.2.2) Since ˆU is linear, the solution of equation Eq. (4.2.2) is given by

|ψ⁰(t+δ)i= ˆP^†U(δ) ˆˆ P|ψ(t)i, (4.2.3) with ˆP :=^Pm−1j=0 |φjihφj| being the projector on the Krylov subspace Km( ˆH,|ψ(t)i). We now want to expand the time evolution operator in the previous equation

|ψ⁰(t+δt)i=^m−1X

j=0

|φjihφj|

∞

X

n=0

(−iδt)ⁿ

n! Hˆⁿ|φ0i

=^m−1X

j=0 m−1

X

n=0

(−iδt)ⁿ

n! |φjihφj|Hˆⁿ|φ0i+O δt^m m!

!

. (4.2.4)

Here we used the fact that the first basis vector|φ0iof the Krylov subspace is equal to|ψi and that for sufficiently largemwe can neglect the higher order terms of the exponential.

If the size of the Krylov subspace is chosen large enough, all vectors ˆHⁿ|φ0i ∈ K^m( ˆH,|ψi) are located in the subspace. Thus, we can replace ˆH with the effective matrix (Heff)j,k = hφj|Hˆ|φki and reintroduce the exponential under the assumption that m is sufficiently large

|ψ⁰(t+δt)i=^m−1X

j=0 m−1

X

n=0

(−iδt)ⁿ

n! |φji(H_effⁿ)j,0+O δt^m m!

!

=^m−1X

j=0

|φji(e^−iHeffδt)j,0 +O δt^m m!

!

. (4.2.5)

Heff can be computed very efficiently during the construction of the Krylov subspace basis and its dimension is in general very small, i.e. in most cases significantly smaller than 12.

This means the exponential can be computed fast through direct diagonalisation.

In the end, the time-evolved state is obtained by summing up the MPSes representing the Krylov vectors with the appropriate pre-factors given by (e^−iHeffδt)j,0. Since the summation of MPSes increases the bond dimension significantly, the time evolved state |ψ⁰(t+δt)i will have a much higher bond dimension than |ψ(t)i. To keep the bond dimension low for the next time step, a couple of SVDs are used to truncate |ψ⁰(t+δt)i while moving through the state from left to right and back. Afterwards, the next time step can be computed with the exact same approach.

Additionally, to the step size and truncation error in the Krylov method we also have to check convergence with respect to the Krylov space dimension. Since the subspace is build iteratively by adding newly calculated basis vectors, there is a good way of doing this.

Computing the effective matrix Heff and the coefficient vector (e^−iHeffδt)j,0 is cheap and can be done easily in every iteration of the creation of the Krylov basis. If the changes of this vector are below a desired threshold, which should be of the order of the truncation threshold, the size of the Krylov subspace is sufficiently large.

4.2.1 Variational Orthogonalisation

Lanczos algorithms and thus also the Krylov method are suffering from orthogonality problems. The iterative construction of the Krylov basis is inherently unstable for finite precision arithmetics. This becomes even worse for MPS arithmetics where the constant addition and subtraction of states leads to an increased bond dimension. It is there-fore crucial to truncate states to the desired precision to keep operations feasible. The truncation does not take the previously achieved orthogonality properties of the basis set into account and tends to destroy them even faster than typically for Lanczos methods.

Simply constantly checking for orthogonality and if necessary reorthogonalising the basis set, will not only increase the numerical effort significantly but also introduce new trun-cation errors, which may often entail the same problem as before. This happens because every orthogonalisation consists of additions and subtractions, which will increase the bond dimension of the MPS significantly. In our experience the better approach is to use a variational procedure to construct the basis vector under the constraint that it is orthogonal to all previous MPSes. This can be done with a previously set bond dimension and thus truncation becomes unnecessary.

To explain how this works, we define a set of M other MPSes {|o^Ai,|o^Bi, . . . ,|o^Mi}

against which we want to orthogonalise a given MPS |ψi on a system of size N to ob-tain a new state |ri. The corresponding tensors describing the MPS are {Ψl}^Nl=1 for |ψi, {{O_l^A}^Nl=1, . . . ,{O^M_l }^Nl=1}for theM MPSes we want to orthogonalise against, and{Rl}^Nl=1

for our resulting state |ri. The following procedure is done iteratively over all sites of the system. We therefore assume that each state is always in the mixed-canonical form with respect to site l. Our aim is to minimise the distance

|| |ψi − |ri||, (4.2.6)

under the constraint that the resulting state |ri is orthogonal to the other given MPSes hr|o^jJi = 0 ∀ J ∈ [A, M]. We can rewrite this problem by using Lagrange multipliers λ^J. Then, we have to minimise

hr|ψi+hr|ri+ ^XM

J=A

λ^Jhr|o^Ji= 0, (4.2.7)

with respect tohr| and λ^J. We now want to focus on optimising a single tensor on site l and differentiate with respect to the tensor R^†_l and λ^J to obtain the equations

−ψE +Rl+ ^XM

j=A

λ^Jo^J_E = 0, (4.2.8)

o^J_E ·R^†_l = 0⇔o^J†_E ·Rl = 0 ∀ J ∈[A, M]. (4.2.9) We define ΨE and o^J_E as the tensors obtained from taking the overlap between |ri and

|ψiand between |ri and |o^Ji respectively and removing the tensor R^†_l. Vectorising Rl as r, ψE asψ and collecting the vectorised tensors o^J_E as columns in the matrix o as well as the Lagrange multipliers in the vectorλ, we can summarise the equations in two compact matrix equation

r=ψ−oλ, (4.2.10)

0 =o^†r. (4.2.11)

Multiplying the first equation from the left witho^† allows us to insert the second equation into the first and solving for λ

λ= (o^†o)⁻¹o^†ψ. (4.2.12)

o^†ψ denotes the vector of overlaps between the vectors we want to orthogonalise against ho^J| and the input state |ψi, while o^†o is the matrix of overlaps between the states |o^Ji. Since we typically have a limited number of constraints, this matrix is in general relatively small and can be inverted exactly. Inserting Eq. (4.2.12) into Eq. (4.2.10) we obtain

r=ψ−o(o^†o)⁻¹o^†ψ, (4.2.13)

which can be reshaped back into the optimised tensor Rl. As in DMRG we move to the next site with help of an SVD and repeat the procedure there, thus optimising the state

|ri iteratively in the whole system. Since this is a variational approach, it is necessary to do several sweeps through the whole system to guarantee a converged and globally optimised result.

There are two issues with this approach that have to be treated carefully:

1. We optimise each tensor inside of the space spanned by a single tensor. Thus, we limit our bond dimension similarly to single site DMRG. This can result in a result-ing state |ri that differs from the result obtained from the normal Gram-Schmidt procedure. However, in general these differences are so small that we can neglect them and even if not, we can adjust the method to optimise two neighbouring sites at the same time similarly to two-site DMRG allowing for an increase of the bond dimension. Since the SVD used to split the merged sites is truncating the MPS again after we orthogonalised the tensor, it is necessary to add a few additional sweeps after the two-site orthogonalisation where we use single-site orthogonalisa-tion. The latter ensures that the orthogonalisation is exact since it does not require truncations.

2. The bond dimension at the boundaries of an MPS is typically very small and thus limits the variational space to e.g. only d parameters. Orthogonalising an MPS against multiple other states {|o^Ai,|o^Bi, . . . ,|o^Mi} in such a small optimisation space can be challenging or impossible. Our solution is to orthogonalise only against a small number of states in the beginning, if the norm of the resulting tensor Rl

does not become too small. While moving towards the centre of the MPS, the bond dimension and therefore the optimisation space become large enough that this criterium can be weakened or abandoned completely. In the following sweeps we give up this criterium completely to ensure that we obtain the solution of our optimisation.

4.2.2 Tensor-Optimised Implementation

There is a crucial difference between numerical linear algebra and tensor networks. For dense matrices the expectation value aXb, with a and b being vectors and X being a matrix only can be computed by first applying X to a and then computing the overlap withb. For tensor networks there is a difference between calculating the expectation value ha|Xˆ|bi represented by a tensor network directly or first evaluating ˆX|bi= |˜bi and then computingha|˜bi. In principle, the same contractions are performed but in different orders.

This makes a huge difference and in fact calculating the overlap is much cheaper than even the application ˆX|bi. To illustrate this, we displayed the contraction order of both cases in Fig. 4.2. Keeping this in mind, it is efficient to reorder the Krylov algorithm such that we need one MPO-MPS application less to calculated the Krylov space. Since we build up around several thousand Krylov subspaces per DMFT iteration, this sums up to a reasonable improvement of the time evolution. The adapted Krylov method is used as follows:

• The initial vector is set as first basis vector |φ0i = |ψi under the assumption that

|ψi is normalised.

• The first matrix element T0,0 =hφ0|Hˆ|φ0i is calculated.

• Any other basis vector |φk+1i for 0 ≤ k ≤ m−1 and the presentation of ˆH in the Krylov space ˆT is constructed from {|φ0i, . . . ,|φki} by

– computing|vki= ˆH|φki,

– orthogonalising the new vector with respect to all previous basis vectors

|wki=|vki −^Pkj=1hvk|φji|φji,

– normalising the new basis vector|φk+1i= _|||w^|wki

ki||, – computingTk+1,k =Tk,k+1 =hφk+1|Hˆ|φki,

– and computing Tk+1,k+1 =hφk+1|Hˆ|φk+1i.

W5

W2

W1 W3 W4 W6 W7 W8 W9 W10

A1 A2 A3 A4 M5 B6 B7 B8 B9 B10

A^∗₁ A^∗₂ A^∗₃ A^∗₄ M₅^∗ B₆^∗ B₇^∗ B₈^∗ B₉^∗ B₁₀^∗

W5

W2

W1 W3 W4 W6 W7 W8 W9 W10

A1 A2 A3 A4 M5 B6 B7 B8 B9 B10

A^∗₁ A^∗₂ A^∗₃ A^∗₄ M₅^∗ B₆^∗ B₇^∗ B₈^∗ B₉^∗ B₁₀^∗ a)

b)

Figure 4.2: Graphical representation of the expectation valuehψ|Hˆ|ψi written in MPS notation.

a) The expectation value can be computed by first applying the Hamiltonian on the state |φi = Hˆ|ψi and then determining the overlap hψ|φi. This means the tensors in the yellow block have to be contracted first and then with theA, M and B tensors on the top. b) The most efficient way of contracting the whole network would be to start with the left-most or right-most tensors (yellow block) and contracting them. The results is contracted with the tensors next to them up to the other end of the system.

The application of an MPO to an MPS becomes significantly costly for large MPOs.

Additionally, it is necessary to orthogonalise the resulting state against at least two other Krylov vectors. Therefore, another important tensor network related improvement is the application of an MPO to an MPS with a simultaneous orthogonalisation to a set of other vectors. This is done in the same manner as the orthogonalisation in section 4.2.1, except that the target environment tensor is not just the input MPS but the contraction of the input MPS|ψi, input MPO ˆH and the MPS |ri that is to be optimised.

This reduces the increase of computation time for each additional Krylov vectors (in line with the exponential growth of the bond dimension during very operator application step) to a much slower rate. It is recommended to introduce a subspace expansion step similarly to the DMRG procedure in section 3.2, to capture the increase of the minimal necessary bond dimension rather than generating an MPS with a large bond dimension and then truncating it.

4.2.3 Reusing the Krylov-Subspace

If we have a set of Krylov vectors that proved sufficient to calculate e^iHˆ(t+δt), that same Krylov subspace can also be sufficiently big to calculatee^iH(t+2δt)ˆ and potentiallye^iH(t+3δt)ˆ , e^iH(t+4δt)ˆ etc. Using the same Krylov subspace to calculate several time steps instead of computing a new Krylov subspace for each time step is desirable since it saves a lot of computational effort.

As described above, we need a certain number of Krylov vectors m to represent the Hamiltonian Heff in the Krylov subspace precisely enough for a good description of the time evolved state|ψ⁰(t+δt)i. From Eq. (4.2.5) it is obvious that this numbermdepends on the time step size δt since

|ψ⁰(t+δt)i=^m−1X

j=0

|φji(e^{−iHeff,exactδt})j,0+O δt^m m!

!

=^m−1X

j=0

|φji(e^−iHeffδte^{−iHerrorsδt})j,0+O δt^m m!

!

. (4.2.14)

The bigger the time step, the bigger the errors introduced by Heff, exact =Heff+Herrors in the exponential, since they get multiplied byδt. Additionally, the time step errorO^δtm!^m

itself is growing with larger time steps. Assume that the Krylov subspace used to deter-mine the time evolved state |ψ⁰(t+nδt)i cannot be used to compute an additional time step |ψ⁰(t+ (n+ 1)δt)i, because the errors would be too large. Adding another Krylov vectorm→m+ 1 can lower the errors enough for determining additional time steps with the same Krylov basis. The advantage is that the numerical effort for computing another Krylov vector can in some cases be smaller than setting up a complete new Krylov basis.

A. Swoboda in collaboration with us developed a heuristic estimator to decide whether building up a new Krylov basis or adding just a new Krylov vector to the already existing basis is better in terms of computation times. Apart from the actual computation cost for a single Krylov vector compared to a complete new Krylov basis, the heuristic also takes into account how many time steps can possibly be calculated with both ansatzes.

The heuristic is based on observations on the last calculations of Krylov vectors and is not simple to tune due to the increasing cost to generate each Krylov vector. Without truncating the generated Krylov vectors, the bond dimension of each new vector will grow by a factor proportional to the MPO bond dimension. That means, even if generating a third basis vector was much cheaper than calculating a new set of Krylov vectors in the previous step, generating a new basis consisting of three vectors can be much cheaper than adding a fourth one in the next iteration. Where exactly the crossover of both approaches is located depends on the initial state, the system and the time step. Thus, the heuristic has to be constantly adapted during the calculation of the time evolutions.

The reuse of the Krylov subspace is obviously more effective for imaginary time than for real time evolutions. The former is basically a projection to the ground state for τ → ∞. For short times the time evolved state will change fast, because of the exponentially fast decaying high energy contributions of the Hamiltonian spectrum, while for longer times

0 10 20 30 40 50 60 70 80 90 τ

0 5 10 15 20 25 30 35

runtime/min

classic

tensor optimised reused

Figure 4.3: Runtimes of Krylov time evolutions with different options for the two-site DCA model in momentum space with interaction strength U = 7t, nearest neighbour hopping t = 1, next-nearest neighbour hopping tp =−0.15t and five bath sites. The classic implementation of the Krylov method (blue line) is the slowest calculation. By adapting the method to be tensor optimised and orthogonalising the Krylov basis states variational (green line) the method can be improved significantly. Reusing the Krylov subspace based on our heuristic (red line) leads to the fastest calculation. Up to τ = 22 the heuristic is not optimally tuned and leads to a bad performance. From there on the Krylov subspace can be used for the rest of the time evolution leading to a very fast computation.

the time evolved state|ψ(t)iapproaches the ground state exponentially slow. Since|ψ(t)i is changing slowly, the Krylov basis will also change only slowly over time. Especially at later times in the time evolution, this allows to reuse the Krylov subspace over long time periods, which generates a massive decrease of computation times.

In contrast, for real-time evolutions the heuristic has problems to predict the number of needed Krylov vectors for the next time step, since the oscillating behaviour of the evo-lution and finite-size effects are hard to predict without taking into account the complete previous time evolution.

4.2.4 Conclusion

With all improvements mentioned, time evolution with the Krylov method becomes rel-atively fast and well-controlled. Furthermore, the Krylov method can be easily applied to all kind of Hamiltonians without any changes of the method and generalises very eas-ily to other tensor network states such as binary tree tensors. The only parts needing adaptation for other network topologies are the calculation of expectation values as well as the operator application with simultaneous orthogonalisation against previous Krylov vectors.

In Fig. 4.3 the runtime changes of the different improvements can be seen for a simple imaginary-time evolution. In this case all improvements lead to a reduction of the runtime by a factor 3.5. Depending on the model and the form of the time evolution the runtime reduction can be even bigger.

Im Dokument Dynamical mean-field theory studies on real materials (Seite 49-57)