Exponential of a skew Hermitian matrix

[jebej]

Would there be interest in having a function for taking the exponential of a skew-Hermitian (anti-Hermitian) matrix in the standard library, ie: expm(i*A) where A is Hermitian or Symmetric?

function expim(A::Union{Symmetric,Hermitian})
    # First decompose A into U*Λ*Uᴴ
    (Λ,U) = eig(A)
    # U must be a complex matrix
    U = complex.(U,0.0)
    # Calculate the imaginary exponential of each eigenvalue and multiply
    # by U to obtain B = U*exp(iΛ)
    B = scale!(copy(U),complex.(cos.(Λ),sin.(Λ)))
    # Finally multiply B by Uᴴ to obtain U*exp(iΛ)*Uᴴ = exp(i*A)
    return B*U'
end

It can be much faster than expm(complex.(0.0,A)):

julia> A = Hermitian(rand(1000,1000));

julia> B = Hermitian(rand(100,100));

julia> C = Hermitian(rand(10,10));

julia> @benchmark expm(complex.(0.0,A))
BenchmarkTools.Trial:
  memory estimate:  686.67 MiB
  allocs estimate:  123
  --------------
  minimum time:     2.558 s (4.62% GC)
  median time:      2.654 s (6.86% GC)
  mean time:        2.654 s (6.86% GC)
  maximum time:     2.750 s (8.94% GC)
  --------------
  samples:          2
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

julia> @benchmark expim(A)
BenchmarkTools.Trial:
  memory estimate:  69.04 MiB
  allocs estimate:  28
  --------------
  minimum time:     489.935 ms (1.02% GC)
  median time:      505.759 ms (0.99% GC)
  mean time:        512.661 ms (3.93% GC)
  maximum time:     566.183 ms (12.66% GC)
  --------------
  samples:          10
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

julia> @benchmark expm(complex.(0.0,B))
BenchmarkTools.Trial:
  memory estimate:  6.42 MiB
  allocs estimate:  117
  --------------
  minimum time:     3.838 ms (0.00% GC)
  median time:      4.272 ms (0.00% GC)
  mean time:        5.037 ms (11.12% GC)
  maximum time:     13.668 ms (0.00% GC)
  --------------
  samples:          977
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

julia> @benchmark expim(B)
BenchmarkTools.Trial:
  memory estimate:  743.11 KiB
  allocs estimate:  27
  --------------
  minimum time:     2.787 ms (0.00% GC)
  median time:      2.912 ms (0.00% GC)
  mean time:        3.015 ms (1.83% GC)
  maximum time:     7.143 ms (0.00% GC)
  --------------
  samples:          1641
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

julia> @benchmark expm(complex.(0.0,C))
BenchmarkTools.Trial:
  memory estimate:  70.88 KiB
  allocs estimate:  72
  --------------
  minimum time:     66.639 μs (0.00% GC)
  median time:      67.809 μs (0.00% GC)
  mean time:        72.095 μs (3.05% GC)
  maximum time:     845.272 μs (85.55% GC)
  --------------
  samples:          10000
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

julia> @benchmark expim(C)
BenchmarkTools.Trial:
  memory estimate:  12.73 KiB
  allocs estimate:  20
  --------------
  minimum time:     59.332 μs (0.00% GC)
  median time:      64.301 μs (0.00% GC)
  mean time:        65.727 μs (0.89% GC)
  maximum time:     1.089 ms (90.74% GC)
  --------------
  samples:          10000
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

If not feel free to close this issue.

[tkelman]

might be worth making a package for skew symmetric matrices and operations on them, if there are many cases of specialized algorithms

[andreasnoack]

Just to follow up on @tkelman's comment. The usual approach in Julia would be to define a type for the matrix and then overload expm for that type.

[jebej]

Yes, I defined it in this way because there isn't a skew-hermitian type.

[stevengj]

Note that we already have an expi function, called cis, for scalars.

Now that exp(::Array) is deprecated, in a future release of Julia (1.0?) it might be reasonable to rename expm to exp, at which point cis(A) could be defined as the exp(iA) function.

[jebej]

Note that we already have an expi function, called cis, for scalars.

Argh, I knew there was such a function but I couldn't remember its name, and searching for "imaginary exponential" did not give anything. Where does the name "cis" come from?

