# Performance of \`exp(A)\` for 9x9 anti-Hermitian matrix: Julia vs. PyTorch vs. MATLAB (CPU & GPU)

**URL:** <https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696>\
**Category:** Performance\
**Tags:** question, performance\
**Created:** [August 19, 2025, 3:57am UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696 "2025-08-19T03:57:33Z")\
**Posts on this page:** 10\
**Page:** 2

<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:** [August 19, 2025, 2:12pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/21 "2025-08-19T14:12:19Z")

</div>

> [@abraemer](#):
>
> You can also try [SkewLinearAlgebra.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/SkewLinearAlgebra).

Note that this package is mainly for `exp(A)` for _real_ anti-Hermitian matrices. For the complex anti-Hermitian case, it doesn’t have any advantage over `cis(Hermitian(-im*A))` IIRC.

---

<div class="post-metadata">

**Author:** ![goerz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerz/32/3269_2.png) [@goerz](https://discourse.julialang.org/u/goerz)\
**Post date:** [August 19, 2025, 4:23pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/22 "2025-08-19T16:23:15Z")

</div>

> [@draftman9](#):
>
> Given that my matrix `A` is always anti-Hermitian (meaning `exp(A)` is unitary), are there specialized algorithms I should be using instead of a general-purpose `expm`?

Have a look at [`QuantumPropagators.Cheby`](https://juliaquantumcontrol.github.io/QuantumPropagators.jl/stable/api/quantumpropagators/#QuantumPropagatorsChebyAPI), and probably use it in combination with `StaticArrays`. If your “batch” has a common spectral envelope (e.g., these are the the same time-dependent Hamiltonian evaluated at points in time), the propagation `|Ψ(t)⟩ = exp(-i H t) |Ψ(0)⟩` boils down to a handful of matrix-vector multiplications, which should be extremely fast.

In general, calculating the _application_ of `U = exp(-i H t)` to a state is more efficient than calculating the matrix `U` (matrix-vector operations vs matrix-matrix operations). For very small matrices, your mileage may vary, as well as if you need to apply the exact same `U` to a lot of different states. But even then, an [expansion into Chebychev polynomials](https://juliaquantumcontrol.github.io/QuantumPropagators.jl/stable/methods/#method_cheby) would be an efficient way to calculate `U`, and it would be fairly straightforward to adapt [the code](https://github.com/JuliaQuantumControl/QuantumPropagators.jl/blob/79667c3cc60e4462a339e56bcb6239bcfe22f0df/src/cheby.jl#L224) to that.

Some benchmarks:

```julia-auto
julia> using BenchmarkTools

julia> using QuantumPropagators.Cheby

julia> using QuantumControlTestUtils.RandomObjects: random_matrix, random_state_vector

julia> Ψ = random_state_vector(9; );

julia> H = random_matrix(9; hermitian=true); dt = 1.0; # spectral radius ≈ 1.0

julia> @btime exp(-1im * $H * $dt) * $Ψ;
  3.724 μs (28 allocations: 14.42 KiB)

julia> A = -1im * H * dt;

julia> @btime exp($A);
  3.412 μs (22 allocations: 11.47 KiB)

julia> @btime ChebyWrk($Ψ, 2.0, -1.0, $dt); # one-time calculations
  371.966 ns (15 allocations: 1.48 KiB)

julia> wrk = ChebyWrk(Ψ, 2.0, -1.0, dt);

julia> @btime cheby!($Ψ, $H, $dt, $wrk);
  938.360 ns (0 allocations: 0 bytes)

```

With `StaticArrays`:

```julia-auto
julia> using StaticArrays

julia> Ψ = SVector{9}(random_state_vector(9; ));

julia> H = SMatrix{9,9}(random_matrix(9; hermitian=true));

julia> @btime exp(-1im * $H * dt) * $Ψ;
  4.613 μs (5 allocations: 4.30 KiB)

julia> @btime ChebyWrk($Ψ, 2.0, -1.0, $dt);
  358.057 ns (9 allocations: 1.27 KiB)

julia> wrk = ChebyWrk(Ψ, 2.0, -1.0, dt);

julia> @btime cheby($Ψ, $H, $dt, $wrk);
  367.823 ns (0 allocations: 0 bytes)

```

But these depend on the spectral range of your `H`, so you should run your own.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [August 19, 2025, 5:01pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/23 "2025-08-19T17:01:50Z")

</div>

`exp` from StaticArrays.jl appears to perform extremely (and unnecessarily) poorly above tiny sizes. It could really use some love.

I got 2x performance on `SMatrix{9, 9, Float32, 81}` just by replacing the lines `(V - U) \ (V + U)` with `lu(V - U) \ (V + U)` (which is to say that `\` might also benefit from some changes or less inlining, since these should be essentially identical).

With some no-op changes to the code – mostly shuffling the control flow, moving the `@evalpoly`s to separate functions, and replacing them with `evalpoly` (`@evalpoly`’s inability to outline seems horrible with `SMatrix`) – I was able to further boost this to 3x while massively reducing the TTFX (30s down to 2.5s – bad but better). I didn’t assess which parts of my changes were actually responsible for improvement, though a proper evaluation should.

Outlining the Pade approximants seems particularly useful because it makes actual algorithmic changes much cleaner to implement. For example, neither `StaticArrays` or `LinearAlgebra` coarsen the approximation for `Float32` although they probably could a little. And there’s no reason for `exp` (`StaticArrays` or generic `LinearAlgebra`) to have a custom inline [Paterson-Stockmeyer polynomial evaluation](https://en.wikipedia.org/wiki/Polynomial_evaluation#Matrix_polynomials) when that could be built into `evalpoly` for all to use. _(EDIT: wrong: there is efficiency in recycling parts of the PS evaluation but I still think having all the Pade’s defined inline in the function is excessive)_

---

<div class="post-metadata">

**Author:** ![draftman9](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/draftman9/32/38128_2.png) [@draftman9](https://discourse.julialang.org/u/draftman9)\
**Post date:** [August 20, 2025, 2:28pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/24 "2025-08-20T14:28:26Z")

</div>

> [@stevengj](#):
>
> For example, `cis(H)` takes advantage of the Hermitian nature of H whereas `exp(im*H)` does not, which makes a big difference for large arrays but for small arrays the overhead dominates

Thank you for this incredibly insightful explanation. This is a huge help and clarifies many of the confusing results I’ve been seeing. I will also add this note into the main poster according to your comments.

> [@stevengj](#):
>
> I tried this with `cis(H)`, implementing a version `mycis!(...)` that lets you pre-allocate all of the workspaces, and I got an allocation-free method, but it is still only comparable on my machine to `exp(im*H)` (which uses a different algorithm, not based on eigenvectors, but which doesn’t exploit symmetry).

Yeah I know. In my experience, in many cases, “pre-allocating” only can obtain a few advantages.

> [@stevengj](#):
>
> Note that this package is mainly for `exp(A)` for _real_ anti-Hermitian matrices. For the complex anti-Hermitian case, it doesn’t have any advantage over `cis(Hermitian(-im*A))` IIRC.

I did some test and found SkewLinearAlgebra.jl can take advantages on large matrix. Here is the testing codes:

```julia
using BenchmarkTools
using SkewLinearAlgebra, LinearAlgebra

num_threads = 1
N = 9

B = rand(N, N)
A_real = B - B'
A_real_s = SkewHermitian(A_real)

A_imag = collect(im * Symmetric(rand(N, N)))
A_imag_s = SkewHermitian(A_imag)

println("\n--- Testing on $(num_threads) thread(s) ---")
LinearAlgebra.BLAS.set_num_threads(num_threads)
println("BLAS threads: ", LinearAlgebra.BLAS.get_num_threads())

@btime exp($A_real)
@btime exp($A_real_s)
@btime exp($A_imag)
@btime exp($A_imag_s)

print()

```

Here is the results for `N = 9` (9x9 matrix):

```bash
--- Testing on 1 thread(s) ---
BLAS threads: 1
  2.178 μs (9 allocations: 4.73 KiB)
  4.543 μs (28 allocations: 11.14 KiB)
  4.243 μs (9 allocations: 8.50 KiB)
  9.700 μs (18 allocations: 16.39 KiB)

```

Here is the results for `N = 900` (900x900 matrix):

```bash
--- Testing on 1 thread(s) ---
BLAS threads: 1
  159.045 ms (15 allocations: 37.09 MiB)
  113.632 ms (40 allocations: 64.79 MiB)
  673.745 ms (15 allocations: 74.17 MiB)
  453.316 ms (24 allocations: 74.91 MiB)

```

---

<div class="post-metadata">

**Author:** ![draftman9](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/draftman9/32/38128_2.png) [@draftman9](https://discourse.julialang.org/u/draftman9)\
**Post date:** [August 20, 2025, 2:44pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/25 "2025-08-20T14:44:16Z")

</div>

> [@stevengj](#):
>
> > [@abraemer](#):
> >
> > You could also try to use [DifferentialEquations.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/DifferentialEquations) to solve this.
> 
> This is definitely worth trying, especially using `SMatrix` to specialize for 9x9, assuming H is varying mostly smoothly with time — I’m guessing a high-order scheme will be much better than using a simple fixed-timestep exponential integrator.

Thank you both! In my real cases, I need solve this differential equation with 2D coordinate space:

\frac{d}{dt}|\psi(\mathbf{R}, t)\rangle = -i \hat H(\mathbf{R}, t)|\psi(\mathbf{R}, t)\rangle

where \mathbf{R} is position vector. I’ve never used DifferentialEquations.jl in the past, so I need some time to think about how to use it.

---

<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:** [August 20, 2025, 4:51pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/26 "2025-08-20T16:51:29Z")

</div>

> [@draftman9](#):
>
> I did some test and found [SkewLinearAlgebra.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/SkewLinearAlgebra) can take advantages on large matrix.

In the complex case, you should compare it to `cis(Hermitian(-im*A))`, which takes advantage of the symmetry too, and gets roughly the same single-threaded performance for large complex skew-Hermitian A:

```julia-auto
julia> using SkewLinearAlgebra, LinearAlgebra, BenchmarkTools

julia> n = 900; A_s = skewhermitian!(randn(ComplexF64, n, n)); A = Matrix(A_s);

julia> BLAS.set_num_threads(1)

julia> @btime exp($A); # does not exploit symmetry
  1.660 s (26 allocations: 74.27 MiB)

julia> @btime cis(Hermitian(-im * $A)); # exploits symmetry
  436.643 ms (38 allocations: 75.03 MiB)

julia> @btime exp($A_s); # exploits symmetry
  453.617 ms (47 allocations: 75.07 MiB)

```

As I said, the advantage of SkewLinearAlgebra is mainly for the _real_ case (where `cis(-im * A)` is much more expensive because it converts real to complex):

```julia-auto
julia> n = 900; A_s = skewhermitian!(randn(n, n)); A = Matrix(A_s); # real skew-Hermitian

julia> exp(A) ≈ real(cis(Hermitian(-im * A))) ≈ exp(A_s)
true

julia> @btime exp($A); # does not exploit symmetry
  376.459 ms (26 allocations: 37.14 MiB)

julia> @btime real(cis(Hermitian(-im * $A))); # exploits symmetry but complexifies
  442.979 ms (41 allocations: 81.22 MiB)

julia> @btime exp($A_s); # exploits symmetry, stays real
  260.892 ms (74 allocations: 64.97 MiB)

```

---

<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:** [August 20, 2025, 5:04pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/27 "2025-08-20T17:04:40Z")

</div>

> [@draftman9](#):
>
> I need solve this differential equation with 2D coordinate space:

Note that DifferentialEquations.jl will be especially advantageous over explicit matrix exponentials if your Hamiltonian operator \hat{H}(t) is large and sparse (e.g. it is -\nabla^2 + V where \nabla^2 is a sparse finite-difference operator and V is diagonal) or otherwise can be multiplied by vectors quickly (e.g. parts of it are diagonal in Fourier space, as in DFT calculations in a planewave basis). Even if you wanted an explicit matrix exponential in such cases, you would be better off computing exponential–vector products by iterative methods.

But you were testing on tiny 9x9 matrices, which is confusingly small if your real problem involves a PDE-like operator. Methods and performance characteristics for tiny dense matrices are _very_ different from those for large sparse matrices!

---

<div class="post-metadata">

**Author:** ![andreasvarga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasvarga/32/11634_2.png) [@andreasvarga](https://discourse.julialang.org/u/andreasvarga)\
**Post date:** [August 27, 2025, 1:16pm UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/28 "2025-08-27T13:16:59Z")

</div>

Is the matrix `H(t)` periodic?

---

<div class="post-metadata">

**Author:** ![draftman9](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/draftman9/32/38128_2.png) [@draftman9](https://discourse.julialang.org/u/draftman9)\
**Post date:** [August 28, 2025, 1:14am UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/29 "2025-08-28T01:14:20Z")

</div>

No. Actually I need solve the differential equation below numerically.

i\frac{\partial}{\partial t} \begin{pmatrix} \chi\_1(\mathbf{R},t) \\ \chi\_2(\mathbf{R},t) \\ \chi\_3(\mathbf{R},t) \\ \vdots \end{pmatrix} = \begin{pmatrix} T+V\_1(\mathbf{R}) & \vec{\mu}\_{12}(\mathbf{R}) \cdot \vec{E}(t) & \vec{\mu}\_{13}(\mathbf{R}) \cdot \vec{E}(t) & \cdots \\ \vec{\mu}\_{12}(\mathbf{R}) \cdot \vec{E}(t) & T+V\_2(\mathbf{R}) & \vec{\mu}\_{23}(\mathbf{R}) \cdot \vec{E}(t) & \cdots \\ \vec{\mu}\_{13}(\mathbf{R}) \cdot \vec{E}(t) & \vec{\mu}\_{23}(\mathbf{R}) \cdot \vec{E}(t) & T+V\_3(\mathbf{R}) & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} \chi\_1(\mathbf{R},t) \\ \chi\_2(\mathbf{R},t) \\ \chi\_3(\mathbf{R},t) \\ \vdots \end{pmatrix}

At every position lattice \mathbf{R} and every time step t, `expm` of 9x9 matrix should apply on the quantum state. Neither the coulping \vec{\mu}\_{13}(\mathbf{R}) nor the \cos^2-envelop pulse \vec{E}(t) is periodic.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [August 28, 2025, 2:25am UTC](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696/30 "2025-08-28T02:25:07Z")

</div>

I’d take a good, long look at @goerz’s QuantumPropagators.jl (I’ve never used it, but that’s just because my work hasn’t been very quantum in the last several years). I would also check out the specialized Lie group ODE solvers in DifferentialEquations.jl, which you can read about here: [https://docs.sciml.ai/DiffEqDocs/stable/solvers/nonautonomous\_linear\_ode/#State-Dependent-Solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/nonautonomous_linear_ode/#State-Dependent-Solvers).

`expm` is great for time-independent problems, but I don’t know if there’s a good reason to prefer it over these specialized integrators when your Hamiltonian is time-dependent.

[Previous page](https://discourse.julialang.org/t/performance-of-exp-a-for-9x9-anti-hermitian-matrix-julia-vs-pytorch-vs-matlab-cpu-gpu/131696.md?page=1)