Keyword `tol` not having desired effect in `exponentiate` (`expintegrator`)

[emstoudenmire]

Using the exponentiate interface to expintegrator, I have been finding that the tol keyword does not behave as expected. Please see the minimal sample code to reproduce below.

I was expecting that by setting tol the number of operations numops would adaptively depend on tol and the total time t. Instead I am finding that the only thing which seems to reduce numops is lowering the krylovdim to a value around 10 or so. The error I find after the code runs, by comparing the returned vector to an exact solution, is usually extremely small (< 1E-15) regardless of how high I make the tol parameter (say if I set it to 1E-3 or even 0.1).

Code to reproduce – the main thing I've been doing is adjusting t from 0.1 up to 1.0 while raising and lowering tol and I see the behavior described above (the amount of linear operators, numops, is something like 22). Lowering krylovdim to 10 does make numops become smaller while still giving an accurate result.

using Random
import KrylovKit: exponentiate
using Printf
using LinearAlgebra

let
  Random.seed!(1)

  # --------------------
  # Adjustable Parameters:

  # Problem definition
  L = 20        # size of matrix
  t = 0.1       # time to evolve to
  ElT = Float64 # element type

  # exponentiate keyword args
  tol = 1E-5
  krylovdim = 30
  maxiter = 100

  #---------------------

  # Make a random Hermitian matrix H
  M = randn(ElT,L,L)
  H = (M + M')/2

  # Random starting vector (t=0)
  x0 = randn(ElT,L)

  # Define linear map A
  count = 0
  function A(x)
    count += 1
    return H*x
  end

  # Compute result from exponentiate
  xt_exp,info = exponentiate(A,-im*t,x0; tol, krylovdim, maxiter, ishermitian=true, verbosity=100)

  # Compute exact result
  xt = exp(-im*H*t)*x0

  println("exponentiate info:")
  @printf("  converged = %s\n",info.converged)
  @printf("  residual  = %.4E\n",info.residual)
  @printf("  normres   = %.4E\n",info.normres)
  @printf("  numiter   = %d\n",info.numiter)
  @printf("  numops    = %d (<- # times linear map A was applied)\n",info.numops)

  @assert count == info.numops # check that exponentiate is reporting this correctly

  @printf("Actual error |xt_exp - xt| = %.4E\n",norm(xt_exp - xt))

  return
end

Typical output for the above settings:

[ Info: Lanczos iteration step 1: normres = 3.549103723766127
[ Info: Lanczos iteration step 2: normres = 2.3614405541788916
[ Info: Lanczos iteration step 3: normres = 3.588689632432383
[ Info: Lanczos iteration step 4: normres = 2.7584942917329247
[ Info: Lanczos iteration step 5: normres = 4.132375600781351
[ Info: Lanczos iteration step 6: normres = 2.622109652985752
[ Info: Lanczos iteration step 7: normres = 2.6603824297597125
[ Info: Lanczos iteration step 8: normres = 3.014748524743459
[ Info: Lanczos iteration step 9: normres = 3.0273258357955317
[ Info: Lanczos iteration step 10: normres = 1.982100738953418
[ Info: Lanczos iteration step 11: normres = 2.409275684096494
[ Info: Lanczos iteration step 12: normres = 1.6056315015783447
[ Info: Lanczos iteration step 13: normres = 1.119638147747482
[ Info: Lanczos iteration step 14: normres = 1.1900999270425054
[ Info: Lanczos iteration step 15: normres = 0.8606215889963659
[ Info: Lanczos iteration step 16: normres = 1.2985006799870673
[ Info: Lanczos iteration step 17: normres = 0.27409751590136305
[ Info: Lanczos iteration step 18: normres = 0.7444688486178209
[ Info: Lanczos iteration step 19: normres = 0.049873174694519816
[ Info: Lanczos iteration step 20: normres = 2.1595912431228573e-31
[ Info: expintegrate finished after 1 iterations: total error = 9.039173795164503e-67
exponentiate info:
  converged = 1
  residual  = 0.0000E+00
  normres   = 9.0392E-67
  numiter   = 1
  numops    = 22 (<- # times linear map A was applied)
Actual error |xt_exp - xt| = 8.0934E-16

[Jutho]

Hi Miles, thanks for the report. I will take a look. Certainly this is very different then how tol works for linear problems and eigenvalues, as for exponentiate there is no notion of a residual. Hence, a particular error estimate is used (the details of which I will need to look up again); that is also why you can obtain unreasonably small values of the (estimated) error, below machine precision.

However, seeing your example, you are observing something else. It is indeed true that I currently don't have a mechanism to adaptively change the number of Krylov steps. Hence, I always compute a full Krylov subspace of the requested krylovdim (unless it terminates prematurely, indicating an invariant subspace). The error model is then used to split the total time step in smaller steps which are accurate up to the requested tolerance, given the current krylov subspace. However, if, as in this case, the whole evolution is captured accurately (more accurately then requested) within that Krylov subspace, then the computation just finishes after this single outer iteration.

The same happens with eigsolve, where by default I compute a full Krylov subspace, before assessing the norms of the residuals. However, there I have added an eager mode to compute the residuals after every single expansion step of the Krylov subspace, which allows to terminate eagerly, at the expense of incurring a small overhead as the Schur decomposition of the Rayleigh quotient needs to be computed every time. Something similar can probably be done for exponentiate.

[Jutho]

Actually, it turns out the eager mode was already implemented. My memory is terrible for these things. So just add eager=true to your keywords and you will see the desired effect. The error is still smaller than the requested tolerance, but that's a consequence of the error model, and the fact that the error can be much bigger with one krylov vector less, and then improve dramatically with one additional Krylov vector.

Note though, that if you set krylovdim and maxiter too low, so that convergence is not obtained, the behaviour of exponentiate will not be to approximate the full time evolution with a bigger error, but rather to evolve up to a smaller time within the requested tolerance. This is explained in the documentation.

[Jutho]

The most recent commit adds documentation for the eager behaviour in exponentiate and expintegrator. I will close this issue; feel free to reopen if you think anything else needs to be done.

[Jutho]

I've reopened this for further discussion.

Regarding the tolerance, also note that (as mentioned in the docs) tol is the requested accuracy per unit time. This means that, in your example with t=0.1, the requested tol on the result is 10 times smaller than the value that you enter. I don't know why I made that choice, but changing it would in principle be a breaking change, which is furthermore hard to have a deprecation/transition rule for.

[emstoudenmire]

Thanks for that helpful clarification. After thinking about it further, I believe that definition of tol is actually the most sensible one. Our code inside C++ ITensor happened to make a different choice (tol being the absolute error) so we will translate between them in the way you said when making any comparisons.

I think this and your earlier comments does address (close) this issue.