Now that exp(::Array) is deprecated, in a future release of Julia (1.0?) it might be reasonable to rename expm to exp, at which point cis(A) could be defined as the exp(iA) function.

That would be ideal, although I would be partial to the name "expi".

[stevengj]

cis = cosine + i sine, and seems to be at least somewhat standard. (I like expi better, too, though I don't have a strong opinion.)

[yuyichao]

Doesn't it need to be cism? Or maybe we don't need to distinguish expm/exp and cism/cis anymore since the broadcast version is deprecated?

[simonbyrne]

Your code assumes A has real coefficients. A more general approach would be:

U * Diagonal(cis.(Λ)) * U'

Another case that I would be interested in is when the skew symmetric matrix itself is real-valued. This would probably require a new type, as suggested above (though I don't know if there are any shortcuts you could use for actual computation, other than converting to complex and back again)

[stevengj]

There are a few specialized algorithms for skew-symmetric and skew-Hermitian matrices, e.g.:

• Eigensystem Computation for Skew-Symmetric and a Class of Symmetric Matrices — like anything in ACM TOMS, the code is unusable in free software, and you shouldn't even look at the code.

• Exponentials of skew-symmetric matrices and logarithms of orthogonal matrices — a Padé method for e^skew.

[stevengj]

However, these skew-specific algorithms seem be specialized enough (there are hardly any papers on this topic) that a SkewSymmetric type and operations on it may be better left to a package.

[jebej]

@simonbyrne

Your code assumes A has real coefficients.

A can be complex, the eigenvalues however will be real.

Another case that I would be interested in is when the skew symmetric matrix itself is real-valued.

Note that the function above does not take directly skew-symmetric or skew-Hermitian matrices as input. If you have one of those, they will need to be multiplied by -i to be made Hermitian, at which point the above can be used.

@stevengj

It might be worth it to have cis (which really is what the code above is, after some improvements of course) in base, and leave more specialized algorithms and types for a package.

[simonbyrne]

What I meant was that:

julia> X = complex(randn(5,5),randn(5,5))
5×5 Array{Complex{Float64},2}:
   1.17223+0.682131im     -2.535-0.0439985im  …  -0.611447+1.37534im    -0.338699-0.234632im
  0.101541+0.452541im    1.03085-0.774498im        1.49028+0.0891255im  -0.945443-0.378772im
  -1.07659+0.582124im  -0.662286+0.463415im      -0.898571-0.794531im     1.13165+0.553503im
 -0.791463+0.354508im   0.663189-1.05251im        0.533912+0.687588im   -0.441968+0.692214im
  -0.51301-0.99587im    -1.20007-0.489075im      -0.286741+0.307531im    0.491975-0.189388im

julia> expim(Hermitian(X'*X))
ERROR: MethodError: no method matching complex(::Complex{Float64}, ::Float64)
Closest candidates are:
  complex(::Real, ::Real) at complex.jl:61
  complex(::Complex{T<:Real}) at complex.jl:63
  complex{T<:Real}(::AbstractArray{T<:Real,N}, ::Real) at complex.jl:842
 in broadcast_t(::Function, ::Type{Any}, ::Array{Complex{Float64},2}, ::Vararg{Any,N}) at ./broadcast.jl:222
 in expim(::Hermitian{Complex{Float64},Array{Complex{Float64},2}}) at ./REPL[1]:5

[jebej]

Oh right yes I changed that when I realized I could simply call 'complex(A)'. Here is the function I have right now. I'm not sure whether it's the best we can do. With complex Hermitian matrices of small sizes it is slower to call the eigendecomposition than to call the regular exp function. We might just need to find a threshold under which to call exp.

function expim(A::Union{Symmetric,Hermitian})
# First decompose A into U*Λ*Uᴴ
(Λ,U) = eig(A)
# U must be a complex matrix
U = complex(U)
# Calculate the imaginary exponential of each eigenvalue and multiply
# each column of U to obtain B = U*exp(iΛ)
n = length(Λ)
B = similar(U)
for j = 1:n
    @inbounds a = Complex(cos(Λ[j]),sin(Λ[j]))
    @simd for i = 1:n
        @inbounds B[i,j] = a * U[i,j]
    end
end
# Finally multiply B by Uᴴ to obtain U*exp(iΛ)*Uᴴ = exp(i*A)
return B*U'
end