From 879dc883f83d63ee782884e5827ec9f6cc2302f4 Mon Sep 17 00:00:00 2001 From: araujoms Date: Sun, 18 May 2025 01:22:36 +0200 Subject: [PATCH 1/2] type stability of matrix functions --- src/dense.jl | 68 +++++++++++++++++++++++++++-------------------- src/symmetric.jl | 4 +-- src/triangular.jl | 2 +- test/dense.jl | 25 +++++++---------- test/symmetric.jl | 53 +++++++++++++++++++++++++++++++++++- 5 files changed, 104 insertions(+), 48 deletions(-) diff --git a/src/dense.jl b/src/dense.jl index c1b99237..da3d422d 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -592,7 +592,7 @@ function schurpow(A::AbstractMatrix, p) end # if A has nonpositive real eigenvalues, retmat is a nonprincipal matrix power. - if isreal(retmat) + if eltype(A) <: Real && isreal(retmat) return real(retmat) else return retmat @@ -602,12 +602,11 @@ function (^)(A::AbstractMatrix{T}, p::Real) where T checksquare(A) # Quicker return if A is diagonal if isdiag(A) - TT = promote_op(^, T, typeof(p)) - retmat = copymutable_oftype(A, TT) - for i in diagind(retmat, IndexStyle(retmat)) - retmat[i] = retmat[i] ^ p + if T <: Real && any(<(0), diagview(A)) + return applydiagonal(x -> complex(x)^p, A) + else + return applydiagonal(x -> x^p, A) end - return retmat end # For integer powers, use power_by_squaring @@ -615,7 +614,7 @@ function (^)(A::AbstractMatrix{T}, p::Real) where T # If possible, use diagonalization if ishermitian(A) - return (Hermitian(A)^p) + return parent(Hermitian(A)^p) end # Otherwise, use Schur decomposition @@ -745,7 +744,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat end return A elseif ishermitian(A) - return copytri!(parent(exp(Hermitian(A))), 'U', true) + return parent(exp(Hermitian(A))) end ilo, ihi, scale = LAPACK.gebal!('B', A) # modifies A nA = opnorm(A, 1) @@ -918,10 +917,10 @@ julia> log(A) function log(A::AbstractMatrix) # If possible, use diagonalization if isdiag(A) && eltype(A) <: Union{Real,Complex} - if eltype(A) <: Real && all(>=(0), diagview(A)) - return applydiagonal(log, A) + if eltype(A) <: Real && any(<(0), diagview(A)) + return applydiagonal(log ∘ complex, A) else - return applydiagonal(log∘complex, A) + return applydiagonal(log, A) end elseif ishermitian(A) logHermA = log(Hermitian(A)) @@ -1004,13 +1003,14 @@ sqrt(::AbstractMatrix) function sqrt(A::AbstractMatrix{T}) where {T<:Union{Real,Complex}} if checksquare(A) == 0 return copy(float(A)) - elseif isdiag(A) && (T <: Complex || all(x -> x ≥ zero(x), diagview(A))) - # Real Diagonal sqrt requires each diagonal element to be positive - return applydiagonal(sqrt, A) + elseif isdiag(A) + if T <: Real && any(<(0), diagview(A)) + return applydiagonal(sqrt ∘ complex, A) + else + return applydiagonal(sqrt, A) + end elseif ishermitian(A) - sqrtHermA = sqrt(Hermitian(A)) - PS = parent(sqrtHermA) - return ishermitian(sqrtHermA) ? copytri_maybe_inplace(PS, 'U', true) : PS + return parent(sqrt(Hermitian(A))) elseif istriu(A) return triu!(parent(sqrt(UpperTriangular(A)))) elseif isreal(A) @@ -1044,7 +1044,7 @@ sqrt(A::TransposeAbsMat) = transpose(sqrt(parent(A))) Computes the real-valued cube root of a real-valued matrix `A`. If `T = cbrt(A)`, then we have `T*T*T ≈ A`, see example given below. -If `A` is symmetric, i.e., of type `HermOrSym{<:Real}`, then ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to find the cube root. Otherwise, a specialized version of the p-th root algorithm [^S03] is utilized, which exploits the real-valued Schur decomposition ([`schur`](@ref)) to compute the cube root. @@ -1077,7 +1077,7 @@ function cbrt(A::AbstractMatrix{<:Real}) elseif isdiag(A) return applydiagonal(cbrt, A) elseif issymmetric(A) - return cbrt(Symmetric(A, :U)) + return copytri!(parent(cbrt(Symmetric(A))), 'U') else S = schur(A) return S.Z * _cbrt_quasi_triu!(S.T) * S.Z' @@ -1118,7 +1118,7 @@ end Compute the matrix cosine of a square matrix `A`. -If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to compute the cosine. Otherwise, the cosine is determined by calling [`exp`](@ref). # Examples @@ -1160,7 +1160,7 @@ end Compute the matrix sine of a square matrix `A`. -If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to compute the sine. Otherwise, the sine is determined by calling [`exp`](@ref). # Examples @@ -1265,7 +1265,7 @@ end Compute the matrix tangent of a square matrix `A`. -If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to compute the tangent. Otherwise, the tangent is determined by calling [`exp`](@ref). # Examples @@ -1357,7 +1357,7 @@ _subadd!!(X, Y) = X - Y, X + Y Compute the inverse matrix cosine of a square matrix `A`. -If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to compute the inverse cosine. Otherwise, the inverse cosine is determined by using [`log`](@ref) and [`sqrt`](@ref). For the theory and logarithmic formulas used to compute this function, see [^AH16_1]. @@ -1391,7 +1391,7 @@ end Compute the inverse matrix sine of a square matrix `A`. -If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to compute the inverse sine. Otherwise, the inverse sine is determined by using [`log`](@ref) and [`sqrt`](@ref). For the theory and logarithmic formulas used to compute this function, see [^AH16_2]. @@ -1425,7 +1425,7 @@ end Compute the inverse matrix tangent of a square matrix `A`. -If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to +If `A` is real-symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to compute the inverse tangent. Otherwise, the inverse tangent is determined by using [`log`](@ref). For the theory and logarithmic formulas used to compute this function, see [^AH16_3]. @@ -1436,8 +1436,8 @@ compute the inverse tangent. Otherwise, the inverse tangent is determined by usi ```julia-repl julia> atan(tan([0.5 0.1; -0.2 0.3])) 2×2 Matrix{ComplexF64}: - 0.5+1.38778e-17im 0.1-2.77556e-17im - -0.2+6.93889e-17im 0.3-4.16334e-17im + 0.5 0.1 + -0.2 0.3 ``` """ function atan(A::AbstractMatrix) @@ -1450,7 +1450,12 @@ function atan(A::AbstractMatrix) SchurF = Schur{Complex}(schur(A)) U = im * UpperTriangular(SchurF.T) R = triu!(parent(log((I + U) / (I - U)) / 2im)) - return SchurF.Z * R * SchurF.Z' + retmat = SchurF.Z * R * SchurF.Z' + if eltype(A) <: Real + return real(retmat) + else + return retmat + end end """ @@ -1493,7 +1498,12 @@ function asinh(A::AbstractMatrix) SchurF = Schur{Complex}(schur(A)) U = UpperTriangular(SchurF.T) R = triu!(parent(log(U + sqrt(I + U^2)))) - return SchurF.Z * R * SchurF.Z' + retmat = SchurF.Z * R * SchurF.Z' + if eltype(A) <: Real + return real(retmat) + else + return retmat + end end """ diff --git a/src/symmetric.jl b/src/symmetric.jl index 9187caf1..1db803c4 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -867,7 +867,7 @@ function ^(A::SymSymTri{<:Complex}, p::Real) return Symmetric(schurpow(A, p)) end -for func in (:cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh, :cbrt) +for func in (:cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :cbrt) @eval begin function ($func)(A::SelfAdjoint) F = eigen(A) @@ -890,7 +890,7 @@ function cis(A::SelfAdjoint) return nonhermitianwrappertype(A)(retmat) end -for func in (:acos, :asin) +for func in (:acos, :asin, :atanh) @eval begin function ($func)(A::SelfAdjoint) F = eigen(A) diff --git a/src/triangular.jl b/src/triangular.jl index 3365ba7c..961bfb83 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -1996,7 +1996,7 @@ function powm!(A0::UpperTriangular, p::Real) A[i, i] = -A[i, i] end # Compute the Padé approximant - c = 0.5 * (p - m) / (2 * m - 1) + c = (p - m) / (4 * m - 2) triu!(A) S = c * A Stmp = similar(S) diff --git a/test/dense.jl b/test/dense.jl index 89bf340b..c2847898 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -639,8 +639,8 @@ end sinA1 = convert(Matrix{elty}, [0.2865568596627417 -1.107751980582015 -0.13772915374386513; -0.6227405671629401 0.2176922827908092 -0.5538759902910078; -0.6227405671629398 -0.6916051440348725 0.3554214365346742]) - @test @inferred(cos(A1)) ≈ cosA1 - @test @inferred(sin(A1)) ≈ sinA1 + @test cos(A1) ≈ cosA1 + @test sin(A1) ≈ sinA1 cosA2 = convert(Matrix{elty}, [-0.6331745163802187 0.12878366262380136 -0.17304181968301532; 0.12878366262380136 -0.5596234510748788 0.5210483146041339; @@ -663,22 +663,22 @@ end # Identities for (i, A) in enumerate((A1, A2, A3, A4, A5)) - @test @inferred(sincos(A)) == (sin(A), cos(A)) + @test sincos(A) == (sin(A), cos(A)) @test cos(A)^2 + sin(A)^2 ≈ Matrix(I, size(A)) @test cos(A) ≈ cos(-A) @test sin(A) ≈ -sin(-A) - @test @inferred(tan(A)) ≈ sin(A) / cos(A) + @test tan(A) ≈ sin(A) / cos(A) @test cos(A) ≈ real(exp(im*A)) @test sin(A) ≈ imag(exp(im*A)) @test cos(A) ≈ real(cis(A)) @test sin(A) ≈ imag(cis(A)) - @test @inferred(cis(A)) ≈ cos(A) + im * sin(A) + @test cis(A) ≈ cos(A) + im * sin(A) - @test @inferred(cosh(A)) ≈ 0.5 * (exp(A) + exp(-A)) - @test @inferred(sinh(A)) ≈ 0.5 * (exp(A) - exp(-A)) - @test @inferred(cosh(A)) ≈ cosh(-A) - @test @inferred(sinh(A)) ≈ -sinh(-A) + @test cosh(A) ≈ 0.5 * (exp(A) + exp(-A)) + @test sinh(A) ≈ 0.5 * (exp(A) - exp(-A)) + @test cosh(A) ≈ cosh(-A) + @test sinh(A) ≈ -sinh(-A) # Some of the following identities fail for A3, A4 because the matrices are singular if i in (1, 2, 5) @@ -687,7 +687,7 @@ end @test @inferred(cot(A)) ≈ inv(tan(A)) @test @inferred(sech(A)) ≈ inv(cosh(A)) @test @inferred(csch(A)) ≈ inv(sinh(A)) - @test @inferred(coth(A)) ≈ inv(@inferred tanh(A)) + @test @inferred(coth(A)) ≈ inv(tanh(A)) end # The following identities fail for A1, A2 due to rounding errors; # probably needs better algorithm for the general case @@ -904,11 +904,6 @@ end end end -@testset "matrix logarithm is type-inferable" for elty in (Float32,Float64,ComplexF32,ComplexF64) - A1 = randn(elty, 4, 4) - @inferred Union{Matrix{elty},Matrix{complex(elty)}} log(A1) -end - @testset "Additional matrix square root tests" for elty in (Float64, ComplexF64) A11 = convert(Matrix{elty}, [3 2; -5 -3]) @test sqrt(A11)^2 ≈ A11 diff --git a/test/symmetric.jl b/test/symmetric.jl index 3f1397c0..934b69db 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -1199,13 +1199,64 @@ end end end -@testset "asin/acos/acosh for matrix outside the real domain" begin +@testset "asin/acos/acosh/tanh 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 + @test tanh(atanh(M2)) ≈ M2 + end +end + +@testset "type inference of matrix functions" begin + for T ∈ (Float32, Float64, ComplexF32, ComplexF64) + a = randn(T, 2, 2) + syma = Symmetric(a) + symtria = SymTridiagonal(syma) + herma = Hermitian(a) + #nasty functions + for f in (x->x^real(T)(0.3), sqrt, log, asin, acos, acosh, atanh) + if T <: Real + @test @inferred Matrix{Complex{T}} f(a) isa Matrix + @test @inferred Symmetric{Complex{T}} f(syma) isa Symmetric + @test @inferred Symmetric{Complex{T}} f(symtria) isa Symmetric + @test @inferred Symmetric{Complex{T}} f(herma) isa Union{Symmetric{Complex{T}}, Hermitian{T}} + else + @test @inferred f(a) isa Matrix{T} + @test @inferred Matrix{T} f(herma) isa Union{Matrix{T}, Hermitian{T}} + end + end + #nice functions + for f in (x->x^2, exp, cos, sin, tan, cosh, sinh, tanh, atan, asinh, cbrt) + if T <: Real + @test @inferred f(a) isa Matrix{T} + @test @inferred f(syma) isa Symmetric{T} + @test @inferred f(symtria) isa Symmetric{T} + @test @inferred f(herma) isa Hermitian{T} + else + f != cbrt && @test @inferred f(a) isa Matrix{T} + @test @inferred f(herma) isa Hermitian{T} + end + end + #special case cis + if T <: Real + @test @inferred cis(a) isa Matrix{Complex{T}} + @test @inferred cis(syma) isa Symmetric{Complex{T}} + @test @inferred cis(symtria) isa Symmetric{Complex{T}} + @test @inferred cis(herma) isa Symmetric{Complex{T}} + else + @test @inferred cis(a) isa Matrix{T} + @test @inferred cis(herma) isa Matrix{T} + end + #special case sincos + if T <: Real + @test @inferred sincos(syma) isa Tuple{Symmetric{T}, Symmetric{T}} + @test @inferred sincos(symtria) isa Tuple{Symmetric{T}, Symmetric{T}} + end + @test @inferred sincos(a) isa Tuple{Matrix{T}, Matrix{T}} + @test @inferred sincos(herma) isa Tuple{Hermitian{T}, Hermitian{T}} end end From bfeb066de82352d551c40058fda3d5ff1459a869 Mon Sep 17 00:00:00 2001 From: araujoms Date: Sun, 18 May 2025 09:31:22 +0200 Subject: [PATCH 2/2] missing test --- test/dense.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/test/dense.jl b/test/dense.jl index c2847898..30b8b93f 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1068,7 +1068,7 @@ end @test lyap(1.0+2.0im, 3.0+4.0im) == -1.5 - 2.0im end -@testset "$elty Matrix to real power" for elty in (Float64, ComplexF64) +@testset "$elty Matrix to real power" for elty in (Float32, Float64, ComplexF32, ComplexF64) # Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014 #Aa : only positive real eigenvalues Aa = convert(Matrix{elty}, [5 4 2 1; 0 1 -1 -1; -1 -1 3 0; 1 1 -1 2]) @@ -1094,7 +1094,13 @@ end ADi += [im 0; 0 im] end - for A in (Aa, Ab, Ac, Ad, Ah, ADi) + #ADin : negative Diagonal Matrix + ADin = convert(Matrix{elty}, [-3 0; 0 3]) + if elty <: LinearAlgebra.BlasComplex + ADin += [im 0; 0 im] + end + + for A in (Aa, Ab, Ac, Ad, Ah, ADi, ADin) @test A^(1/2) ≈ sqrt(A) @test A^(-1/2) ≈ inv(sqrt(A)) @test A^(3/4) ≈ sqrt(A) * sqrt(sqrt(A)) @@ -1107,7 +1113,7 @@ end end Tschurpow = Union{Matrix{real(elty)}, Matrix{complex(elty)}} - @test (@inferred Tschurpow LinearAlgebra.schurpow(Aa, 2.0)) ≈ Aa^2 + @test (@inferred Tschurpow LinearAlgebra.schurpow(Aa, real(elty)(2.0))) ≈ Aa^2 end @testset "BigFloat triangular real power" begin