fix: make exponential hermitian to the last bit for complex input

`MatrixFunctionViaEigh` built the real case as the symmetric product
`VexpD * transpose(VexpD)`, exact by construction, but the complex case as
a plain `VexpD * V'` with no closing projection, so the result was only
approximately hermitian. `squareroot`, `logarithm` and `power` are all
exact for both scalar types, leaving `exponential` the odd one out.

Since `eigh_full!` returns real eigenvalues, `exp(τD/2)` is self-adjoint
whenever `τ` is real, and `exp(τA) = (V exp(τD/2)) * (V exp(τD/2))'` is
then hermitian by construction. Route the complex case through
`_mul_herm!`, the same mechanism the other three use.

This needed a second change to take effect for `exponential(A)` itself:
the implicit scalar was `one(eltype(A))`, i.e. `1.0 + 0.0im` for a complex
matrix, which sent a plain `exp(A)` down the complex-`τ` branch even though
the scalar was real. Use `one(real(eltype(A)))` instead. The same edit for
`MatrixFunctionViaEig` is behaviour-neutral -- its branch condition already
excluded a complex `A` -- but avoids needless complex-scalar arithmetic.

A genuinely complex `τ` is still left unprojected, which is not an
oversight: `exp(τA)` of a hermitian `A` is hermitian only for real `τ`, and
projecting would silently return a different matrix.

Verified against `LinearAlgebra.exp(Hermitian(A))` for `Float32`,
`Float64`, `ComplexF32` and `ComplexF64` with both `QRIteration` and
`DivideAndConquer`, and for `Diagonal` with complex storage.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: lkdvos

src/implementations/exponential.jl

@@ -62,8 +62,10 @@ function exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaLA)
     return expA
 end
 
-exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaEigh) = exponential!((one(eltype(A)), A), expA, alg)
-exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaEig) = exponential!((one(eltype(A)), A), expA, alg)
+# the implicit scalar is real even for a complex matrix, which keeps `exp(A)` of a hermitian `A`
+# on the branch that produces a hermitian result
+exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaEigh) = exponential!((one(real(eltype(A))), A), expA, alg)
+exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaEig) = exponential!((one(real(eltype(A))), A), expA, alg)
 
 function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::AbstractAlgorithm)
     expA .= A .* τ
@@ -74,18 +76,21 @@ function exponential!((τ, A)::Tuple{Number, AbstractMatrix}, expA, alg::MatrixF
     check_input(exponential!, (τ, A), expA, alg)
     D, V = eigh_full!(A, alg.eigh_alg)
     if eltype(A) <: Real
+        # `exp(τA) = (V exp(τD/2)) * transpose(V exp(τD/2))` is symmetric by construction;
+        # `transpose` rather than `'` keeps this valid for a complex `τ` as well
         if eltype(τ) <: Real
             VexpD = rmul!(V, exponential!((τ / 2, D), D))
         else
             VexpD = V * exponential((τ / 2, D))
         end
         return mul!(expA, VexpD, transpose(VexpD))
+    elseif eltype(τ) <: Real
+        # `D` is real, so `exp(τD/2)` is self-adjoint for a real `τ` and
+        # `exp(τA) = (V exp(τD/2)) * (V exp(τD/2))'` is hermitian by construction
+        return _mul_herm!(expA, rmul!(V, exponential!((τ / 2, D), D)))
     else
-        if eltype(τ) <: Real
-            VexpD = V * exponential!((τ, D), D)
-        else
-            VexpD = V * exponential((τ, D))
-        end
+        # a complex `τ` makes `exp(τA)` normal but not hermitian, so there is nothing to project
+        VexpD = V * exponential((τ, D))
         return mul!(expA, VexpD, V')
     end
 end

test/testsuite/matrixfunctions/exponential.jl

@@ -80,16 +80,22 @@ function test_exponential_hermitian(
         @test logarithm(expA, alg) ≈ A
         @test A == Ac
 
-        # `exponential` forms the real case as a symmetric product `VexpD * transpose(VexpD)`,
-        # which is hermitian to the last bit, but the complex case as a plain `VexpD * V'` with
-        # no closing projection, so there it is only approximately hermitian. This differs from
-        # `squareroot`, `logarithm` and `power`, which are exact for both scalar types.
-        if exact_hermiticity && eltype(T) <: Real
+        # `exp(A)` is built as a product of a factor with its own adjoint (transpose, in the real
+        # case), so it comes out hermitian to the last bit for both scalar types.
+        if exact_hermiticity
             @test ishermitian(expA)
             test_spectrum && @test eigh_vals(expA) ≈ exp.(eigh_vals(A))
         else
             @test ishermitian(expA; rtol = precision(T))
         end
+
+        # the scaled entry point keeps that guarantee for a real `τ`, since `exp(τD/2)` is then
+        # still self-adjoint. A complex `τ` makes `exp(τA)` normal but not hermitian, which
+        # `test_exponential_scaled` covers, so only correctness is asserted there.
+        τ = randn(rng, real(eltype(T)))
+        expτA = @testinferred exponential((τ, A), alg)
+        @test expτA ≈ exponential(τ * A, alg)
+        exact_hermiticity && @test ishermitian(expτA)
     end
 end