From 54d5d10c7b94b3e41807feb8d1f042a52ec262e6 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 26 Jun 2026 17:36:20 +0200 Subject: [PATCH 1/7] Exponential via MatrixAlgebraKit --- src/tensors/linalg.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index 40a0e0024..c45c9f4f6 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -421,7 +421,16 @@ 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)) + MatrixAlgebraKit.exponential!(b, b) + end + return t +end + +function exp!((τ, t)::Tuple{Number, TensorMap}) + domain(t) == codomain(t) || + error("Exponential of a tensor only exist when domain == codomain.") + for (c, b) in blocks(t) + MatrixAlgebraKit.exponential!((τ, b), b) end return t end From caab0e48fc066979c915035c3ac9991d4db4caf9 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 26 Jun 2026 17:58:15 +0200 Subject: [PATCH 2/7] include case where t is real and tau is complex --- src/tensors/linalg.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index c45c9f4f6..f07f1b174 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -429,6 +429,13 @@ end function exp!((τ, t)::Tuple{Number, TensorMap}) domain(t) == codomain(t) || error("Exponential of a tensor only exist when domain == codomain.") + if eltype(t) <: Real && eltype(τ) <: Complex + t_complex = complex(t) + for ((cr, br), (cc, bc)) in zip(blocks(t), blocks(t_complex)) + MatrixAlgebraKit.exponential!((τ, br), bc) + end + return t_complex + end for (c, b) in blocks(t) MatrixAlgebraKit.exponential!((τ, b), b) end From 863ca6d28b709f16caafa96ab866b652d53c8214 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 26 Jun 2026 18:43:27 +0200 Subject: [PATCH 3/7] Add test --- test/tensors/exponential.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 test/tensors/exponential.jl diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl new file mode 100644 index 000000000..cb5befb4e --- /dev/null +++ b/test/tensors/exponential.jl @@ -0,0 +1,25 @@ +using Test, TestExtras +using TensorKit +using Random + +spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1//2 => 1)] +scalartypes = [Float64, ComplexF64] + +@timedtestset "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space = spaces, st1 = scalartypes, st2 = scalartypes + A = randn(st1, space, space) + τ = rand(st2) + + @test exp!(copy(A)) == exp!((1.0, copy(A))) + + A2 = exp!((τ,A)) + if st1 <: Real && st2 <: Complex + @test objectid(A2) != objectid(A) + else + @test objectid(A2) == objectid(A) + end + + expτA = exp!((τ,copy(A))) + expminτA = exp!((-τ,copy(A))) + @test expτA * expminτA ≈ id(scalartype(expτA), space) + @test expτA ≈ inv(expminτA) +end From bf8c3af6a3c254d21dbd599456c3d730668b00fc Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 26 Jun 2026 18:44:06 +0200 Subject: [PATCH 4/7] Fix formatting --- test/tensors/exponential.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl index cb5befb4e..51fe1072b 100644 --- a/test/tensors/exponential.jl +++ b/test/tensors/exponential.jl @@ -2,24 +2,24 @@ using Test, TestExtras using TensorKit using Random -spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1//2 => 1)] +spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)] scalartypes = [Float64, ComplexF64] -@timedtestset "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space = spaces, st1 = scalartypes, st2 = scalartypes +@timedtestset "exp!(τ,A) 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 exp!(copy(A)) == exp!((1.0, copy(A))) - A2 = exp!((τ,A)) + A2 = exp!((τ, A)) if st1 <: Real && st2 <: Complex @test objectid(A2) != objectid(A) else @test objectid(A2) == objectid(A) end - expτA = exp!((τ,copy(A))) - expminτA = exp!((-τ,copy(A))) + expτA = exp!((τ, copy(A))) + expminτA = exp!((-τ, copy(A))) @test expτA * expminτA ≈ id(scalartype(expτA), space) @test expτA ≈ inv(expminτA) end From 9db5a963cb986e26a5642a7011ed934daf890b4b Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 3 Jul 2026 15:37:16 +0200 Subject: [PATCH 5/7] Update correspondence with MAK and extend tests --- src/factorizations/matrixalgebrakit.jl | 34 ++++++++++++++++ test/tensors/exponential.jl | 55 ++++++++++++++++++++++---- 2 files changed, 81 insertions(+), 8 deletions(-) diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 83d53fb0b..87d20c042 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/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(complex(t)) diff --git a/test/tensors/exponential.jl b/test/tensors/exponential.jl index 51fe1072b..fbb2d73e9 100644 --- a/test/tensors/exponential.jl +++ b/test/tensors/exponential.jl @@ -1,25 +1,64 @@ using Test, TestExtras -using TensorKit +using TensorKit, MatrixAlgebraKit using Random -spaces = [ℂ^4, Vect[U1Irrep](0 => 1, 1 => 2), Vect[SU2Irrep](0 => 1, 1 // 2 => 1)] -scalartypes = [Float64, ComplexF64] +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 "exp!(τ,A) for $space, scalartype(A) = $st1, scalartype(τ) = $st2" for space in spaces, st1 in scalartypes, st2 in scalartypes +@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) - @test exp!(copy(A)) == exp!((1.0, copy(A))) + expA = @constinferred exponential(A) + expτA = @constinferred exponential((τ, A)) - A2 = exp!((τ, 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 = exp!((τ, copy(A))) - expminτA = exp!((-τ, copy(A))) + 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 From 6e5d094618676743d4031b687ae5b45535a7c8a7 Mon Sep 17 00:00:00 2001 From: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com> Date: Fri, 3 Jul 2026 17:35:40 +0200 Subject: [PATCH 6/7] Update src/factorizations/matrixalgebrakit.jl Co-authored-by: Lukas Devos --- src/factorizations/matrixalgebrakit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 87d20c042..7de5cb740 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -254,4 +254,4 @@ MAK.initialize_output(::typeof(project_isometric!), tsrc::AbstractTensorMap, ::A # ---------------- 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(complex(t)) +MAK.initialize_output(::typeof(exponential!), (τ, t)::Tuple{T1, AbstractTensorMap{T2}}, ::AbstractAlgorithm) where {T1 <: Complex, T2 <: Real} = similar(t, complex(eltype(t))) From 121de729da6deba590024571bad09665c7c87057 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 3 Jul 2026 17:37:04 +0200 Subject: [PATCH 7/7] deprecate `exp!` --- src/tensors/linalg.jl | 26 +------------------------- 1 file changed, 1 insertion(+), 25 deletions(-) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index f07f1b174..0852d9c27 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -416,31 +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) - MatrixAlgebraKit.exponential!(b, b) - end - return t -end - -function exp!((τ, t)::Tuple{Number, TensorMap}) - domain(t) == codomain(t) || - error("Exponential of a tensor only exist when domain == codomain.") - if eltype(t) <: Real && eltype(τ) <: Complex - t_complex = complex(t) - for ((cr, br), (cc, bc)) in zip(blocks(t), blocks(t_complex)) - MatrixAlgebraKit.exponential!((τ, br), bc) - end - return t_complex - end - for (c, b) in blocks(t) - MatrixAlgebraKit.exponential!((τ, b), b) - end - return t -end +@deprecate exp!(t) exponential!(t) # Sylvester equation with TensorMap objects: function LinearAlgebra.sylvester(A::AbstractTensorMap, B::AbstractTensorMap, C::AbstractTensorMap)