I’m trying to use Julia’s ExponentialUtilities.jl package to perform matrix-exponential-matrix multiplication. I have a sparse NxN Hamiltonian H, a small time step dt, and an NxN matrix A. The goal is to compute exp(-idt*H)A. In order to do this using the ExponentialUtilities.jl package, I first pre-allocate a Krylov subspace Ks. Then, I loop over the columns aj of A. At each iteration, I use the arnoldi!() function, H, and aj to fill the pre-allocated Krylov subspace. Then, I call expv!() to write the output into a temporary workspace vector, and copy this workspace vector into the jth column of the output NxN matrix.
My issue is that the expv!() function seems to be allocating a lot of memory, about 100 KiB to perform each matrix-exponential-vector multiplication for a 997x997 matrix and a vector of length 997. Is this something that I should expect? I was under the impression that if I preallocated the Krylov subspace, I wouldn’t have to allocate a lot of extra memory. Thank you in advance for any guidance on this matter.
Hi Chris, thank you for your response. I ran the allocation profiler and it’s indicating that the allocations are due to expv!(). Each time it’s called, it allocates 3 arrays. Moreover, there are internal calls to exponential!() that also allocate memory as well. As far as I can tell, all of the memory allocation is due to something that’s going on with expv!().
I noticed in the ExponentialUtilities.jl documentation that there’s an expv!() function that takes a cache as an argument as well:
Could this be used to reduce the number of allocations? I would like to try it out, but I cannot find instructions in the documentation on how to use the cache. Thanks again!
Hi @random_guy, did you get this resolved?
I am similarly interested in expv!
I was able to mostly replicate your code, and I think I’ve figured out the expv! cache argument, this reduces the allocations by 1, and reduces memory a bit. I don’t know if any more reductions are possible, but would be interested to find out.
using ExponentialUtilities
N = 997
Ks = KrylovSubspace{ComplexF64}(N)
A = rand(ComplexF64, N, N)
my_cache = ExponentialUtilities.ExpvCache{ComplexF64}(N)
workspace = Vector{ComplexF64}(undef, N)
dt = 0.01
for j in 1:N
aj = @view A[:,j]
arnoldi!(Ks, A, aj, ishermitian=false)
@time expv!(workspace, -1im*dt, Ks, cache = my_cache)
copyto!(aj, workspace)
end
This changed the allocations of each call to 28, with memory 86.203 KiB
If you ran that as posted then the remaining allocations very likely come from using the top-level scope.
Always wrap your code in a function to allow Julia to compile it:
using ExponentialUtilities
function benchmark(N)
Ks = KrylovSubspace{ComplexF64}(N)
A = rand(ComplexF64, N, N)
my_cache = ExponentialUtilities.ExpvCache{ComplexF64}(N)
workspace = Vector{ComplexF64}(undef, N)
dt = 0.01
for j in 1:N
aj = @view A[:,j]
arnoldi!(Ks, A, aj, ishermitian=false)
@time expv!(workspace, -1im*dt, Ks, cache = my_cache)
copyto!(aj, workspace)
end
return A
end
benchmark(997)
Thank you for your reply! Wrapping the code in a function using the code you posted reduced the allocations from 28 to 25, and reduced memory allocated from 86.203 KiB to 86.141 KiB, so it seems there is another source of allocations. Do you have another suggestion?
I took a quick look through the performance tips, but didn’t see anything obviously applicable.
Thanks Isaac! Sorry for the late response, I’ve been traveling. I will run some experiments with the cache and let you know if I can get the allocations down even further.
I don’t know about you, but an issue I ran into when using expv!() is that the pre-allocated Krylov subspace somehow accumulates NaNs when I start with a vector and a matrix to be exponentiated, and then iterate the matrix-exponential-vector multiplication around five or more times, which is problematic for my use case. I performed the exact same iterated multiplication using expv() and did not run into any issues with NaN accumulation so I just ran with that temporarily, but at some point if I want to increase the performance characteristics of my simulator I’ll need to figure out why using expv!() with the same Krylov subspace several times causes NaNs to mysteriously show up.
