Exponential

ext/MatrixAlgebraKitGenericSchurExt.jl

@@ -22,4 +22,9 @@ function MatrixAlgebraKit.eig_vals!(A::AbstractMatrix, D, ::GS_QRIteration)
     return GenericSchur.eigvals!(A)
 end
 
+function MatrixAlgebraKit.default_exponential_algorithm(E::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}}
+    eig_alg = MatrixAlgebraKit.default_eig_algorithm(E; kwargs...)
+    return MatrixFunctionViaEig(eig_alg)
+end
+
 end

src/MatrixAlgebraKit.jl

@@ -30,13 +30,15 @@ export left_polar, right_polar
 export left_polar!, right_polar!
 export left_orth, right_orth, left_null, right_null
 export left_orth!, right_orth!, left_null!, right_null!
+export exponential, exponential!, exponentiali, exponentiali!
 
 export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
     LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations,
     LAPACK_DivideAndConquer, LAPACK_Jacobi
 export GLA_HouseholderQR, GLA_QRIteration, GS_QRIteration
 export LQViaTransposedQR
 export PolarViaSVD, PolarNewton
+export MatrixFunctionViaLA, MatrixFunctionViaEig, MatrixFunctionViaEigh
 export DiagonalAlgorithm
 export NativeBlocked
 export CUSOLVER_Simple, CUSOLVER_HouseholderQR, CUSOLVER_QRIteration, CUSOLVER_SVDPolar,
@@ -81,9 +83,12 @@ include("common/matrixproperties.jl")
 
 include("yalapack.jl")
 include("algorithms.jl")
+
 include("interface/projections.jl")
 include("interface/decompositions.jl")
 include("interface/truncation.jl")
+include("interface/matrixfunctions.jl")
+
 include("interface/qr.jl")
 include("interface/lq.jl")
 include("interface/svd.jl")
@@ -93,6 +98,7 @@ include("interface/gen_eig.jl")
 include("interface/schur.jl")
 include("interface/polar.jl")
 include("interface/orthnull.jl")
+include("interface/exponential.jl")
 
 include("common/gauge.jl") # needs to be defined after the functions are
 
@@ -107,6 +113,7 @@ include("implementations/gen_eig.jl")
 include("implementations/schur.jl")
 include("implementations/polar.jl")
 include("implementations/orthnull.jl")
+include("implementations/exponential.jl")
 
 include("pullbacks/qr.jl")
 include("pullbacks/lq.jl")

src/common/view.jl

@@ -20,6 +20,26 @@ See also [`diagview`](@ref).
 
 diagonal(v::AbstractVector) = Diagonal(v)
 
+"""
+    map_diagonal!(f, dst, src...)
+
+Map the scalar function `f` over all elements of the diagonal of `src...`, returning
+a diagonal result.
+
+See also [`map_diagonal!`](@ref).
+"""
+map_diagonal(f, src, srcs...) = diagonal(f.(diagview(src), map(diagview, srcs)...))
+
+"""
+    map_diagonal!(f, dst, src...)
+
+Map the scalar function `f` over all elements of the diagonal of `src...`,
+into the diagonal elements of destination `dst`.
+
+See also [`map_diagonal`](@ref).
+"""
+map_diagonal!(f, dst, src, srcs...) = (diagview(dst) .= f.(diagview(src), map(diagview, srcs)...); dst)
+
 # triangularind
 function lowertriangularind(A::AbstractMatrix)
     Base.require_one_based_indexing(A)

src/implementations/exponential.jl

