Revert "Cleanup and generalize functions of Hermitian matrices (#1340)" (#1357)

This reverts commit 611395d.

Unfortunately, this PR introduced some method ambiguities that weren't
caught by the CI, as the ambiguity test was broken. After the fix in
#1349, this now causes test failures on master
(https://buildkite.com/julialang/linearalgebra-dot-jl/builds/740#0196c86e-2d8c-42c2-b9ae-a2143427f43e).

Please re-land the PR after fixing the ambiguities.

CC: @araujoms

Author: jishnub

src/symmetric.jl

@@ -224,15 +224,11 @@ const RealHermSymSymTri{T<:Real} = Union{RealHermSym{T}, SymTridiagonal{T}}
 const RealHermSymComplexHerm{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}
 const RealHermSymComplexSym{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Symmetric{Complex{T},S}}
 const RealHermSymSymTriComplexHerm{T<:Real} = Union{RealHermSymComplexSym{T}, SymTridiagonal{T}}
-const SelfAdjoint = Union{SymTridiagonal{<:Real}, Symmetric{<:Real}, Hermitian}
+const SelfAdjoint = Union{Symmetric{<:Real}, Hermitian{<:Number}}
 
 wrappertype(::Union{Symmetric, SymTridiagonal}) = Symmetric
 wrappertype(::Hermitian) = Hermitian
 
-nonhermitianwrappertype(::SymSymTri{<:Real}) = Symmetric
-nonhermitianwrappertype(::Hermitian{<:Real}) = Symmetric
-nonhermitianwrappertype(::Hermitian) = identity
-
 size(A::HermOrSym) = size(A.data)
 axes(A::HermOrSym) = axes(A.data)
 @inline function Base.isassigned(A::HermOrSym, i::Int, j::Int)
@@ -838,74 +834,119 @@ end
 ^(A::Symmetric{<:Complex}, p::Integer) = sympow(A, p)
 ^(A::SymTridiagonal{<:Real}, p::Integer) = sympow(A, p)
 ^(A::SymTridiagonal{<:Complex}, p::Integer) = sympow(A, p)
-^(A::Hermitian, p::Integer) = sympow(A, p)
-function sympow(A, p::Integer)
+function sympow(A::SymSymTri, p::Integer)
+    if p < 0
+        return Symmetric(Base.power_by_squaring(inv(A), -p))
+    else
+        return Symmetric(Base.power_by_squaring(A, p))
+    end
+end
+for hermtype in (:Symmetric, :SymTridiagonal)
+    @eval begin
+        function ^(A::$hermtype{<:Real}, p::Real)
+            isinteger(p) && return integerpow(A, p)
+            F = eigen(A)
+            if all(λ -> λ ≥ 0, F.values)
+                return Symmetric((F.vectors * Diagonal((F.values).^p)) * F.vectors')
+            else
+                return Symmetric((F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors')
+            end
+        end
+        function ^(A::$hermtype{<:Complex}, p::Real)
+            isinteger(p) && return integerpow(A, p)
+            return Symmetric(schurpow(A, p))
+        end
+    end
+end
+function ^(A::Hermitian, p::Integer)
     if p < 0
         retmat = Base.power_by_squaring(inv(A), -p)
     else
         retmat = Base.power_by_squaring(A, p)
     end
-    return wrappertype(A)(retmat)
+    return Hermitian(retmat)
 end
-function ^(A::SelfAdjoint, p::Real)
+function ^(A::Hermitian{T}, p::Real) where T
     isinteger(p) && return integerpow(A, p)
     F = eigen(A)
     if all(λ -> λ ≥ 0, F.values)
         retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors'
-        return wrappertype(A)(retmat)
+        return Hermitian(retmat)
     else
-        retmat = (F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors'
-        return nonhermitianwrappertype(A)(retmat)
+        retmat = (F.vectors * Diagonal((complex.(F.values).^p))) * F.vectors'
+        if T <: Real
+            return Symmetric(retmat)
+        else
+            return retmat
+        end
     end
 end
-function ^(A::SymSymTri{<:Complex}, p::Real)
-    isinteger(p) && return integerpow(A, p)
-    return Symmetric(schurpow(A, p))
-end
 
 for func in (:exp, :cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh, :cbrt)
     @eval begin
-        function ($func)(A::SelfAdjoint)
+        function ($func)(A::RealHermSymSymTri)
+            F = eigen(A)
+            return wrappertype(A)((F.vectors * Diagonal(($func).(F.values))) * F.vectors')
+        end
+        function ($func)(A::Hermitian{<:Complex})
             F = eigen(A)
             retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
-            return wrappertype(A)(retmat)
+            return Hermitian(retmat)
         end
     end
 end
 
-function cis(A::SelfAdjoint)
+function cis(A::RealHermSymSymTri)
     F = eigen(A)
-    retmat = F.vectors .* cis.(F.values') * F.vectors'
-    return nonhermitianwrappertype(A)(retmat)
+    return Symmetric(F.vectors .* cis.(F.values') * F.vectors')
 end
+function cis(A::Hermitian{<:Complex})
+    F = eigen(A)
+    return F.vectors .* cis.(F.values') * F.vectors'
+end
+
 
 for func in (:acos, :asin)
     @eval begin
-        function ($func)(A::SelfAdjoint)
+        function ($func)(A::RealHermSymSymTri)
+            F = eigen(A)
+            if all(λ -> -1 ≤ λ ≤ 1, F.values)
+                return wrappertype(A)((F.vectors * Diagonal(($func).(F.values))) * F.vectors')
+            else
+                return Symmetric((F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors')
+            end
+        end
+        function ($func)(A::Hermitian{<:Complex})
             F = eigen(A)
             if all(λ -> -1 ≤ λ ≤ 1, F.values)
                 retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
-                return wrappertype(A)(retmat)
+                return Hermitian(retmat)
             else
-                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
-                return nonhermitianwrappertype(A)(retmat)
+                return (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
             end
         end
     end
 end
 
-function acosh(A::SelfAdjoint)
+function acosh(A::RealHermSymSymTri)
+    F = eigen(A)
+    if all(λ -> λ ≥ 1, F.values)
+        return wrappertype(A)((F.vectors * Diagonal(acosh.(F.values))) * F.vectors')
+    else
+        return Symmetric((F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors')
+    end
+end
+function acosh(A::Hermitian{<:Complex})
     F = eigen(A)
     if all(λ -> λ ≥ 1, F.values)
         retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors'
-        return wrappertype(A)(retmat)
+        return Hermitian(retmat)
     else
-        retmat = (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors'
-        return nonhermitianwrappertype(A)(retmat)
+        return (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors'
     end
 end
 
-function sincos(A::SelfAdjoint)
+function sincos(A::RealHermSymSymTri)
     n = checksquare(A)
     F = eigen(A)
     T = float(eltype(F.values))
@@ -915,28 +956,49 @@ function sincos(A::SelfAdjoint)
     end
     return wrappertype(A)((F.vectors * S) * F.vectors'), wrappertype(A)((F.vectors * C) * F.vectors')
 end
-
-function log(A::SelfAdjoint)
+function sincos(A::Hermitian{<:Complex})
+    n = checksquare(A)
     F = eigen(A)
-    if all(λ -> λ ≥ 0, F.values)
-        retmat = (F.vectors * Diagonal(log.(F.values))) * F.vectors'
-        return wrappertype(A)(retmat)
-    else
-        retmat = (F.vectors * Diagonal(log.(complex.(F.values)))) * F.vectors'
-        return nonhermitianwrappertype(A)(retmat)
+    T = float(eltype(F.values))
+    S, C = Diagonal(similar(A, T, (n,))), Diagonal(similar(A, T, (n,)))
+    for i in eachindex(S.diag, C.diag, F.values)
+        S.diag[i], C.diag[i] = sincos(F.values[i])
+    end
+    retmatS, retmatC = (F.vectors * S) * F.vectors', (F.vectors * C) * F.vectors'
+    for i in diagind(retmatS, IndexStyle(retmatS))
+        retmatS[i] = real(retmatS[i])
+        retmatC[i] = real(retmatC[i])
     end
+    return Hermitian(retmatS), Hermitian(retmatC)
 end
 
-# sqrt has rtol kwarg to handle matrices that are semidefinite up to roundoff errors
-function sqrt(A::SelfAdjoint; rtol = eps(real(float(eltype(A)))) * size(A, 1))
-    F = eigen(A)
-    λ₀ = -maximum(abs, F.values) * rtol # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff
-    if all(λ -> λ ≥ λ₀, F.values)
-        retmat = (F.vectors * Diagonal(sqrt.(max.(0, F.values)))) * F.vectors'
-        return wrappertype(A)(retmat)
-    else
-        retmat = (F.vectors * Diagonal(sqrt.(complex.(F.values)))) * F.vectors'
-        return nonhermitianwrappertype(A)(retmat)
+
+for func in (:log, :sqrt)
+    # sqrt has rtol arg to handle matrices that are semidefinite up to roundoff errors
+    rtolarg = func === :sqrt ? Any[Expr(:kw, :(rtol::Real), :(eps(real(float(one(T))))*size(A,1)))] : Any[]
+    rtolval = func === :sqrt ? :(-maximum(abs, F.values) * rtol) : 0
+    @eval begin
+        function ($func)(A::RealHermSymSymTri{T}; $(rtolarg...)) where {T<:Real}
+            F = eigen(A)
+            λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff
+            if all(λ -> λ ≥ λ₀, F.values)
+                return wrappertype(A)((F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors')
+            else
+                return Symmetric((F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors')
+            end
+        end
+        function ($func)(A::Hermitian{T}; $(rtolarg...)) where {T<:Complex}
+            n = checksquare(A)
+            F = eigen(A)
+            λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff
+            if all(λ -> λ ≥ λ₀, F.values)
+                retmat = (F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors'
+                return Hermitian(retmat)
+            else
+                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
+                return retmat
+            end
+        end
     end
 end
 

test/symmetric.jl

@@ -1199,14 +1199,4 @@ end
     end
 end
 
-@testset "asin/acos/acosh for matrix outside the real domain" begin
-    M = [0 2;2 0] #eigenvalues are ±2
-    for T ∈ (Float32, Float64, ComplexF32, ComplexF64)
-        M2 = Hermitian(T.(M))
-        @test sin(asin(M2)) ≈ M2
-        @test cos(acos(M2)) ≈ M2
-        @test cosh(acosh(M2)) ≈ M2
-    end
-end
-
 end # module TestSymmetric