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

[seadra]

Imaginary exponentiation of a Hermitian matrix often appears in quantum mechanics and quantum computing. It'd be nice to have this built-in, because the current implementation without the specialization is very slow (compared to Mathematica's MatrixExp).

[stevengj]

See my comment above — the way is now clear to define cis(A) = exp(im*A), so we don't need a special matrix type (except for the relatively uncommon case of real skew-Hermitian).

A PR should be pretty easy and short since it could just handle the A::Hermitian case and fall back to exp otherwise.

[seadra]

I'm not sure if you were responding to me, but that's what I said

Imaginary exponentiation of a Hermitian matrix

which is cis(A) where A is Hermitian.

Regarding being "uncommon", again, this is a very common and pretty much a fundamental operation in quantum mechanics and quantum computing, but I can't speak for fields other than physics.

[stevengj]

Regarding being "uncommon", again, this is a very common and pretty much a fundamental operation in quantum mechanics and quantum computing

It's common to compute cis(A) where A is purely imaginary Hermitian, so that cis(A) is real? Doesn't quantum computing mostly work with complex quantities? (One of the three Pauli matrices is indeed purely imaginary, but for exponentiation with Pauli matrices there are analytical formulas.)

(Maybe you misunderstood me: I wasn't saying that Hermitian cis was rare, I was simply talking about whether we need a specialized exp(A) for real skew-Hermitian A, corresponding to cis(A) for imaginary Hermitian A.)

[simonbyrne]

I have had reason to do it, e.g. it appears in the geodesics of the Steifel manifold (p 310 of http://math.mit.edu/~edelman/publications/geometry_of_algorithms.pdf)

[seadra]

@stevengj Oh, no, what I mean is, in quantum mechanics, we indeed have exp(i A) where A is Hermitian ---or equivalently cis(A). Which indeed results in complex numbers in general. exp(i A) belongs to the special unitary group, which is typically represented using NxN complex matrices.

Regarding cis(A), I couldn't be sure what you meant by "so we don't need a special matrix type", since there already is a Hermitian type. I guess I'm on the same page now :)

For two-level quantum systems, as you mention, Pauli matrices can be used to expand the exponent. But it is very rarely the case that quantum systems are two-levels. For N-level systems, there is no such analytical formula (even for three-level systems for which Gell-Mann matrices replace Pauli matrices).

In quantum computing, two-level systems appear as qubits. Entangling operations involve at least four-level systems (a pair of qubits), for which SU(4) operations can't fit in the SU(2)xSU(2) subgroup.

[stevengj]

@seadra, I'm familiar with quantum mechanics (I have a PhD in physics).

All I'm saying is that we can leave the specific case of exp(real anti-Hermitian) = cis(purely imaginary Hermitian), which would require a new matrix type, to external packages. A cis(A) function with a special case for A Hermitian would be enough for the LinearAlgebra stdlib.

[There are applications for such real anti-Hermitian exponentials, but they are not generally quantum mechanics. (e.g. real anti-Hermitian exponentials arise in time propagators for dissipation-free classical wave equations, though in practice classical wave simulations typically use other methods).]

[seadra]

OK, and sure, I personally don't imagine I'd ever need cis of a purely imaginary Hermitian.

I was only talking about cis of a Hermitian, which is currently very slow in Julia compared to Mathematica.

[stevengj]

If someone wants to put together a pull request to add:

cis(A::AbstractMatrix) = exp!([float(-imag(a), real(a)) for a in A]) # fallback
function cis(A::RealHermSymComplexHerm)
    # optimized Hermitian case
end

etcetera to the LinearAlgebra stdlib, that would be great.

[StefanKarpinski]

Good to have a clear action to be taken here 🙂

[jebej]

I was testing again, and it seems that currently (new computer, Julia v1.1.1) the results aren't as clear cut as before. The exp matrix exponential function is much faster. Test script: https://gist.github.com/jebej/edd0abaf0e8c43de3967bc671d15b4b3

[stevengj]

@jebej, I'm not sure how you got these results. I got a speedup of about 2x for complex A and 4x for real A on my machine (Julia 1.3). I used the following implementation:

using BenchmarkTools, LinearAlgebra

Base.cis(A::AbstractMatrix) = exp(im*A) # fallback
function Base.cis(A::LinearAlgebra.RealHermSymComplexHerm)
    F = eigen(A)
    return F.vectors * Diagonal(cis.(F.values)) * F.vectors'
end
Base.cis(A::Diagonal) = Diagonal(cis.(A.diag))

compared to the naive implementation cis_naive(A) = exp(im*A) and benchmarked for both real and complex Hermitian matrices using BenchmarkTools.

Getting a factor of 2–4 speedup on a common operation with ~5 lines of code (+ docs and tests) seems worth including in LinearAlgebra.

[olof3]

(I like expi better, too, though I don't have a strong opinion.)

I would prefer expim as well (for the scalar method too), I had never heard of the cis function before today. If the name cis is kept it would be nice if it at least is mentioned in the docs for exp.

Somewhat unrelated, but cis(z) and exp(im*z) behave differently for extended complex numbers such as z = Inf + im, with exp(im*z) = NaN + NaN*im while cis(z) throws ERROR: DomainError with Inf: sincos(x) is only defined for finite x.

[stevengj]

Seems like this is closed by #40194.

The remaining task, defining specialized eigensolvers and exp for real anti-symmetric matrices, seems best left to a package as discussed above.