Convergence issue of exponentiate function

[Qiaoyi-Li]

Dear Prof. Jutho,

Since the major time cost in MPS-TDVP comes from the multiplication of linear operator on vector, the computational time should be proportional to the krylovdim . As a result, a better convergence v.s. krylovdim is vital. However, I found KrylovKit.exponentiate seems to show bad convergence. The following code gives a simple comparison to the naive Lanczos implementation.

using LinearAlgebra
using KrylovKit
using CairoMakie

N = 100
KMax = 32
tol = 1e-8
A = rand(N, N)
A += A'
t = -1

v = rand(N)
normalize!(v)

v_ex = exp(t * A) * v

err_KrylovKit = map(1:KMax) do K
     v_k, info = exponentiate(A, t, v, Lanczos(
          krylovdim=K,
          maxiter=1,
          tol=tol,
          orth=ModifiedGramSchmidt(),
          eager=true,
          verbosity=0))
     norm(v_k - v_ex) / norm(v_k)
end


T = zeros(KMax + 1, KMax + 1)
w = Vector{typeof(v)}(undef, KMax + 1)
err_naive = zeros(KMax)
for k = 1:KMax
     if k == 1
          w[k] = v / norm(v)
     end
     w[k+1] = A * w[k]
     T[k, k] = dot(w[k+1], w[k])

     w[k+1] -= T[k, k] * w[k]
     if k > 1
          w[k+1] -= T[k-1, k] * w[k-1]
     end
     T[k, k+1] = T[k+1, k] = norm(w[k+1])

     normalize!(w[k+1])

     expT = exp(t * T[1:k, 1:k])
     v_k = mapreduce(i -> expT[i, 1] * w[i], +, 1:k)

     err_naive[k] = norm(v_k - v_ex) / norm(v_k)
end

The result is plotted below, it clearly shows that naive Lanczos implementation is better.

I'm not sure whether exponentiate is used in the expected way here. Or this function should be further optimized for the special usage such as v -> exp(tA)v?

[Jutho]

Thanks for reporting this. What is going on is that KrylovKit's exponentiate is not behaving as you might expect when the values of krylovdim or maxiter are not high enough to reach the required tolerance.

What KrylovKit does is to guarantee that, within the constraints given by krylovdim and maxiter, it does compute an approxmation to exp(t̃ * A) v₀ for some time t̃ with abs(t̃) ≤ abs(t) such that the requested tolerance (per unit time) is achieved, i.e. such that

norm(exp(t̃ * A) * v₀ - krylovkitresult) < tol * abs(t̃)

So what you see in your plot is that it does not compute the complete exponentiation for t=-1 up to the value K=21, where suddenly the error drops to a value below the requested tolerance (10^-8). Also, because you set eager=true, it also stops when it can compute the full t=-1 exponentiation with the requested tolerance, and so the results for K>21 are just identical (explaining the flat line at that point); the larger Krylov subspaces are never built. If you set eager=false (default value), the curve would go down exactly as with your simple implementation.

Now for the first part of the curve, and this specific termination behaviour of KrylovKit, I agree that

1. It is not extremely well documented (it is in the doc string, but should probably be stressed a lot more strongly)

2. It is probably not very intuitive and useful.

I do not remember what drove me to design exponentiate/expintegrator in this way, but I am certainly happy to change this behaviour to the more traditional behaviour. However, as this is technically a breaking change, it will require a version update to 0.7.

I am going to leave this issue open, so that I certainly do not forget about it.