Exponential

src/factorizations/matrixalgebrakit.jl

@@ -8,6 +8,7 @@ for f in
         :eig_full, :eig_vals, :eigh_full, :eigh_vals,
         :left_polar, :right_polar,
         :project_hermitian, :project_antihermitian, :project_isometric,
+        :exponential,
     ]
     f! = Symbol(f, :!)
     @eval function MAK.default_algorithm(::typeof($f!), ::Type{T}; kwargs...) where {T <: AbstractTensorMap}
@@ -49,6 +50,7 @@ for f! in (
         :qr_null!, :lq_null!,
         :svd_vals!, :eig_vals!, :eigh_vals!,
         :project_hermitian!, :project_antihermitian!, :project_isometric!,
+        :exponential!,
     )
     @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm)
         $(f! in (:eig_vals!, :eigh_vals!, :project_hermitian!, :project_antihermitian!) && :(LinearAlgebra.checksquare(t)))
@@ -62,6 +64,31 @@ for f! in (
     end
 end
 
+# Exponential with Tuple
+function MAK.exponential!((τ, t)::Tuple{E, T}, N, alg::AbstractAlgorithm) where {E <: Number, T <: AbstractTensorMap}
+    LinearAlgebra.checksquare(t)
+    foreachblock(t, N) do _, (tblock, Nblock)
+        Nblock′ = exponential!((τ, tblock), Nblock, alg)
+        # deal with the case where the output is not the same as the input
+        Nblock === Nblock′ || copy!(Nblock, Nblock′)
+        return nothing
+    end
+    return N
+end
+
+# Default algorithm for exponential with Tuple
+MAK.exponential!((τ, t)::Tuple{E, T}) where {E <: Number, T <: DiagonalTensorMap} = MAK.exponential!((τ, t), DefaultAlgorithm())
+MAK.exponential!(t::Tuple{E, T}, ::DefaultAlgorithm) where {E <: Number, T <: Diagonal} = MAK.exponential!(t, DiagonalAlgorithm())
+MAK.exponential!(t::Tuple{E, T1}, out::T2, ::DefaultAlgorithm) where {E <: Number, T1 <: Diagonal, T2 <: Diagonal} = MAK.exponential!(t, out, DiagonalAlgorithm())
+
+function MAK.default_algorithm(::typeof(exponential!), ::Type{Tuple{E, T}}; kwargs...) where {E <: Number, T}
+    return MAK.default_algorithm(exponential!, blocktype(T); kwargs...)
+end
+
+function MAK.copy_input(::typeof(exponential), (τ, t)::Tuple{E, T}) where {E <: Number, T <: AbstractTensorMap}
+    return (τ, copy_oftype(t, factorisation_scalartype(exponential, t)))
+end
+
 MAK.zero!(t::AbstractTensorMap) = zerovector!(t)
 
 # Default algorithm
@@ -76,6 +103,7 @@ for f in [
         :left_polar, :right_polar,
         :left_orth, :right_orth, :left_null, :right_null,
         :project_hermitian, :project_antihermitian, :project_isometric,
+        :exponential,
     ]
     f! = Symbol(f, :!)
     @eval MAK.$f!(t::AbstractTensorMap, alg::DefaultAlgorithm) =
@@ -221,3 +249,9 @@ MAK.initialize_output(::typeof(project_antihermitian!), tsrc::AbstractTensorMap,
     tsrc
 MAK.initialize_output(::typeof(project_isometric!), tsrc::AbstractTensorMap, ::AbstractAlgorithm) =
     similar(tsrc)
+
+# Exponential
+# ----------------
+MAK.initialize_output(::typeof(exponential!), t::AbstractTensorMap, ::AbstractAlgorithm) = t
+MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{Number, AbstractTensorMap}, ::AbstractAlgorithm) = t
+MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{T1, AbstractTensorMap{T2}}, ::AbstractAlgorithm) where {T1 <: Complex, T2 <: Real} = similar(t, complex(eltype(t)))

src/tensors/linalg.jl

@@ -416,15 +416,7 @@ function Base.:(/)(t1::AbstractTensorMap, t2::AbstractTensorMap)
     return t
 end
 
-# TensorMap exponentation:
-function exp!(t::TensorMap)
-    domain(t) == codomain(t) ||
-        error("Exponential of a tensor only exist when domain == codomain.")
-    for (c, b) in blocks(t)
-        copy!(b, LinearAlgebra.exp!(b))
-    end
-    return t
-end
+@deprecate exp!(t) exponential!(t)
 
 # Sylvester equation with TensorMap objects:
 function LinearAlgebra.sylvester(A::AbstractTensorMap, B::AbstractTensorMap, C::AbstractTensorMap)

test/tensors/exponential.jl

@@ -0,0 +1,64 @@
+using Test, TestExtras
+using TensorKit, MatrixAlgebraKit
+using Random
+
+spaces = [ℂ^4] #, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)]
+scalartypes = [Float64, ComplexF32] #, ComplexF64]
+algorithms = [MatrixFunctionViaLA(), MatrixFunctionViaEig(DefaultAlgorithm()), MatrixFunctionViaEigh(DefaultAlgorithm())]
+
+@timedtestset "exponential for Hermitian matrices with $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
+    A = randn(st1, space, space)
+    A = project_hermitian!(A)
+    τ = rand(st2)
+
+    expA = @constinferred exponential(A)
+    expτA = @constinferred exponential((τ, A))
+
+    for alg in algorithms
+        expA2 = @constinferred exponential(A, alg)
+        expτA2 = @constinferred exponential((τ, A), alg)
+
+        @test expA ≈ expA2
+        @test expτA ≈ expτA2
+    end
+end
+
+@timedtestset "exponential! for general matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
+    A = randn(st1, space, space)
+    τ = rand(st2)
+
+    @test exponential!(copy(A)) == exponential!((1.0, copy(A)))
+
+    A2 = exponential!((τ, A))
+    if st1 <: Real && st2 <: Complex
+        @test objectid(A2) != objectid(A)
+    else
+        @test objectid(A2) == objectid(A)
+    end
+
+    expτA = exponential!((τ, copy(A)))
+    expminτA = exponential!((-τ, copy(A)))
+    @test expτA * expminτA ≈ id(scalartype(expτA), space)
+    @test expτA ≈ inv(expminτA)
+end
+
+@timedtestset "exponential! for diagonal matrices for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes
+    A = DiagonalTensorMap(rand(st1, reduceddim(space)), space)
+    τ = rand(st2)
+
+    exponential!(copy(A))
+    @test exponential!(copy(A)) ≈ exponential!((1.0, copy(A))) #, DiagonalAlgorithm())
+
+    A2 = @constinferred exponential!((τ, A))
+    @test A2 isa DiagonalTensorMap
+    if st1 <: Real && st2 <: Complex
+        @test objectid(A2) != objectid(A)
+    else
+        @test objectid(A2) == objectid(A)
+    end
+
+    expτA = exponential!((τ, copy(A)))
+    expminτA = exponential!((-τ, copy(A)))
+    @test expτA * expminτA ≈ id(scalartype(expτA), space)
+    @test expτA ≈ inv(expminτA)
+end