exp for SkewHermitian slower than exp for skew-hermitian Matrix

[sethaxen]

Since this packaged specializes exp for SkewHermitian matrices, I expected that it would be faster than for the same skew-Hermitian matrix but stored as a Matrix, but I'm seeing that it is consistently slower:

julia> using Random, LinearAlgebra, SkewLinearAlgebra, BenchmarkTools

julia> skew(x) = (x - x') / 2;

julia> Random.seed!(42);

julia> A = skew(randn(100, 100));

julia> B = SkewHermitian(copy(A));

julia> exp(A) ≈ exp(B)
true

julia> @btime exp($A); @btime exp($B);
  915.888 μs (24 allocations: 471.25 KiB)
  1.215 ms (65 allocations: 847.39 KiB)

julia> Random.seed!(98);

julia> A = skew(randn(500, 500));

julia> B = SkewHermitian(copy(A));

julia> @btime exp($A); @btime exp($B);
  39.125 ms (26 allocations: 11.45 MiB)
  75.937 ms (72 allocations: 19.98 MiB)

Environment

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, tigerlake)
Threads: 8 default, 0 interactive, 4 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_THREADS = auto

julia> using Pkg; Pkg.status()
Status `/tmp/jl_AqfuIB/Project.toml`
  [5c889d49] SkewLinearAlgebra v1.0.0

[ararslan]

It looks like there's an efficient algorithm published recently(ish): https://arxiv.org/abs/2103.10132. They also have their code online here, but in Matlab with no license.

[stevengj]

I think what you're seeing here is mainly the consequence of multi-threading — in single-threaded operation, I'm seeing that exp is faster for SkewHermitian:

julia> using LinearAlgebra, SkewLinearAlgebra, BenchmarkTools

julia> B = skewhermitian!(randn(1000,1000)); A = Matrix(B);

julia> LinearAlgebra.BLAS.get_num_threads()
4

julia> @btime exp($A); @btime exp($B); # multi-threaded exp(A)
  150.355 ms (26 allocations: 45.95 MiB)
  240.712 ms (74 allocations: 80.32 MiB)

julia> LinearAlgebra.BLAS.set_num_threads(1)

julia> @btime exp($A); @btime exp($B); # single-threaded exp(A)
  514.780 ms (26 allocations: 45.95 MiB)
  375.475 ms (74 allocations: 80.32 MiB)

Making this package faster in the multi-threaded case requires work to further parallelize its eigensolver. However, if you calling exp from multiple Julia threads, you will want to set BLAS to single-threaded mode.