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

cisv(t, A, b) to compute exp(im*A*t)*b for Hermitian A? #176

Open stevengj opened 4 months ago

stevengj commented 4 months ago

It is pretty common to want to compute unitary matrix exponentials of the form $e^{iAt}b$ for Hermitian matrices $A$. It would be nice to have a method that takes advantage of $A$ being Hermitian here, and right now this doesn't seem to exist:

The simplest approach to improving this would just be to allow t to be complex (or even an arbitrary Number type). Offhand, I don't see anything in the algorithm that requires it to be Real? That way, you could immediately take advantage of A being Hermitian in using a tridiagonal eigensolver as well as in the Krylov procedure (via Lanczos).

For dense matrix exponentials, Julia exports an optimized Hermitian method for cis(A) = exp(im*A). It would be natural to expose an analogous cisv(t, A, b) that computes exp(im*A*t)*b as well, and which under the hood calls exp(im*t, A, b) to exploit Hermitian A.

(An alternative is to have optimized expv methods for skew-Hermitian matrices, ala SkewLinearAlgebra.jl, but this seems like it would take more work. The main advantage of this would be to speed up the case where $iA$ is purely real, i.e. $e^{Bt} b$ where $B = -B^T$ is real skew-symmetric.)

ChrisRackauckas commented 3 months ago

I'm definitely not opposed to it. Someone just needs to implement it. It's only missing because we haven't ran into it yet, but it's a good self-contained project.

goerz commented 6 days ago

Allowing t to be a complex Number would be a prerequisite for using ExponentialUtilities as a backend for QuantumPropagators:

https://github.com/JuliaQuantumControl/QuantumPropagators.jl/issues/86

Absorbing 1im in A (the Hamiltonian) isn't really an option in that context

stevengj commented 6 days ago

Seems like it would be just a few-line change (+ tests) to change t::Real to t::Number.

goerz commented 6 days ago

Yeah, I'll probably try KrylovKit first, and then when I get around to ExponentialUtilities, I'll make a PR here (unless someone beats me to it)