search

SciML / ExponentialUtilities.jl

Fast and differentiable implementations of matrix exponentials, Krylov exponential matrix-vector multiplications ("expmv"), KIOPS, ExpoKit functions, and more. All your exponential needs in SciML form.
https://docs.sciml.ai/ExponentialUtilities/stable/
Other
96 stars 31 forks source link

mode=:error_estimate gives wrong results for long times #160

Closed axsk closed 9 months ago

axsk commented 10 months ago

Describe the bug 🐞

expv(t, A, b, mode=:error_estimate) gives vastly different results for longer time-steps.

Expected behavior

I would expect error_estimate and happy_breakdown results to at least have the same magnitudes

Minimal Reproducible Example 👇

A = I - sprand(100,100,.1) * .1
b = rand(100)
@assert expv(10, A,b) ≈ expv(10, A,b ,mode=:error_estimate)

Error & Stacktrace ⚠️

julia>  expv(10, A,b)
100-element Vector{Float64}:
  19372.458459100264
 -57030.13477800453
   3852.973450336655
  -7684.450057227177

julia>  expv(10, A,b ,mode=:error_estimate)
100-element Vector{Float64}:
 -3.1939130045290945e15
  4.314219301303509e15
 -1.6738861310296682e15
 -4.276422508317541e15

Environment (please complete the following information):

julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 64 × AMD EPYC 7502P 32-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 11 on 64 virtual cores
Environment:
  JULIA_HISTORY = ./.history.jl_
daviehh commented 9 months ago

It's been a while since I last read through the code, but iirc the error estimation routine only implements the Lanczos algorithm, i.e. mostly for Hermitian (or, in your real-valued case, symmetric) matrices.

https://github.com/SciML/ExponentialUtilities.jl/blob/18e5161ad28e49e2b92c5f4879495dbd55a0b52e/src/krylov_phiv_error_estimate.jl#L3

https://github.com/SciML/ExponentialUtilities.jl/blob/18e5161ad28e49e2b92c5f4879495dbd55a0b52e/src/krylov_phiv_error_estimate.jl#L86-L101

Note with your example, with a symmetric A A += A' the results do agree w/

expv(10, A, b) - expv(10, A, b, mode=:error_estimate)

-- more details: the non error-estimation expv call uses the arnoldi method which only switches to the Lanczos step based on the isherimtian argument

https://github.com/SciML/ExponentialUtilities.jl/blob/18e5161ad28e49e2b92c5f4879495dbd55a0b52e/src/arnoldi.jl#L246-L247

however, the implementation of the error-estimation just checks the KrylovSubspace type

https://github.com/SciML/ExponentialUtilities.jl/blob/18e5161ad28e49e2b92c5f4879495dbd55a0b52e/src/krylov_phiv_error_estimate.jl#L37-L39

ChrisRackauckas commented 9 months ago

So if you do :error_estimate and are non Hermition, we should throw an error saying that it's an inappropriate choice of error estimator?

daviehh commented 9 months ago

@ChrisRackauckas Yes, would make sense since it assumes Lanczos/tridagonal KrylovSubspace.H to use with the stegr call

https://github.com/SciML/ExponentialUtilities.jl/blob/18e5161ad28e49e2b92c5f4879495dbd55a0b52e/src/krylov_phiv_error_estimate.jl#L29

correction in PR. Thanks!