@@ -0,0 +1,136 @@
+# Inputs
+# ------
+function copy_input(::typeof(exponential), A::AbstractMatrix)
+    return copy!(similar(A, float(eltype(A))), A)
+end
+
+copy_input(::typeof(exponential), A::Diagonal) = copy(A)
+
+function copy_input(::typeof(exponentiali), τ::Number, A::AbstractMatrix)
+    return τ, copy!(similar(A, complex(eltype(A))), A)
+end
+
+copy_input(::typeof(exponentiali), τ::Number, A::Diagonal) = τ, copy(A)
+
+function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, alg::AbstractAlgorithm)
+    m, n = size(A)
+    m == n || throw(DimensionMismatch("square input matrix expected. Got ($m,$n)"))
+    @check_size(expA, (m, m))
+    return @check_scalar(expA, A)
+end
+
+function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEigh)
+    if !ishermitian(A)
+        throw(DomainError(A, "Hermitian matrix was expected. Use `project_hermitian` to project onto the nearest hermitian matrix)"))
+    end
+    m, n = size(A)
+    m == n || throw(DimensionMismatch("square input matrix expected. Got ($m,$n)"))
+    @check_size(expA, (m, m))
+    return @check_scalar(expA, A)
+end
+
+function check_input(::typeof(exponential!), A::AbstractMatrix, expA::AbstractMatrix, ::DiagonalAlgorithm)
+    m, n = size(A)
+    @assert m == n && isdiag(A)
+    @assert expA isa Diagonal
+    @check_size(expA, (m, m))
+    @check_scalar(expA, A)
+    return nothing
+end
+
+function check_input(::typeof(exponentiali!), A::AbstractMatrix, expA::AbstractMatrix, alg::AbstractAlgorithm)
+    m, n = size(A)
+    m == n || throw(DimensionMismatch("square input matrix expected. Got ($m,$n)"))
+    return @check_size(expA, (m, m))
+end
+
+function check_input(::typeof(exponentiali!), A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEigh)
+    if !ishermitian(A)
+        throw(DomainError(A, "Hermitian matrix was expected. Use `project_hermitian` to project onto the nearest hermitian matrix)"))
+    end
+    m, n = size(A)
+    m == n || throw(DimensionMismatch("square input matrix expected. Got ($m,$n)"))
+    return @check_size(expA, (m, m))
+end
+
+function check_input(::typeof(exponentiali!), A::AbstractMatrix, expA::AbstractMatrix, ::DiagonalAlgorithm)
+    m, n = size(A)
+    @assert m == n && isdiag(A)
+    @assert expA isa Diagonal
+    return @check_size(expA, (m, m))
+end
+
+# Outputs
+# -------
+initialize_output(::typeof(exponential!), A::AbstractMatrix, ::AbstractAlgorithm) = A
+initialize_output(::typeof(exponentiali!), τ::Number, A::AbstractMatrix, ::AbstractAlgorithm) =
+    complex(A)
+
+# Implementation
+# --------------
+function exponential!(A, expA, alg::MatrixFunctionViaLA)
+    check_input(exponential!, A, expA, alg)
+    return LinearAlgebra.exp!(A)
+end
+
+function exponential!(A, expA, alg::MatrixFunctionViaEigh)
+    check_input(exponential!, A, expA, alg)
+    D, V = eigh_full!(A, alg.eigh_alg)
+    expD = map_diagonal!(x -> exp(x / 2), D, D)
+    VexpD = rmul!(V, expD)
+    return mul!(expA, VexpD, V')
+end
+
+function exponential!(A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEig)
+    check_input(exponential!, A, expA, alg)
+    D, V = eig_full!(A, alg.eig_alg)
+    expD = map_diagonal!(exp, D, D)
+    iV = inv(V)
+    VexpD = rmul!(V, expD)
+    if eltype(A) <: Real
+        expA .= real.(VexpD * iV)
+    else
+        mul!(expA, VexpD, iV)
+    end
+    return expA
+end
+
+function exponentiali!(τ::Number, A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaLA)
+    check_input(exponentiali!, A, expA, alg)
+    expA .= A .* (im * τ)
+    return LinearAlgebra.exp!(expA)
+end
+
+function exponentiali!(τ::Number, A::AbstractMatrix, expA::AbstractMatrix, alg::MatrixFunctionViaEigh)
+    check_input(exponentiali!, A, expA, alg)
+    D, V = eigh_full!(A, alg.eigh_alg)
+    expD = map_diagonal(x -> exp(x * im * τ), D)
+    if eltype(A) <: Real
+        VexpD = V * expD
+        return expA .= real.(VexpD * V')
+    else
+        VexpD = rmul!(V, expD)
+        return mul!(expA, VexpD, V')
+    end
+end
+
+function exponentiali!(τ::Number, A, expA, alg::MatrixFunctionViaEig)
+    check_input(exponentiali!, A, expA, alg)
+    D, V = eig_full!(A, alg.eig_alg)
+    expD = map_diagonal!(x -> exp(x * im * τ), D, D)
+    iV = inv(V)
+    VexpD = rmul!(V, expD)
+    return mul!(expA, VexpD, iV)
+end
+
+# Diagonal logic
+# --------------
+function exponential!(A, expA, alg::DiagonalAlgorithm)
+    check_input(exponential!, A, expA, alg)
+    return map_diagonal!(exp, expA, A)
+end
+
+function exponentiali!(τ::Number, A, expA, alg::DiagonalAlgorithm)
+    check_input(exponentiali!, A, expA, alg)
+    return map_diagonal!(x -> exp(x * im * τ), expA, A)
+end

src/interface/exponential.jl

@@ -0,0 +1,60 @@
+# Exponential functions
+# --------------
+
+"""
+    exponential(A; kwargs...) -> expA
+    exponential(A, alg::AbstractAlgorithm) -> expA
+    exponential!(A, [expA]; kwargs...) -> expA
+    exponential!(A, [expA], alg::AbstractAlgorithm) -> expA
+
+Compute the exponential of the square matrix `A`,
+
+!!! note
+    The bang method `exponential!` optionally accepts the output structure and
+    possibly destroys the input matrix `A`. Always use the return value of the function
+    as it may not always be possible to use the provided `expA` as output.
+
+See also [`exponentiali(!)`](@ref exponentiali).
+"""
+@functiondef exponential
+
+"""
+    exponentiali(τ, A; kwargs...) -> expiτA
+    exponentiali(τ, A, alg::AbstractAlgorithm) -> expiτA
+    exponentiali!(τ, A, [expiτA]; kwargs...) -> expiτA
+    exponentiali!(τ, A, [expiτA], alg::AbstractAlgorithm) -> expiτA
+
+Compute the exponential of `i*τ*A`, where `i` is the imaginary unit, `τ` is a scalar, and `A` is a square matrix.
+This allows the user to use the hermitian eigendecomposition when `A` is hermitian, even when `i*τ*A` is not.
+
+!!! note
+    The bang method `exponentiali!` optionally accepts the output structure and
+    possibly destroys the input matrix `A`.
+    Always use the return value of the function as it may not always be
+    possible to use the provided `expiτA` as output.
+
+See also [`exponential(!)`](@ref exponential).
+"""
+@functiondef n_args = 2 exponentiali
+
+# Algorithm selection
+# -------------------
+default_exponential_algorithm(A; kwargs...) = default_exponential_algorithm(typeof(A); kwargs...)
+function default_exponential_algorithm(T::Type; kwargs...)
+    return MatrixFunctionViaLA(; kwargs...)
+end
+function default_exponential_algorithm(::Type{T}; kwargs...) where {T <: Diagonal}
+    return DiagonalAlgorithm(; kwargs...)
+end
+
+for f in (:exponential!,)
+    @eval function default_algorithm(::typeof($f), ::Type{A}; kwargs...) where {A}
+        return default_exponential_algorithm(A; kwargs...)
+    end
+end
+
+for f in (:exponentiali!,)
+    @eval function default_algorithm(::typeof($f), ::Tuple{A, B}; kwargs...) where {A, B}
+        return default_exponential_algorithm(B; kwargs...)
+    end
+end

src/interface/matrixfunctions.jl

@@ -0,0 +1,39 @@
+# ================================
+# EXPONENTIAL ALGORITHMS
+# ================================
+"""
+    MatrixFunctionViaLA()
+
+Algorithm type to denote finding the exponential of `A` via the implementation of `LinearAlgebra`.
+"""
+@algdef MatrixFunctionViaLA
+
+"""
+    MatrixFunctionViaEigh()
+
+Algorithm type to denote finding the exponential `A` by computing the hermitian eigendecomposition of `A`.
+The `eigh_alg` specifies which hermitian eigendecomposition implementation to use.
+"""
+struct MatrixFunctionViaEigh{A <: AbstractAlgorithm} <: AbstractAlgorithm
+    eigh_alg::A
+end
+function Base.show(io::IO, alg::MatrixFunctionViaEigh)
+    print(io, "MatrixFunctionViaEigh(")
+    _show_alg(io, alg.eigh_alg)
+    return print(io, ")")
+end
+
+"""
+    MatrixFunctionViaEig()
+
+Algorithm type to denote finding the exponential `A` by computing the eigendecomposition of `A`.
+The `eig_alg` specifies which eigendecomposition implementation to use.
+"""
+struct MatrixFunctionViaEig{A <: AbstractAlgorithm} <: AbstractAlgorithm
+    eig_alg::A
+end
+function Base.show(io::IO, alg::MatrixFunctionViaEig)
+    print(io, "MatrixFunctionViaEig(")
+    _show_alg(io, alg.eig_alg)
+    return print(io, ")")
+end