# Fast Exponential of a large matrix

**URL:** <https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998>\
**Category:** Performance\
**Created:** [July 13, 2024, 1:27pm UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998 "2024-07-13T13:27:08Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 13, 2024, 1:27pm UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/1 "2024-07-13T13:27:08Z")

</div>

I have to perform the exponential of a matrix. I have compared the built-in `exp` function to the `exponential!` function from the [ExponentialUtilites](https://github.com/SciML/ExponentialUtilities.jl) package. The latter is supposed to be much faster than the built-in function. However, I don’t see this!

```julia
A = randn(2^12, 2^12)
@time exp(im*A);

 71.589389 seconds (2.02 M allocations: 1.880 GiB, 0.10% gc time, 2.65% compilation time)

A = randn(2^12, 2^12)
@time exponential!(im*A);

71.524002 seconds (3.93 M allocations: 1.741 GiB, 0.25% gc time, 5.79% compilation time)

```

Why is this so? And is there any other way to calculate the exponential of a large matrix efficiently, like using multi-threading?

---

<div class="post-metadata">

**Author:** ![tom-plaa](https://avatars.discourse-cdn.com/v4/letter/t/d26b3c/32.png) [@tom-plaa](https://discourse.julialang.org/u/tom-plaa)\
**Post date:** [July 13, 2024, 3:01pm UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/2 "2024-07-13T15:01:35Z")

</div>

If you’re willing to sacrifice accuracy maybe an idea could be to truncate its series expansion: [Matrix exponential - Wikipedia](https://en.wikipedia.org/wiki/Matrix_exponential)

---

<div class="post-metadata">

**Author:** ![freeman](https://avatars.discourse-cdn.com/v4/letter/f/ec9cab/32.png) [@freeman](https://discourse.julialang.org/u/freeman)\
**Post date:** [July 13, 2024, 3:13pm UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/3 "2024-07-13T15:13:00Z")

</div>

Wait, can’t this be done without allocating?

I’m not an expert, but… if we just think about a taylor expansion:

x + x^2 + x^3 + x^4 times some coefficients.

This sum can be split into an accumulator where the final result will be stored and some chunk of memory were the calculation is being done. The calculation is you just keep multiplying by the same matrix and can be done iteratively: each step is multiplying some NxN matrix corresponding to the output of the previous step, times x which is also an NxN matrix. This has known size and can be allocated beforehand.

Alternatively, isn’t this the kind of thing that BLAS should be really good at doing? Why not just use that?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 13, 2024, 3:41pm UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/4 "2024-07-13T15:41:54Z")

</div>

> [@devanshu](#):
>
> `exp(im*A)`

If you have a Hermitian `A` you should use `cis(A)` instead.

Do you need the whole matrix, or do you just need to multiply it by a vector? In the latter case there are faster possibilities.

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 14, 2024, 9:19am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/5 "2024-07-14T09:19:29Z")

</div>

Hi, @stevengj! Yes, the matrix is Hermitian. And yes, you are right; I just need to multiply `exp(im*A)` by a vector. Could you please tell me what other possibilities you are referring to? Thanks!

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 14, 2024, 9:23am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/6 "2024-07-14T09:23:36Z")

</div>

Hi @tom-plaa! I am not sure how much accuracy I will lose by doing a series truncation. But I would prefer to calculate `exp(im*A)` since it’s important to preserve the unitarity.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [July 14, 2024, 9:28am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/7 "2024-07-14T09:28:13Z")

</div>

`expv` which uses krylov methods is likely much better here.

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 14, 2024, 9:34am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/8 "2024-07-14T09:34:23Z")

</div>

@Oscar_Smith this is wonderful! Thanks!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 14, 2024, 11:18am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/9 "2024-07-14T11:18:10Z")

</div>

> [@freeman](#):
>
> I’m not an expert, but… if we just think about a taylor expansion:
> 
> x + x^2 + x^3 + x^4 times some coefficients.

It’s well known that Taylor series are a terrible way to compute matrix exponentials in general: [Taking gradients of a matrix exponential - #8 by stevengj](https://discourse.julialang.org/t/taking-gradients-of-a-matrix-exponential/107865/8)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 14, 2024, 11:42am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/10 "2024-07-14T11:42:49Z")

</div>

> [@devanshu](#):
>
> Yes, the matrix is Hermitian.

Unfortunately, it looks like ExponentialUtilities.jl can’t currently exploit this: [cisv(t, A, b) to compute exp(im\*A\*t)\*b for Hermitian A? · Issue #176 · SciML/ExponentialUtilities.jl · GitHub](https://github.com/SciML/ExponentialUtilities.jl/issues/176)

It might be worth trying KrylovKit.jl instead, since it [supports this case in its `exponentiate` function](https://jutho.github.io/KrylovKit.jl/stable/man/matfun/), by passing an imaginary t and a Hermitian A.

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 15, 2024, 4:44am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/11 "2024-07-15T04:44:13Z")

</div>

@stevengj it seems `expv` in `ExponentialUtilities` is much faster than `exponentiate` in `KrylovKit`

```julia
@time expv(im, B, phi);
0.491426 seconds (46 allocations: 1.065 MiB)

@time exponentiate(B, im, phi);
11.420562 seconds (973 allocations: 26.749 MiB, 0.10% gc time)

```

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [July 15, 2024, 6:19am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/12 "2024-07-15T06:19:13Z")

</div>

Check the accuracy of the result. In my experience ExponentialUtilities is faster than KrylovKit with default settings. However KrylovKit tries very hard to give you the most accurate result (by using restart IIRC) while ExponentialUtilities just gives you less accurate results.

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 15, 2024, 6:34am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/13 "2024-07-15T06:34:01Z")

</div>

@Oscar_Smith How to reduce the number of allocations for `expv`? So in my work, I have to repeatedly apply the `exp(im*H)`, where H is Hermitian, on a state, which gets updated in each run. Because of large memory allocations, I end up getting much slower computation with `expv`, than with `exponential!`. Here is a minimal working code for the latter

```julia
function run_sim(t_final, U, initial_state, final_state)
    for t in 1:t_final
        final_state .= U * initial_state
    end
    final_state
end

function main(n)

    d = 2^n
    initial_state = rand(ComplexF64, d)
    final_state = initial_state

    t_final = 10
    
    A = rand(ComplexF64, (d, d))
    H = (A .+ A') / 2
    
    U = exponential!(im*H)

    for _ in 1:10
        final_state .= run_sim(t_final, U, initial_state, final_state)
    end

end

@time main(10);
  1.326368 seconds (122 allocations: 145.607 MiB)

```

Here is the same code using `expv`

```julia
function run_sim(t_final, H, initial_state, final_state)
    for t in 1:t_final
        final_state .= expv(im, H, initial_state)
    end
    final_state
end

function main(n)

    d = 2^n
    initial_state = rand(ComplexF64, d)
    final_state = initial_state

    t_final = 10
    
    A = rand(ComplexF64, (d, d))
    H = (A .+ A') / 2

    for _ in 1:10
        final_state .= run_sim(t_final, H, initial_state, final_state)
    end

end

@time main(10);
  5.373086 seconds (4.91 k allocations: 101.364 MiB)

```

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [July 15, 2024, 7:08am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/14 "2024-07-15T07:08:54Z")

</div>

> [@devanshu](#):
>
> How to reduce the number of allocations for `expv`?

You can reuse some memory by using caches but unfortunately the [docs are broken](https://docs.sciml.ai/ExponentialUtilities/stable/expv/#Caches) for that. Here is how some of it works:

```julia
state = rand(ComplexF64, d)
krylovdim = 30 # default
# 1. create a Krylov subspace cache
KS = ExponentialUtilities.KrylovSubspace{ComplexF64}(length(state), min(krylovdim, length(state)))
# 2. populate Krylov space
ExponentialUtilities.arnoldi!(KS, H, state; ishermitian=true)
# 3. compute exponential, stores result back in state
ExponentialUtilities.expv!(state, -im*Δt, KS)
# now repeat 2. and 3.

```

I think this still allocates some memory but you cannot really work around it because you have complex valued states. But I am not quite sure about this.

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [July 15, 2024, 7:21am UTC](https://discourse.julialang.org/t/fast-exponential-of-a-large-matrix/116998/15 "2024-07-15T07:21:47Z")

</div>

It does reduce the allocations but is still large.

```julia
  5.003996 seconds (1.51 k allocations: 58.343 MiB)

```