From 17024d6cc8b08d6c8a7073acb5e1b9bfaa2e65be Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 16:23:13 -0400 Subject: [PATCH 01/22] feat: add `squareroot`, `logarithm` and `power` matrix functions Following the `exponential` template, each function supports `MatrixFunctionViaLA` (wrapping the LinearAlgebra implementation), `MatrixFunctionViaEig`/`MatrixFunctionViaEigh` (eigendecomposition-based, also serving generic eltypes via GenericSchur/GenericLinearAlgebra), and `DiagonalAlgorithm`. `power(A, p)` uses the two-argument `@functiondef` form with an exponent-type-based split between integer and fractional powers. The output scalar type always matches the input: out-of-domain eigenvalues (negative real axis for real input, zero eigenvalues for `logarithm` and negative fractional powers) throw a `DomainError`, while eigenvalues that violate the domain within a tolerance `domain_atol` (new keyword of `MatrixFunctionViaEig(h)`, defaulting to the new `default_domain_atol`) are clamped as rounding artifacts. --- ext/MatrixAlgebraKitGenericSchurExt.jl | 12 +++ src/MatrixAlgebraKit.jl | 10 +++ src/common/defaults.jl | 9 +++ src/implementations/logarithm.jl | 88 +++++++++++++++++++++ src/implementations/matrixfunctions.jl | 91 ++++++++++++++++++++++ src/implementations/power.jl | 102 +++++++++++++++++++++++++ src/implementations/squareroot.jl | 86 +++++++++++++++++++++ src/interface/logarithm.jl | 41 ++++++++++ src/interface/matrixfunctions.jl | 30 ++++++-- src/interface/power.jl | 42 ++++++++++ src/interface/squareroot.jl | 40 ++++++++++ 11 files changed, 545 insertions(+), 6 deletions(-) create mode 100644 src/implementations/logarithm.jl create mode 100644 src/implementations/matrixfunctions.jl create mode 100644 src/implementations/power.jl create mode 100644 src/implementations/squareroot.jl create mode 100644 src/interface/logarithm.jl create mode 100644 src/interface/power.jl create mode 100644 src/interface/squareroot.jl diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index c34f1abc4..93c8f5823 100644 --- a/ext/MatrixAlgebraKitGenericSchurExt.jl +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -21,6 +21,18 @@ function MatrixAlgebraKit.default_exponential_algorithm( return MatrixFunctionViaEig(eig_alg) end +for default_f_algorithm in ( + :default_squareroot_algorithm, :default_logarithm_algorithm, + :default_power_algorithm, + ) + @eval function MatrixAlgebraKit.$default_f_algorithm( + type::Type{T}; domain_atol::Union{Nothing, Real} = nothing, kwargs... + ) where {T <: StridedMatrix{<:GSFloat}} + eig_alg = MatrixAlgebraKit.default_eig_algorithm(type; kwargs...) + return MatrixFunctionViaEig(eig_alg; domain_atol) + end +end + function geev!(::GS, A::AbstractMatrix, Dd::AbstractVector, V::AbstractMatrix; kwargs...) D, Vmat = GenericSchur.eigen!(A) copyto!(Dd, D) diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 4d4e0084e..3b6165349 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -31,6 +31,9 @@ export left_polar!, right_polar! export left_orth, right_orth, left_null, right_null export left_orth!, right_orth!, left_null!, right_null! export exponential, exponential! +export squareroot, squareroot! +export logarithm, logarithm! +export power, power! export Householder, Native_HouseholderQR, Native_HouseholderLQ export DivideAndConquer, SafeDivideAndConquer, QRIteration, Bisection, Jacobi, SVDViaPolar @@ -115,6 +118,9 @@ include("interface/schur.jl") include("interface/polar.jl") include("interface/orthnull.jl") include("interface/exponential.jl") +include("interface/squareroot.jl") +include("interface/logarithm.jl") +include("interface/power.jl") include("implementations/projections.jl") include("implementations/truncation.jl") @@ -128,6 +134,10 @@ include("implementations/schur.jl") include("implementations/polar.jl") include("implementations/orthnull.jl") include("implementations/exponential.jl") +include("implementations/matrixfunctions.jl") +include("implementations/squareroot.jl") +include("implementations/logarithm.jl") +include("implementations/power.jl") include("common/gauge.jl") # needs to be defined after the functions are diff --git a/src/common/defaults.jl b/src/common/defaults.jl index c64c24ad4..008d2ca6b 100644 --- a/src/common/defaults.jl +++ b/src/common/defaults.jl @@ -43,6 +43,15 @@ Default tolerance for deciding to warn if the provided `A` is not hermitian. """ default_hermitian_tol(A) = eps(norm(A, Inf))^(3 / 4) +""" + default_domain_atol(λ) + +Default absolute tolerance for deciding when the eigenvalues `λ` should be considered +to lie outside of the domain of a matrix function, e.g. on the negative real axis for +[`squareroot`](@ref) and [`logarithm`](@ref) of a real matrix. +""" +default_domain_atol(λ) = defaulttol(λ) * maximum(abs, λ; init = abs(zero(eltype(λ)))) + const DEFAULT_FIXGAUGE = Ref(true) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl new file mode 100644 index 000000000..0625323c3 --- /dev/null +++ b/src/implementations/logarithm.jl @@ -0,0 +1,88 @@ +# Inputs +# ------ +function copy_input(::typeof(logarithm), A::AbstractMatrix) + return copy!(similar(A, float(eltype(A))), A) +end +copy_input(::typeof(logarithm), A::Diagonal) = Diagonal(float.(diagview(A))) + +function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, alg::AbstractAlgorithm) + m = LinearAlgebra.checksquare(A) + @check_size(logA, (m, m)) + @check_scalar(logA, A) + return nothing +end + +function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, ::DiagonalAlgorithm) + m = LinearAlgebra.checksquare(A) + @assert isdiag(A) + @assert logA isa Diagonal + @check_size(logA, (m, m)) + @check_scalar(logA, A) + return nothing +end + +# Algorithm selection +# ------------------- +logarithm!(A::AbstractMatrix, alg::DefaultAlgorithm) = logarithm!(A, select_algorithm(logarithm!, A, nothing; alg.kwargs...)) +logarithm!(A::AbstractMatrix, out, alg::DefaultAlgorithm) = logarithm!(A, out, select_algorithm(logarithm!, A, nothing; alg.kwargs...)) + +# Outputs +# ------- +initialize_output(::typeof(logarithm!), A::AbstractMatrix, ::AbstractAlgorithm) = A + +# Implementation +# -------------- +function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaLA) + check_input(logarithm!, A, logA, alg) + isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `logarithm`")) + result = LinearAlgebra.log(A) + _copy_result!(logarithm!, logA, result) + return logA +end + +function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh) + check_input(logarithm!, A, logA, alg) + D, V = eigh_full!(A, alg.eigh_alg) + λ = diagview(D) + atol = something(alg.domain_atol, default_domain_atol(λ)) + _check_nonzero_eigenvalues(λ, atol) + _clamp_domain_eigenvalues!(λ, atol) + λ .= log.(λ) + VD = V * D + mul!(logA, VD, V') + return project_hermitian!(logA) +end + +function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig) + check_input(logarithm!, A, logA, alg) + D, V = eig_full!(A, alg.eig_alg) + λ = diagview(D) + atol = something(alg.domain_atol, default_domain_atol(λ)) + _check_nonzero_eigenvalues(λ, atol) + if eltype(A) <: Real + _clamp_domain_eigenvalues!(λ, atol) + λ .= log.(λ) + VD = V * D + logAc = rdiv!(VD, LinearAlgebra.lu!(V)) + return logA .= real.(logAc) + else + λ .= log.(λ) + logA .= V .* transpose(λ) + return rdiv!(logA, LinearAlgebra.lu!(V)) + end +end + +# Diagonal logic +# -------------- +function logarithm!(A::AbstractMatrix, logA, alg::DiagonalAlgorithm) + check_input(logarithm!, A, logA, alg) + λ = diagview(logA) + copyto!(λ, diagview(A)) + atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) + _check_nonzero_eigenvalues(λ, atol) + if eltype(λ) <: Real + _clamp_domain_eigenvalues!(λ, atol) + end + λ .= log.(λ) + return logA +end diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl new file mode 100644 index 000000000..b6ba5a8e8 --- /dev/null +++ b/src/implementations/matrixfunctions.jl @@ -0,0 +1,91 @@ +# Shared helpers for matrix functions with a restricted domain +# ------------------------------------------------------------- + +# Clamp real eigenvalues that are negative within `atol` (rounding artifacts) to zero, +# and throw a `DomainError` for eigenvalues that are genuinely negative, since then the +# result cannot be expressed with the same (real) scalar type. +function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) + for i in eachindex(λ) + x = λ[i] + if x < -atol + throw( + DomainError( + x, + "The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " * + "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." + ) + ) + elseif x < 0 + λ[i] = zero(x) + end + end + return λ +end + +# Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative +# real axis obstruct a real result; complex-conjugate pairs do not. +function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real) + for i in eachindex(λ) + x = λ[i] + if abs(imag(x)) <= atol && real(x) < 0 + if real(x) < -atol + throw( + DomainError( + x, + "The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " * + "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." + ) + ) + else + λ[i] = zero(x) + end + end + end + return λ +end + +# Reject (numerically) zero eigenvalues for functions that are undefined there, +# e.g. `logarithm` and `power` with a negative fractional power. +function _check_nonzero_eigenvalues(λ, atol::Real) + for x in λ + if abs(x) <= atol + throw( + DomainError( + x, + "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined." + ) + ) + end + end + return λ +end + +# For `MatrixFunctionViaLA`, domain violations surface as a complex result from +# `LinearAlgebra`; detect this after the fact to preserve type-stable output. +# `LinearAlgebra` sometimes computes a real result in complex arithmetic (e.g. fractional +# powers via a complex Schur decomposition), so rounding-level imaginary components do +# not signal a domain violation and are simply discarded. +function _copy_result!(f, out::AbstractMatrix, result::AbstractMatrix) + if eltype(result) <: Complex && !(eltype(out) <: Complex) + # base the tolerance on the working precision, which is the lower of the two: + # `LinearAlgebra` may internally promote, e.g. `Float32` fractional powers are + # computed with `Float64`-eltype results but `Float32`-level accuracy + atol = max(defaulttol(out), defaulttol(result)) * norm(result, Inf) + for x in result + if abs(imag(x)) > atol + throw( + DomainError( + f, + "The result of this matrix function applied to the given real matrix is complex (eigenvalues on the negative real axis). " * + "Pass a complex matrix to obtain the principal value, or use `MatrixFunctionViaEigh`/`MatrixFunctionViaEig` with a suitable " * + "`domain_atol` if the offending eigenvalues are rounding artifacts." + ) + ) + end + end + out .= real.(result) + else + copy!(out, result) + end + return out +end diff --git a/src/implementations/power.jl b/src/implementations/power.jl new file mode 100644 index 000000000..c73a18dd3 --- /dev/null +++ b/src/implementations/power.jl @@ -0,0 +1,102 @@ +# Inputs +# ------ +function copy_input(::typeof(power), A::AbstractMatrix, p::Real) + return copy!(similar(A, float(eltype(A))), A), p +end +copy_input(::typeof(power), A::Diagonal, p::Real) = Diagonal(float.(diagview(A))), p + +function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, alg::AbstractAlgorithm) + m = LinearAlgebra.checksquare(A) + @check_size(powA, (m, m)) + @check_scalar(powA, A) + return nothing +end + +function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, ::DiagonalAlgorithm) + m = LinearAlgebra.checksquare(A) + @assert isdiag(A) + @assert powA isa Diagonal + @check_size(powA, (m, m)) + @check_scalar(powA, A) + return nothing +end + +# Algorithm selection +# ------------------- +power!(A::AbstractMatrix, p::Real, alg::DefaultAlgorithm) = power!(A, p, select_algorithm(power!, (A, p), nothing; alg.kwargs...)) +power!(A::AbstractMatrix, p::Real, out, alg::DefaultAlgorithm) = power!(A, p, out, select_algorithm(power!, (A, p), nothing; alg.kwargs...)) + +# Outputs +# ------- +initialize_output(::typeof(power!), A::AbstractMatrix, p::Real, ::AbstractAlgorithm) = A + +# Implementation +# -------------- +function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaLA) + check_input(power!, A, p, powA, alg) + isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `power`")) + result = A^p + _copy_result!(power!, powA, result) + return powA +end + +function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) + check_input(power!, A, p, powA, alg) + D, V = eigh_full!(A, alg.eigh_alg) + λ = diagview(D) + if isinteger(p) + p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) + λ .= λ .^ p + VD = V * D + mul!(powA, VD, V') + else + atol = something(alg.domain_atol, default_domain_atol(λ)) + p < 0 && _check_nonzero_eigenvalues(λ, atol) + _clamp_domain_eigenvalues!(λ, atol) + # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction + λ .= λ .^ (p / 2) + Vs = rmul!(V, D) + mul!(powA, Vs, Vs') + end + return project_hermitian!(powA) +end + +function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) + check_input(power!, A, p, powA, alg) + D, V = eig_full!(A, alg.eig_alg) + λ = diagview(D) + if isinteger(p) + p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) + else + atol = something(alg.domain_atol, default_domain_atol(λ)) + p < 0 && _check_nonzero_eigenvalues(λ, atol) + eltype(A) <: Real && _clamp_domain_eigenvalues!(λ, atol) + end + if eltype(A) <: Real + λ .= λ .^ p + VD = V * D + powAc = rdiv!(VD, LinearAlgebra.lu!(V)) + return powA .= real.(powAc) + else + λ .= λ .^ p + powA .= V .* transpose(λ) + return rdiv!(powA, LinearAlgebra.lu!(V)) + end +end + +# Diagonal logic +# -------------- +function power!(A::AbstractMatrix, p::Real, powA, alg::DiagonalAlgorithm) + check_input(power!, A, p, powA, alg) + λ = diagview(powA) + copyto!(λ, diagview(A)) + if isinteger(p) + p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) + else + atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) + p < 0 && _check_nonzero_eigenvalues(λ, atol) + eltype(λ) <: Real && _clamp_domain_eigenvalues!(λ, atol) + end + λ .= λ .^ p + return powA +end diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl new file mode 100644 index 000000000..fe687e2c4 --- /dev/null +++ b/src/implementations/squareroot.jl @@ -0,0 +1,86 @@ +# Inputs +# ------ +function copy_input(::typeof(squareroot), A::AbstractMatrix) + return copy!(similar(A, float(eltype(A))), A) +end +copy_input(::typeof(squareroot), A::Diagonal) = Diagonal(float.(diagview(A))) + +function check_input(::typeof(squareroot!), A::AbstractMatrix, sqrtA, alg::AbstractAlgorithm) + m = LinearAlgebra.checksquare(A) + @check_size(sqrtA, (m, m)) + @check_scalar(sqrtA, A) + return nothing +end + +function check_input(::typeof(squareroot!), A::AbstractMatrix, sqrtA, ::DiagonalAlgorithm) + m = LinearAlgebra.checksquare(A) + @assert isdiag(A) + @assert sqrtA isa Diagonal + @check_size(sqrtA, (m, m)) + @check_scalar(sqrtA, A) + return nothing +end + +# Algorithm selection +# ------------------- +squareroot!(A::AbstractMatrix, alg::DefaultAlgorithm) = squareroot!(A, select_algorithm(squareroot!, A, nothing; alg.kwargs...)) +squareroot!(A::AbstractMatrix, out, alg::DefaultAlgorithm) = squareroot!(A, out, select_algorithm(squareroot!, A, nothing; alg.kwargs...)) + +# Outputs +# ------- +initialize_output(::typeof(squareroot!), A::AbstractMatrix, ::AbstractAlgorithm) = A + +# Implementation +# -------------- +function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaLA) + check_input(squareroot!, A, sqrtA, alg) + isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `squareroot`")) + result = LinearAlgebra.sqrt(A) + _copy_result!(squareroot!, sqrtA, result) + return sqrtA +end + +function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEigh) + check_input(squareroot!, A, sqrtA, alg) + D, V = eigh_full!(A, alg.eigh_alg) + λ = diagview(D) + atol = something(alg.domain_atol, default_domain_atol(λ)) + _clamp_domain_eigenvalues!(λ, atol) + # `sqrt(A) = (V * D^(1/4)) * (V * D^(1/4))'` is hermitian by construction + λ .= sqrt.(sqrt.(λ)) + Vs = rmul!(V, D) + mul!(sqrtA, Vs, Vs') + return project_hermitian!(sqrtA) +end + +function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig) + check_input(squareroot!, A, sqrtA, alg) + D, V = eig_full!(A, alg.eig_alg) + λ = diagview(D) + if eltype(A) <: Real + atol = something(alg.domain_atol, default_domain_atol(λ)) + _clamp_domain_eigenvalues!(λ, atol) + λ .= sqrt.(λ) + VD = V * D + sqrtAc = rdiv!(VD, LinearAlgebra.lu!(V)) + return sqrtA .= real.(sqrtAc) + else + λ .= sqrt.(λ) + sqrtA .= V .* transpose(λ) + return rdiv!(sqrtA, LinearAlgebra.lu!(V)) + end +end + +# Diagonal logic +# -------------- +function squareroot!(A::AbstractMatrix, sqrtA, alg::DiagonalAlgorithm) + check_input(squareroot!, A, sqrtA, alg) + λ = diagview(sqrtA) + copyto!(λ, diagview(A)) + if eltype(λ) <: Real + atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) + _clamp_domain_eigenvalues!(λ, atol) + end + λ .= sqrt.(λ) + return sqrtA +end diff --git a/src/interface/logarithm.jl b/src/interface/logarithm.jl new file mode 100644 index 000000000..39acb133c --- /dev/null +++ b/src/interface/logarithm.jl @@ -0,0 +1,41 @@ +# Logarithm +# --------- + +""" + logarithm(A; kwargs...) -> logA + logarithm(A, alg::AbstractAlgorithm) -> logA + logarithm!(A, [logA]; kwargs...) -> logA + logarithm!(A, [logA], alg::AbstractAlgorithm) -> logA + +Compute the principal logarithm `logA` of the square matrix `A`, i.e. the logarithm +whose eigenvalues have imaginary part in `(-π, π]`. + +The scalar type of the output matches that of the input. +As a consequence, a real matrix with eigenvalues on the negative real axis, for which +the principal logarithm is complex, leads to a `DomainError`; pass a complex matrix +to obtain the principal value. +A matrix with (numerically) zero eigenvalues has no logarithm and also leads to a +`DomainError`. +Both checks use a tolerance `domain_atol`, which can be specified for the algorithms +that support it and defaults to [`default_domain_atol`](@ref). + +!!! note + The bang method `logarithm!` optionally accepts the output structure and + possibly destroys the input matrix `A`. Always use the return value of the function + as it may not always be possible to use the provided `logA` as output. +""" +@functiondef logarithm + +# Algorithm selection +# ------------------- +default_logarithm_algorithm(A; kwargs...) = default_logarithm_algorithm(typeof(A); kwargs...) +function default_logarithm_algorithm(T::Type; kwargs...) + return MatrixFunctionViaLA(; kwargs...) +end +function default_logarithm_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} + return DiagonalAlgorithm(; kwargs...) +end + +function default_algorithm(::typeof(logarithm!), ::Type{A}; kwargs...) where {A} + return default_logarithm_algorithm(A; kwargs...) +end diff --git a/src/interface/matrixfunctions.jl b/src/interface/matrixfunctions.jl index b0ea57da8..a0e0a7304 100644 --- a/src/interface/matrixfunctions.jl +++ b/src/interface/matrixfunctions.jl @@ -1,10 +1,10 @@ # ================================ -# EXPONENTIAL ALGORITHMS +# MATRIX FUNCTION ALGORITHMS # ================================ """ MatrixFunctionViaLA() -Algorithm type to denote finding the exponential of `A` via the implementation of `LinearAlgebra`. +Algorithm type to denote computing a function of a matrix `A` via the implementation of `LinearAlgebra`. """ @algdef MatrixFunctionViaLA @@ -29,31 +29,49 @@ As this algorithm requires no LAPACK support, it also applies at arbitrary preci @algdef MatrixFunctionViaTaylor """ - MatrixFunctionViaEigh(eigh_alg) + MatrixFunctionViaEigh(eigh_alg; domain_atol=nothing) Algorithm type for computing a function of a matrix by computing its hermitian eigenvalue decomposition and applying the function to the eigenvalues. The `eigh_alg` specifies which hermitian eigendecomposition implementation to use. +For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`logarithm`](@ref)), +`domain_atol` specifies the absolute tolerance below which out-of-domain eigenvalues are treated +as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default +tolerance [`default_domain_atol`](@ref). """ -struct MatrixFunctionViaEigh{A <: AbstractAlgorithm} <: AbstractAlgorithm +struct MatrixFunctionViaEigh{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm eigh_alg::A + domain_atol::T +end +function MatrixFunctionViaEigh(eigh_alg::AbstractAlgorithm; domain_atol::Union{Nothing, Real} = nothing) + return MatrixFunctionViaEigh(eigh_alg, domain_atol) end function Base.show(io::IO, alg::MatrixFunctionViaEigh) print(io, "MatrixFunctionViaEigh(") _show_alg(io, alg.eigh_alg) + isnothing(alg.domain_atol) || print(io, "; domain_atol=", alg.domain_atol) return print(io, ")") end """ - MatrixFunctionViaEig(eig_alg) + MatrixFunctionViaEig(eig_alg; domain_atol=nothing) Algorithm type for computing a function of a matrix by computing its eigenvalue decomposition and applying the function to the eigenvalues. The `eig_alg` specifies which eigendecomposition implementation to use. +For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`logarithm`](@ref)), +`domain_atol` specifies the absolute tolerance below which out-of-domain eigenvalues are treated +as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default +tolerance [`default_domain_atol`](@ref). """ -struct MatrixFunctionViaEig{A <: AbstractAlgorithm} <: AbstractAlgorithm +struct MatrixFunctionViaEig{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm eig_alg::A + domain_atol::T +end +function MatrixFunctionViaEig(eig_alg::AbstractAlgorithm; domain_atol::Union{Nothing, Real} = nothing) + return MatrixFunctionViaEig(eig_alg, domain_atol) end function Base.show(io::IO, alg::MatrixFunctionViaEig) print(io, "MatrixFunctionViaEig(") _show_alg(io, alg.eig_alg) + isnothing(alg.domain_atol) || print(io, "; domain_atol=", alg.domain_atol) return print(io, ")") end diff --git a/src/interface/power.jl b/src/interface/power.jl new file mode 100644 index 000000000..6da9b42e8 --- /dev/null +++ b/src/interface/power.jl @@ -0,0 +1,42 @@ +# Power +# ----- + +""" + power(A, p::Real; kwargs...) -> powA + power(A, p::Real, alg::AbstractAlgorithm) -> powA + power!(A, p::Real, [powA]; kwargs...) -> powA + power!(A, p::Real, [powA], alg::AbstractAlgorithm) -> powA + +Compute the matrix power `powA = A^p` of the square matrix `A`. +For integer `p` this is defined for any square matrix (invertible if `p < 0`); +for fractional `p` the principal power `exp(p * log(A))` is computed. + +The scalar type of the output matches that of the input. +As a consequence, for fractional `p`, a real matrix with eigenvalues on the negative +real axis, for which the principal power is complex, leads to a `DomainError`; pass a +complex matrix to obtain the principal value. +Real eigenvalues that are negative within a tolerance `domain_atol` are treated as +rounding artifacts and clamped to zero, where `domain_atol` can be specified for the +algorithms that support it and defaults to [`default_domain_atol`](@ref). +For negative fractional `p`, (numerically) zero eigenvalues also lead to a `DomainError`. + +!!! note + The bang method `power!` optionally accepts the output structure and + possibly destroys the input matrix `A`. Always use the return value of the function + as it may not always be possible to use the provided `powA` as output. +""" +@functiondef n_args = 2 power + +# Algorithm selection +# ------------------- +default_power_algorithm(A; kwargs...) = default_power_algorithm(typeof(A); kwargs...) +function default_power_algorithm(T::Type; kwargs...) + return MatrixFunctionViaLA(; kwargs...) +end +function default_power_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} + return DiagonalAlgorithm(; kwargs...) +end + +function default_algorithm(::typeof(power!), ::Tuple{A, P}; kwargs...) where {A, P} + return default_power_algorithm(A; kwargs...) +end diff --git a/src/interface/squareroot.jl b/src/interface/squareroot.jl new file mode 100644 index 000000000..295cb4cbe --- /dev/null +++ b/src/interface/squareroot.jl @@ -0,0 +1,40 @@ +# Square root +# ----------- + +""" + squareroot(A; kwargs...) -> sqrtA + squareroot(A, alg::AbstractAlgorithm) -> sqrtA + squareroot!(A, [sqrtA]; kwargs...) -> sqrtA + squareroot!(A, [sqrtA], alg::AbstractAlgorithm) -> sqrtA + +Compute the principal square root `sqrtA` of the square matrix `A`, i.e. the square root +whose eigenvalues have nonnegative real part. + +The scalar type of the output matches that of the input. +As a consequence, a real matrix with eigenvalues on the negative real axis, for which +the principal square root is complex, leads to a `DomainError`; pass a complex matrix +to obtain the principal value. +Real eigenvalues that are negative within a tolerance `domain_atol` are treated as +rounding artifacts and clamped to zero, where `domain_atol` can be specified for the +algorithms that support it and defaults to [`default_domain_atol`](@ref). + +!!! note + The bang method `squareroot!` optionally accepts the output structure and + possibly destroys the input matrix `A`. Always use the return value of the function + as it may not always be possible to use the provided `sqrtA` as output. +""" +@functiondef squareroot + +# Algorithm selection +# ------------------- +default_squareroot_algorithm(A; kwargs...) = default_squareroot_algorithm(typeof(A); kwargs...) +function default_squareroot_algorithm(T::Type; kwargs...) + return MatrixFunctionViaLA(; kwargs...) +end +function default_squareroot_algorithm(::Type{T}; kwargs...) where {T <: Diagonal} + return DiagonalAlgorithm(; kwargs...) +end + +function default_algorithm(::typeof(squareroot!), ::Type{A}; kwargs...) where {A} + return default_squareroot_algorithm(A; kwargs...) +end From 080ee51517d1e7181bad06bf59622bc8c4bab15e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 16:23:13 -0400 Subject: [PATCH 02/22] test: add tests for `squareroot`, `logarithm` and `power` Mirrors the `exponential` test layout: BLAS floats against the LinearAlgebra reference, hermitian fast paths, domain-error and tolerance-clamping edge cases, `Diagonal` inputs over generic float types, and BigFloat coverage via the GenericSchur and GenericLinearAlgebra test directories. --- test/genericlinearalgebra/logarithm.jl | 28 ++++++ test/genericlinearalgebra/power.jl | 29 ++++++ test/genericlinearalgebra/squareroot.jl | 30 ++++++ test/genericschur/logarithm.jl | 32 ++++++ test/genericschur/power.jl | 31 ++++++ test/genericschur/squareroot.jl | 32 ++++++ test/logarithm.jl | 102 +++++++++++++++++++ test/power.jl | 126 ++++++++++++++++++++++++ test/squareroot.jl | 106 ++++++++++++++++++++ 9 files changed, 516 insertions(+) create mode 100644 test/genericlinearalgebra/logarithm.jl create mode 100644 test/genericlinearalgebra/power.jl create mode 100644 test/genericlinearalgebra/squareroot.jl create mode 100644 test/genericschur/logarithm.jl create mode 100644 test/genericschur/power.jl create mode 100644 test/genericschur/squareroot.jl create mode 100644 test/logarithm.jl create mode 100644 test/power.jl create mode 100644 test/squareroot.jl diff --git a/test/genericlinearalgebra/logarithm.jl b/test/genericlinearalgebra/logarithm.jl new file mode 100644 index 000000000..4cd69d342 --- /dev/null +++ b/test/genericlinearalgebra/logarithm.jl @@ -0,0 +1,28 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using GenericLinearAlgebra + +GenericFloats = (BigFloat, Complex{BigFloat}) + +@testset "logarithm! for T = $T" for T in GenericFloats + rng = StableRNG(123) + m = 24 + + X = randn(rng, T, m, m) + A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I + D, V = @constinferred eigh_full(A) + + algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) + @testset "algorithm $alg" for alg in algs + logA = @constinferred logarithm!(copy(A); alg) + logA2 = @constinferred logarithm(A; alg) + @test logA2 ≈ logA + + Dlog, Vlog = @constinferred eigh_full(logA) + @test diagview(Dlog) ≈ log.(diagview(D)) + end +end diff --git a/test/genericlinearalgebra/power.jl b/test/genericlinearalgebra/power.jl new file mode 100644 index 000000000..1170be498 --- /dev/null +++ b/test/genericlinearalgebra/power.jl @@ -0,0 +1,29 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using GenericLinearAlgebra + +GenericFloats = (BigFloat, Complex{BigFloat}) + +@testset "power! for T = $T" for T in GenericFloats + rng = StableRNG(123) + m = 24 + + X = randn(rng, T, m, m) + A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I + D, V = @constinferred eigh_full(A) + + algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) + @testset "algorithm $alg" for alg in algs + @test power(A, 3, alg) ≈ A * A * A + @test power(A, -1, alg) * A ≈ LinearAlgebra.I + @test power(A, big"0.5", alg) ≈ squareroot(A; alg) + + powA = @constinferred power(A, big"0.5", alg) + Dpow, Vpow = @constinferred eigh_full(powA) + @test diagview(Dpow) ≈ sqrt.(diagview(D)) + end +end diff --git a/test/genericlinearalgebra/squareroot.jl b/test/genericlinearalgebra/squareroot.jl new file mode 100644 index 000000000..1f9ee93af --- /dev/null +++ b/test/genericlinearalgebra/squareroot.jl @@ -0,0 +1,30 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using GenericLinearAlgebra + +GenericFloats = (BigFloat, Complex{BigFloat}) + +@testset "squareroot! for T = $T" for T in GenericFloats + rng = StableRNG(123) + m = 24 + + X = randn(rng, T, m, m) + A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I + D, V = @constinferred eigh_full(A) + + algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) + @testset "algorithm $alg" for alg in algs + sqrtA = @constinferred squareroot!(copy(A); alg) + sqrtA2 = @constinferred squareroot(A; alg) + @test sqrtA2 ≈ sqrtA + @test sqrtA * sqrtA ≈ A + @test LinearAlgebra.ishermitian(sqrtA) + + Dsqrt, Vsqrt = @constinferred eigh_full(sqrtA) + @test diagview(Dsqrt) ≈ sqrt.(diagview(D)) + end +end diff --git a/test/genericschur/logarithm.jl b/test/genericschur/logarithm.jl new file mode 100644 index 000000000..aeaf2c7d4 --- /dev/null +++ b/test/genericschur/logarithm.jl @@ -0,0 +1,32 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using GenericSchur + +GenericFloats = (BigFloat, Complex{BigFloat}) + +@testset "logarithm! for T = $T" for T in GenericFloats + rng = StableRNG(123) + m = 24 + + # spectrum inside a disk around 1, away from the negative real axis and zero + A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) + D, V = @constinferred eig_full(A) + + logA = @constinferred logarithm(A) + @test eltype(logA) == T + @test exponential(logA) ≈ A + + algs = (MatrixFunctionViaEig(GS_QRIteration()),) + @testset "algorithm $alg" for alg in algs + logA2 = @constinferred logarithm(A; alg) + @test logA2 ≈ logA + + Dlog, Vlog = @constinferred eig_full(logA2) + by = x -> (real(x), imag(x)) + @test sort(diagview(Dlog); by) ≈ sort(log.(diagview(D)); by) + end +end diff --git a/test/genericschur/power.jl b/test/genericschur/power.jl new file mode 100644 index 000000000..f4b4a907b --- /dev/null +++ b/test/genericschur/power.jl @@ -0,0 +1,31 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using GenericSchur + +GenericFloats = (BigFloat, Complex{BigFloat}) + +@testset "power! for T = $T" for T in GenericFloats + rng = StableRNG(123) + m = 24 + + # spectrum inside a disk around 1, away from the negative real axis and zero + A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) + + powA = @constinferred power(A, 3) + @test powA ≈ A * A * A + @test eltype(powA) == T + + powA = @constinferred power(A, big"0.5") + @test powA ≈ squareroot(A) + + algs = (MatrixFunctionViaEig(GS_QRIteration()),) + @testset "algorithm $alg" for alg in algs + @test power(A, 3, alg) ≈ A * A * A + @test power(A, -1, alg) * A ≈ LinearAlgebra.I + @test power(A, big"0.5", alg) ≈ squareroot(A) + end +end diff --git a/test/genericschur/squareroot.jl b/test/genericschur/squareroot.jl new file mode 100644 index 000000000..33a2140e4 --- /dev/null +++ b/test/genericschur/squareroot.jl @@ -0,0 +1,32 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using GenericSchur + +GenericFloats = (BigFloat, Complex{BigFloat}) + +@testset "squareroot! for T = $T" for T in GenericFloats + rng = StableRNG(123) + m = 24 + + # spectrum inside a disk around 1, away from the negative real axis + A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) + D, V = @constinferred eig_full(A) + + sqrtA = @constinferred squareroot(A) + @test sqrtA * sqrtA ≈ A + @test eltype(sqrtA) == T + + algs = (MatrixFunctionViaEig(GS_QRIteration()),) + @testset "algorithm $alg" for alg in algs + sqrtA2 = @constinferred squareroot(A; alg) + @test sqrtA2 ≈ sqrtA + + Dsqrt, Vsqrt = @constinferred eig_full(sqrtA2) + by = x -> (real(x), imag(x)) + @test sort(diagview(Dsqrt); by) ≈ sort(sqrt.(diagview(D)); by) + end +end diff --git a/test/logarithm.jl b/test/logarithm.jl new file mode 100644 index 000000000..2307141a5 --- /dev/null +++ b/test/logarithm.jl @@ -0,0 +1,102 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra +using LinearAlgebra: exp + +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) + +@testset "logarithm! for T = $T" for T in BLASFloats + rng = StableRNG(123) + m = 54 + + # spectrum inside a disk around 1, away from the negative real axis and zero + A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) + Ac = copy(A) + logA = LinearAlgebra.log(A) + + logA2 = @constinferred logarithm(A) + @test logA ≈ logA2 + @test exp(logA2) ≈ A + @test A == Ac + + algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) + @testset "algorithm $alg" for alg in algs + logA2 = @constinferred logarithm(A, alg) + @test logA ≈ logA2 + @test A == Ac + end + + @test_throws DomainError logarithm(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) + + # roundtrip with exponential + @test logarithm(exponential(logA)) ≈ logA +end + +@testset "logarithm! for hermitian T = $T" for T in BLASFloats + rng = StableRNG(123) + m = 54 + + X = randn(rng, T, m, m) + A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) + Ac = copy(A) + logA = LinearAlgebra.log(LinearAlgebra.Hermitian(A)) + + algs = ( + MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()), + MatrixFunctionViaEig(LAPACK_Simple()), + ) + @testset "algorithm $alg" for alg in algs + logA2 = @constinferred logarithm(A, alg) + @test logA ≈ logA2 + @test A == Ac + end +end + +@testset "logarithm! domain handling for T = $T" for T in (Float32, Float64) + rng = StableRNG(123) + m = 4 + + X = randn(rng, T, m, m) + V = Matrix(LinearAlgebra.qr(X).Q) + A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V' + A = (A + A') / 2 + + # negative eigenvalue: DomainError for real input, complex principal value otherwise + @test_throws DomainError logarithm(A) + @test_throws DomainError logarithm(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) + @test_throws DomainError logarithm(A; alg = MatrixFunctionViaEig(LAPACK_Simple())) + + logA = @constinferred logarithm(complex.(A)) + @test exp(logA) ≈ A + + # (numerically) zero eigenvalue: no logarithm exists + Asing = V * LinearAlgebra.Diagonal(T[0, 1, 2, 3]) * V' + Asing = (Asing + Asing') / 2 + @test_throws DomainError logarithm(Asing; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) + @test_throws DomainError logarithm(Asing; alg = MatrixFunctionViaEig(LAPACK_Simple())) + @test_throws DomainError logarithm(complex.(Asing); alg = MatrixFunctionViaEig(LAPACK_Simple())) +end + +@testset "logarithm! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) + rng = StableRNG(123) + m = 54 + + data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) + A = Diagonal(data) + Ac = copy(A) + logA = LinearAlgebra.log(A) + + logA2 = @constinferred logarithm(A) + @test logA2 isa Diagonal + @test logA ≈ logA2 + @test A == Ac + + if T <: Real + @test_throws DomainError logarithm(Diagonal(-data)) + end + @test_throws DomainError logarithm(Diagonal(zeros(T, m))) +end diff --git a/test/power.jl b/test/power.jl new file mode 100644 index 000000000..99cf659b1 --- /dev/null +++ b/test/power.jl @@ -0,0 +1,126 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra + +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) + +@testset "power! for T = $T" for T in BLASFloats + rng = StableRNG(123) + m = 54 + + # spectrum inside a disk around 1, away from the negative real axis and zero + A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) + Ac = copy(A) + + # integer powers + @testset "integer p = $p" for p in (0, 1, 3, -2) + powA = A^p + powA2 = @constinferred power(A, p) + @test powA ≈ powA2 + @test A == Ac + + algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) + @testset "algorithm $alg" for alg in algs + powA2 = @constinferred power(A, p, alg) + @test powA ≈ powA2 + @test A == Ac + end + end + + # fractional powers + @testset "fractional p = $p" for p in (0.5, 0.75, -0.25, 3.5) + powA = A^p + powA2 = @constinferred power(A, p) + @test powA ≈ powA2 + @test A == Ac + + algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) + @testset "algorithm $alg" for alg in algs + powA2 = @constinferred power(A, p, alg) + @test powA ≈ powA2 + @test A == Ac + end + end + + @test power(A, 0.5) ≈ squareroot(A) + @test power(A, 3.0) ≈ power(A, 3) + @test_throws DomainError power(A, 0.5; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) +end + +@testset "power! for hermitian T = $T" for T in BLASFloats + rng = StableRNG(123) + m = 54 + + X = randn(rng, T, m, m) + A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) + Ac = copy(A) + + algs = ( + MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()), + MatrixFunctionViaEig(LAPACK_Simple()), + ) + @testset "p = $p, algorithm $alg" for p in (2, -1, 0.5, -0.5), alg in algs + powA = LinearAlgebra.Hermitian(A)^p + powA2 = @constinferred power(A, p, alg) + @test powA ≈ powA2 + @test A == Ac + end + + @test LinearAlgebra.ishermitian(power(A, 0.5, MatrixFunctionViaEigh(LAPACK_QRIteration()))) +end + +@testset "power! domain handling for T = $T" for T in (Float32, Float64) + rng = StableRNG(123) + m = 4 + + X = randn(rng, T, m, m) + V = Matrix(LinearAlgebra.qr(X).Q) + A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V' + A = (A + A') / 2 + + # integer powers of matrices with negative eigenvalues are fine + @test power(A, 3) ≈ A * A * A + @test power(A, 3, MatrixFunctionViaEigh(LAPACK_QRIteration())) ≈ A * A * A + @test power(A, 3, MatrixFunctionViaEig(LAPACK_Simple())) ≈ A * A * A + + # fractional powers require staying off the negative real axis for real input + @test_throws DomainError power(A, 0.5) + @test_throws DomainError power(A, 0.5, MatrixFunctionViaEigh(LAPACK_QRIteration())) + @test_throws DomainError power(A, 0.5, MatrixFunctionViaEig(LAPACK_Simple())) + powA = @constinferred power(complex.(A), 0.5) + @test powA * powA ≈ A + + # (numerically) singular matrices with negative fractional powers + Asing = V * LinearAlgebra.Diagonal(T[0, 1, 2, 3]) * V' + Asing = (Asing + Asing') / 2 + @test_throws DomainError power(Asing, -0.5, MatrixFunctionViaEigh(LAPACK_QRIteration())) + @test_throws DomainError power(Asing, -0.5, MatrixFunctionViaEig(LAPACK_Simple())) +end + +@testset "power! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) + rng = StableRNG(123) + m = 54 + + data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) + A = Diagonal(data) + Ac = copy(A) + + @testset "p = $p" for p in (2, -1, 0.5) + powA = A^p + powA2 = @constinferred power(A, p) + @test powA2 isa Diagonal + @test powA ≈ powA2 + @test A == Ac + end + + if T <: Real + @test_throws DomainError power(Diagonal(-data), 0.5) + @test power(Diagonal(-data), 2) ≈ Diagonal(data .^ 2) + end + @test_throws SingularException power(Diagonal(zeros(T, m)), -1) + @test_throws DomainError power(Diagonal(zeros(T, m)), -0.5) +end diff --git a/test/squareroot.jl b/test/squareroot.jl new file mode 100644 index 000000000..0d02a2f5a --- /dev/null +++ b/test/squareroot.jl @@ -0,0 +1,106 @@ +using MatrixAlgebraKit +using Test +using TestExtras +using StableRNGs +using MatrixAlgebraKit: diagview +using LinearAlgebra + +BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) +GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) + +@testset "squareroot! for T = $T" for T in BLASFloats + rng = StableRNG(123) + m = 54 + + # spectrum inside a disk around 1, away from the negative real axis + A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) + Ac = copy(A) + sqrtA = LinearAlgebra.sqrt(A) + + sqrtA2 = @constinferred squareroot(A) + @test sqrtA ≈ sqrtA2 + @test sqrtA2 * sqrtA2 ≈ A + @test A == Ac + + algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) + @testset "algorithm $alg" for alg in algs + sqrtA2 = @constinferred squareroot(A, alg) + @test sqrtA ≈ sqrtA2 + @test A == Ac + end + + @test_throws DomainError squareroot(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) +end + +@testset "squareroot! for hermitian T = $T" for T in BLASFloats + rng = StableRNG(123) + m = 54 + + X = randn(rng, T, m, m) + A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) + Ac = copy(A) + sqrtA = LinearAlgebra.sqrt(LinearAlgebra.Hermitian(A)) + + algs = ( + MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()), + MatrixFunctionViaEig(LAPACK_Simple()), + ) + @testset "algorithm $alg" for alg in algs + sqrtA2 = @constinferred squareroot(A, alg) + @test sqrtA ≈ sqrtA2 + @test A == Ac + end + + @test LinearAlgebra.ishermitian(squareroot(A, MatrixFunctionViaEigh(LAPACK_QRIteration()))) +end + +@testset "squareroot! domain handling for T = $T" for T in (Float32, Float64) + rng = StableRNG(123) + m = 4 + + # genuinely negative eigenvalue: DomainError for real input, complex principal value otherwise + X = randn(rng, T, m, m) + V = Matrix(LinearAlgebra.qr(X).Q) + A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V' + A = (A + A') / 2 + + @test_throws DomainError squareroot(A) + @test_throws DomainError squareroot(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) + @test_throws DomainError squareroot(A; alg = MatrixFunctionViaEig(LAPACK_Simple())) + + sqrtA = @constinferred squareroot(complex.(A)) + @test sqrtA * sqrtA ≈ A + + # roundoff-scale negative eigenvalue: clamped to zero instead of throwing + Aclamp = V * LinearAlgebra.Diagonal(T[-10 * eps(T), 1, 2, 3]) * V' + Aclamp = (Aclamp + Aclamp') / 2 + @testset "algorithm $alg" for alg in ( + MatrixFunctionViaEigh(LAPACK_QRIteration()), + MatrixFunctionViaEig(LAPACK_Simple()), + ) + sqrtA2 = @constinferred squareroot(Aclamp, alg) + @test eltype(sqrtA2) == T + @test sqrtA2 * sqrtA2 ≈ Aclamp atol = sqrt(eps(T)) + end +end + +@testset "squareroot! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) + rng = StableRNG(123) + m = 54 + + data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : randn(rng, T, m) + A = Diagonal(data) + Ac = copy(A) + sqrtA = LinearAlgebra.sqrt(A) + + sqrtA2 = @constinferred squareroot(A) + @test sqrtA2 isa Diagonal + @test sqrtA ≈ sqrtA2 + @test A == Ac + + if T <: Real + @test_throws DomainError squareroot(Diagonal(-data)) + # roundoff-scale negative entries are clamped + @test squareroot(Diagonal(T[-eps(T), 1])) ≈ Diagonal(T[0, 1]) + end +end From ebe7d4c9e2bfc67c6ced25173b689b520f8ce8c5 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 16:23:14 -0400 Subject: [PATCH 03/22] docs: document new matrix functions and domain handling --- docs/src/changelog.md | 4 +++ docs/src/user_interface/algorithms.md | 6 ++-- docs/src/user_interface/matrix_functions.md | 39 +++++++++++++++++++++ 3 files changed, 46 insertions(+), 3 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 58f704b6e..18e579c0b 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -22,6 +22,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added +- New matrix functions `squareroot`, `logarithm` and `power` (with both integer and fractional exponents), supporting the `MatrixFunctionViaLA`, `MatrixFunctionViaEig`, `MatrixFunctionViaEigh` and `DiagonalAlgorithm` algorithms. +- The scalar type of the output of these matrix functions matches that of the input, and out-of-domain eigenvalues (e.g. on the negative real axis for a real input) throw a `DomainError`; eigenvalues that violate the domain within a tolerance `domain_atol` (defaulting to `default_domain_atol`) are treated as rounding artifacts and clamped to the domain boundary. +- `MatrixFunctionViaEig` and `MatrixFunctionViaEigh` accept a `domain_atol` keyword argument to control this tolerance. + ### Changed ### Deprecated diff --git a/docs/src/user_interface/algorithms.md b/docs/src/user_interface/algorithms.md index d90b76074..dcb93f9ad 100644 --- a/docs/src/user_interface/algorithms.md +++ b/docs/src/user_interface/algorithms.md @@ -103,9 +103,9 @@ The following algorithms for matrix functions are available. | Algorithm | Applicable matrix functions | Key keyword arguments | |:----------|:--------------------------|:----------------------| | [`MatrixFunctionViaTaylor`](@ref) | exponential | `tol`, `balance` | -| [`MatrixFunctionViaLA`](@ref) | exponential | | -| [`MatrixFunctionViaEig`](@ref) | exponential | `eig_alg` | -| [`MatrixFunctionViaEigh`](@ref) | exponential | `eigh_alg` | +| [`MatrixFunctionViaLA`](@ref) | exponential, squareroot, logarithm, power | | +| [`MatrixFunctionViaEig`](@ref) | exponential, squareroot, logarithm, power | `eig_alg`, `domain_atol` | +| [`MatrixFunctionViaEigh`](@ref) | exponential, squareroot, logarithm, power | `eigh_alg`, `domain_atol` | For full docstring details on each algorithm type, see the corresponding section in [Decompositions](@ref). diff --git a/docs/src/user_interface/matrix_functions.md b/docs/src/user_interface/matrix_functions.md index 6153f2e6a..71570ba3c 100644 --- a/docs/src/user_interface/matrix_functions.md +++ b/docs/src/user_interface/matrix_functions.md @@ -37,3 +37,42 @@ MatrixAlgebraKit.MatrixFunctionViaLA MatrixAlgebraKit.MatrixFunctionViaEig MatrixAlgebraKit.MatrixFunctionViaEigh ``` + +## Domain considerations + +The functions below ([`squareroot`](@ref), [`logarithm`](@ref) and [`power`](@ref) with fractional powers) are only defined for matrices whose eigenvalues avoid (part of) the negative real axis, and their principal values are complex whenever eigenvalues on that axis are present. +In MatrixAlgebraKit, we aim to keep type stability, and thus the scalar type of the output always matches that of the input. +As such, a real matrix with eigenvalues on the negative real axis leads to a `DomainError`, you should pass a complex matrix instead to obtain the complex principal value. +To avoid spurious errors for eigenvalues that lie on the negative real axis only because of rounding errors (e.g. a positive semidefinite matrix with a tiny negative eigenvalue), eigenvalues within an absolute tolerance `domain_atol` of the domain boundary are clamped onto it. +This tolerance defaults to [`default_domain_atol`](@ref) and can be specified explicitly for the algorithms that support it, e.g. `MatrixFunctionViaEigh(eigh_alg; domain_atol=...)`. + +```@docs; canonical=false +MatrixAlgebraKit.default_domain_atol +``` + +## Square root + +The principal [square root](https://en.wikipedia.org/wiki/Square_root_of_a_matrix) of a square matrix `A` is the unique square root whose eigenvalues have nonnegative real part. +It is computed by the function [`squareroot`](@ref), where the default algorithm [`MatrixFunctionViaLA`](@ref) wraps the Schur-based implementation of `LinearAlgebra`, and the eigenvalue-decomposition-based algorithms [`MatrixFunctionViaEig`](@ref) and [`MatrixFunctionViaEigh`](@ref) are available as well. + +```@docs; canonical=false +squareroot +``` + +## Logarithm + +The principal [logarithm](https://en.wikipedia.org/wiki/Logarithm_of_a_matrix) of a square matrix `A` is the unique logarithm whose eigenvalues have imaginary part in `(-π, π]`, and exists for matrices without (numerically) zero eigenvalues that satisfy the domain considerations above. +It is computed by the function [`logarithm`](@ref), with the same algorithm choices as [`squareroot`](@ref). + +```@docs; canonical=false +logarithm +``` + +## Power + +Matrix powers `A^p` for real `p` are computed by the function [`power`](@ref), which takes the exponent as a second positional argument. +Integer powers are defined for any square matrix (invertible for negative powers) and reduce to repeated multiplication, while fractional powers are principal powers subject to the domain considerations above. + +```@docs; canonical=false +power +``` From 793a697e76ad624c397d62c8eca94340f6660307 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 16:39:30 -0400 Subject: [PATCH 04/22] refactor: reuse existing code patterns in matrix-function kernels Move `_mul_herm!` from `implementations/polar.jl` to `implementations/projections.jl` and use it for the symmetric-product kernels of `squareroot` and fractional `power` via `MatrixFunctionViaEigh`, so hermitian output is guaranteed through the same mechanism as the polar decomposition (BLAS `herk`/`syrk`, explicit projection otherwise, GPU overrides). Rewrite the `MatrixFunctionViaLA` kernels in the style of the `exponential` implementation: an explicit eltype branch with a `real.(...)` assignment instead of the bespoke `_copy_result!` helper. `squareroot`/`logarithm` use a strict complex-result check (LinearAlgebra returns real results whenever the principal value is real), while `power` tolerates rounding-level imaginary components since fractional powers are computed in complex arithmetic upstream. --- src/implementations/logarithm.jl | 9 +++++-- src/implementations/matrixfunctions.jl | 35 ++++++-------------------- src/implementations/polar.jl | 13 ---------- src/implementations/power.jl | 19 +++++++++++--- src/implementations/projections.jl | 13 ++++++++++ src/implementations/squareroot.jl | 12 ++++++--- 6 files changed, 51 insertions(+), 50 deletions(-) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index 0625323c3..fcb2e857f 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -35,8 +35,13 @@ initialize_output(::typeof(logarithm!), A::AbstractMatrix, ::AbstractAlgorithm) function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaLA) check_input(logarithm!, A, logA, alg) isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `logarithm`")) - result = LinearAlgebra.log(A) - _copy_result!(logarithm!, logA, result) + # `LinearAlgebra.log` of a real matrix is real whenever the principal logarithm is, + # so a complex result with a real output signals a genuine domain violation + logAc = LinearAlgebra.log(A) + if eltype(logAc) <: Complex && !(eltype(logA) <: Complex) + throw(_realness_domainerror(logarithm!)) + end + copy!(logA, logAc) return logA end diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index b6ba5a8e8..df524dceb 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -61,31 +61,12 @@ function _check_nonzero_eigenvalues(λ, atol::Real) end # For `MatrixFunctionViaLA`, domain violations surface as a complex result from -# `LinearAlgebra`; detect this after the fact to preserve type-stable output. -# `LinearAlgebra` sometimes computes a real result in complex arithmetic (e.g. fractional -# powers via a complex Schur decomposition), so rounding-level imaginary components do -# not signal a domain violation and are simply discarded. -function _copy_result!(f, out::AbstractMatrix, result::AbstractMatrix) - if eltype(result) <: Complex && !(eltype(out) <: Complex) - # base the tolerance on the working precision, which is the lower of the two: - # `LinearAlgebra` may internally promote, e.g. `Float32` fractional powers are - # computed with `Float64`-eltype results but `Float32`-level accuracy - atol = max(defaulttol(out), defaulttol(result)) * norm(result, Inf) - for x in result - if abs(imag(x)) > atol - throw( - DomainError( - f, - "The result of this matrix function applied to the given real matrix is complex (eigenvalues on the negative real axis). " * - "Pass a complex matrix to obtain the principal value, or use `MatrixFunctionViaEigh`/`MatrixFunctionViaEig` with a suitable " * - "`domain_atol` if the offending eigenvalues are rounding artifacts." - ) - ) - end - end - out .= real.(result) - else - copy!(out, result) - end - return out +# `LinearAlgebra` while the output should remain real. +function _realness_domainerror(f) + return DomainError( + f, + "The result of this matrix function applied to the given real matrix is complex (eigenvalues on the negative real axis). " * + "Pass a complex matrix to obtain the principal value, or use `MatrixFunctionViaEigh`/`MatrixFunctionViaEig` with a suitable " * + "`domain_atol` if the offending eigenvalues are rounding artifacts." + ) end diff --git a/src/implementations/polar.jl b/src/implementations/polar.jl index eb9b21d43..04738c5e1 100644 --- a/src/implementations/polar.jl +++ b/src/implementations/polar.jl @@ -81,19 +81,6 @@ function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarViaSVD) return (P, Wᴴ) end -# Implement `mul!(C, A', A)` and guarantee the result is hermitian. -# For BLAS calls that dispatch to `syrk` or `herk` this works automatically -# for GPU this currently does not seem to be guaranteed so we manually project -function _mul_herm!(C, A) - mul!(C, A, A') - project_hermitian!(C) - return C -end -function _mul_herm!(C::YALAPACK.BlasMat{T}, A::YALAPACK.BlasMat{T}) where {T <: YALAPACK.BlasFloat} - mul!(C, A, A') - return C -end - # Implementation via Newton # -------------------------- function left_polar!(A::AbstractMatrix, WP, alg::PolarNewton) diff --git a/src/implementations/power.jl b/src/implementations/power.jl index c73a18dd3..880885d1f 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -35,8 +35,19 @@ initialize_output(::typeof(power!), A::AbstractMatrix, p::Real, ::AbstractAlgori function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaLA) check_input(power!, A, p, powA, alg) isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `power`")) - result = A^p - _copy_result!(power!, powA, result) + powAc = A^p + if eltype(powAc) <: Complex && !(eltype(powA) <: Complex) + # `LinearAlgebra` computes fractional powers of real matrices in complex + # arithmetic and only casts back to real when the result is exactly real, + # so rounding-level imaginary components do not signal a domain violation. + # The tolerance is based on the working precision, which may be lower than + # the result eltype suggests (e.g. `Float32` input promotes to `ComplexF64`). + atol = defaulttol(powA) * norm(powAc, Inf) + all(x -> abs(imag(x)) <= atol, powAc) || throw(_realness_domainerror(power!)) + powA .= real.(powAc) + return powA + end + copy!(powA, powAc) return powA end @@ -49,6 +60,7 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) λ .= λ .^ p VD = V * D mul!(powA, VD, V') + return project_hermitian!(powA) else atol = something(alg.domain_atol, default_domain_atol(λ)) p < 0 && _check_nonzero_eigenvalues(λ, atol) @@ -56,9 +68,8 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction λ .= λ .^ (p / 2) Vs = rmul!(V, D) - mul!(powA, Vs, Vs') + return _mul_herm!(powA, Vs) end - return project_hermitian!(powA) end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) diff --git a/src/implementations/projections.jl b/src/implementations/projections.jl index c7f1f5006..b3c2aa567 100644 --- a/src/implementations/projections.jl +++ b/src/implementations/projections.jl @@ -140,3 +140,16 @@ end _imimag(x::Real) = zero(x) _imimag(x::Complex) = im * imag(x) + +# Implement `mul!(C, A, A')` and guarantee the result is hermitian. +# For BLAS calls that dispatch to `syrk` or `herk` this works automatically +# for GPU this currently does not seem to be guaranteed so we manually project +function _mul_herm!(C, A) + mul!(C, A, A') + project_hermitian!(C) + return C +end +function _mul_herm!(C::YALAPACK.BlasMat{T}, A::YALAPACK.BlasMat{T}) where {T <: YALAPACK.BlasFloat} + mul!(C, A, A') + return C +end diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index fe687e2c4..289086cce 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -35,8 +35,13 @@ initialize_output(::typeof(squareroot!), A::AbstractMatrix, ::AbstractAlgorithm) function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaLA) check_input(squareroot!, A, sqrtA, alg) isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `squareroot`")) - result = LinearAlgebra.sqrt(A) - _copy_result!(squareroot!, sqrtA, result) + # `LinearAlgebra.sqrt` of a real matrix is real whenever the principal square root is, + # so a complex result with a real output signals a genuine domain violation + sqrtAc = LinearAlgebra.sqrt(A) + if eltype(sqrtAc) <: Complex && !(eltype(sqrtA) <: Complex) + throw(_realness_domainerror(squareroot!)) + end + copy!(sqrtA, sqrtAc) return sqrtA end @@ -49,8 +54,7 @@ function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEigh) # `sqrt(A) = (V * D^(1/4)) * (V * D^(1/4))'` is hermitian by construction λ .= sqrt.(sqrt.(λ)) Vs = rmul!(V, D) - mul!(sqrtA, Vs, Vs') - return project_hermitian!(sqrtA) + return _mul_herm!(sqrtA, Vs) end function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig) From a40b26c7bf19600f6f2b0acc2cb24caf1a9129ad Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 19:21:08 -0400 Subject: [PATCH 05/22] refactor: delegate eigenvalue transforms to the DiagonalAlgorithm kernels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In the `MatrixFunctionViaEig(h)` implementations of `squareroot`, `logarithm` and `power`, apply the scalar function to the eigenvalue matrix by calling the function's own `DiagonalAlgorithm` kernel (which already contains the domain checks and tolerance clamping), mirroring how `exponential!` recurses into `exponential!((τ, D), D)`. The eig kernels only retain the negative-real-axis scan of the complex eigenvalues for real input, which has no `Diagonal` counterpart. --- src/implementations/logarithm.jl | 22 +++++++-------------- src/implementations/power.jl | 33 +++++++++++-------------------- src/implementations/squareroot.jl | 22 +++++++++------------ 3 files changed, 27 insertions(+), 50 deletions(-) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index fcb2e857f..b5c2daa7e 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -48,12 +48,7 @@ end function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh) check_input(logarithm!, A, logA, alg) D, V = eigh_full!(A, alg.eigh_alg) - λ = diagview(D) - atol = something(alg.domain_atol, default_domain_atol(λ)) - _check_nonzero_eigenvalues(λ, atol) - _clamp_domain_eigenvalues!(λ, atol) - λ .= log.(λ) - VD = V * D + VD = V * logarithm!(D, D, DiagonalAlgorithm(; domain_atol = alg.domain_atol)) mul!(logA, VD, V') return project_hermitian!(logA) end @@ -61,18 +56,15 @@ end function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig) check_input(logarithm!, A, logA, alg) D, V = eig_full!(A, alg.eig_alg) - λ = diagview(D) - atol = something(alg.domain_atol, default_domain_atol(λ)) - _check_nonzero_eigenvalues(λ, atol) + diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real - _clamp_domain_eigenvalues!(λ, atol) - λ .= log.(λ) - VD = V * D - logAc = rdiv!(VD, LinearAlgebra.lu!(V)) + atol = something(alg.domain_atol, default_domain_atol(diagview(D))) + _clamp_domain_eigenvalues!(diagview(D), atol) + VlogD = V * logarithm!(D, D, diag_alg) + logAc = rdiv!(VlogD, LinearAlgebra.lu!(V)) return logA .= real.(logAc) else - λ .= log.(λ) - logA .= V .* transpose(λ) + logA .= V .* transpose(diagview(logarithm!(D, D, diag_alg))) return rdiv!(logA, LinearAlgebra.lu!(V)) end end diff --git a/src/implementations/power.jl b/src/implementations/power.jl index 880885d1f..a85701ec7 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -54,20 +54,14 @@ end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) check_input(power!, A, p, powA, alg) D, V = eigh_full!(A, alg.eigh_alg) - λ = diagview(D) + diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if isinteger(p) - p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) - λ .= λ .^ p - VD = V * D + VD = V * power!(D, p, D, diag_alg) mul!(powA, VD, V') return project_hermitian!(powA) else - atol = something(alg.domain_atol, default_domain_atol(λ)) - p < 0 && _check_nonzero_eigenvalues(λ, atol) - _clamp_domain_eigenvalues!(λ, atol) # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction - λ .= λ .^ (p / 2) - Vs = rmul!(V, D) + Vs = rmul!(V, power!(D, p / 2, D, diag_alg)) return _mul_herm!(powA, Vs) end end @@ -75,22 +69,17 @@ end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) check_input(power!, A, p, powA, alg) D, V = eig_full!(A, alg.eig_alg) - λ = diagview(D) - if isinteger(p) - p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) - else - atol = something(alg.domain_atol, default_domain_atol(λ)) - p < 0 && _check_nonzero_eigenvalues(λ, atol) - eltype(A) <: Real && _clamp_domain_eigenvalues!(λ, atol) - end + diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real - λ .= λ .^ p - VD = V * D - powAc = rdiv!(VD, LinearAlgebra.lu!(V)) + if !isinteger(p) + atol = something(alg.domain_atol, default_domain_atol(diagview(D))) + _clamp_domain_eigenvalues!(diagview(D), atol) + end + VpD = V * power!(D, p, D, diag_alg) + powAc = rdiv!(VpD, LinearAlgebra.lu!(V)) return powA .= real.(powAc) else - λ .= λ .^ p - powA .= V .* transpose(λ) + powA .= V .* transpose(diagview(power!(D, p, D, diag_alg))) return rdiv!(powA, LinearAlgebra.lu!(V)) end end diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index 289086cce..1d8e9a136 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -48,29 +48,25 @@ end function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEigh) check_input(squareroot!, A, sqrtA, alg) D, V = eigh_full!(A, alg.eigh_alg) - λ = diagview(D) - atol = something(alg.domain_atol, default_domain_atol(λ)) - _clamp_domain_eigenvalues!(λ, atol) + diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) # `sqrt(A) = (V * D^(1/4)) * (V * D^(1/4))'` is hermitian by construction - λ .= sqrt.(sqrt.(λ)) - Vs = rmul!(V, D) + sqrtD = squareroot!(D, D, diag_alg) + Vs = rmul!(V, squareroot!(sqrtD, sqrtD, diag_alg)) return _mul_herm!(sqrtA, Vs) end function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig) check_input(squareroot!, A, sqrtA, alg) D, V = eig_full!(A, alg.eig_alg) - λ = diagview(D) + diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real - atol = something(alg.domain_atol, default_domain_atol(λ)) - _clamp_domain_eigenvalues!(λ, atol) - λ .= sqrt.(λ) - VD = V * D - sqrtAc = rdiv!(VD, LinearAlgebra.lu!(V)) + atol = something(alg.domain_atol, default_domain_atol(diagview(D))) + _clamp_domain_eigenvalues!(diagview(D), atol) + VsqrtD = V * squareroot!(D, D, diag_alg) + sqrtAc = rdiv!(VsqrtD, LinearAlgebra.lu!(V)) return sqrtA .= real.(sqrtAc) else - λ .= sqrt.(λ) - sqrtA .= V .* transpose(λ) + sqrtA .= V .* transpose(diagview(squareroot!(D, D, diag_alg))) return rdiv!(sqrtA, LinearAlgebra.lu!(V)) end end From 922948e44bea2b53ba3a5e5dd162d6b0037c7eb9 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 19:36:20 -0400 Subject: [PATCH 06/22] refactor: fold the eigenvalue domain scan into a single helper call Add a `_clamp_domain_eigenvalues!(D::Diagonal, domain_atol)` convenience method that derives the default tolerance itself, so the real-input branch of the `MatrixFunctionViaEig` kernels reduces to a one-line preamble before delegating to the `DiagonalAlgorithm` kernel. --- src/implementations/logarithm.jl | 3 +-- src/implementations/matrixfunctions.jl | 8 ++++++++ src/implementations/power.jl | 5 +---- src/implementations/squareroot.jl | 3 +-- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index b5c2daa7e..f86b11d58 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -58,8 +58,7 @@ function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig) D, V = eig_full!(A, alg.eig_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real - atol = something(alg.domain_atol, default_domain_atol(diagview(D))) - _clamp_domain_eigenvalues!(diagview(D), atol) + _clamp_domain_eigenvalues!(D, alg.domain_atol) VlogD = V * logarithm!(D, D, diag_alg) logAc = rdiv!(VlogD, LinearAlgebra.lu!(V)) return logA .= real.(logAc) diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index df524dceb..ffe15fa2b 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -22,6 +22,14 @@ function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) return λ end +# Convenience method for the eigenvalues of a decomposition, deriving the default +# tolerance from the eigenvalues themselves when `domain_atol` is `nothing`. +function _clamp_domain_eigenvalues!(D::Diagonal, domain_atol::Union{Nothing, Real}) + λ = diagview(D) + atol = something(domain_atol, default_domain_atol(λ)) + return _clamp_domain_eigenvalues!(λ, atol) +end + # Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative # real axis obstruct a real result; complex-conjugate pairs do not. function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real) diff --git a/src/implementations/power.jl b/src/implementations/power.jl index a85701ec7..c0e0b7cb5 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -71,10 +71,7 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) D, V = eig_full!(A, alg.eig_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real - if !isinteger(p) - atol = something(alg.domain_atol, default_domain_atol(diagview(D))) - _clamp_domain_eigenvalues!(diagview(D), atol) - end + isinteger(p) || _clamp_domain_eigenvalues!(D, alg.domain_atol) VpD = V * power!(D, p, D, diag_alg) powAc = rdiv!(VpD, LinearAlgebra.lu!(V)) return powA .= real.(powAc) diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index 1d8e9a136..aa5661b81 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -60,8 +60,7 @@ function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig) D, V = eig_full!(A, alg.eig_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real - atol = something(alg.domain_atol, default_domain_atol(diagview(D))) - _clamp_domain_eigenvalues!(diagview(D), atol) + _clamp_domain_eigenvalues!(D, alg.domain_atol) VsqrtD = V * squareroot!(D, D, diag_alg) sqrtAc = rdiv!(VsqrtD, LinearAlgebra.lu!(V)) return sqrtA .= real.(sqrtAc) From 22e6b53e974ab92d773af63f6fa05b9cfea4eb39 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 8 Jul 2026 19:57:05 -0400 Subject: [PATCH 07/22] refactor: make domain-clamp helpers GPU compatible Replace the scalar-indexing loops in `_clamp_domain_eigenvalues!` and `_check_nonzero_eigenvalues` with reductions and broadcasts, so the helpers work on GPU arrays. The `DomainError` now reports the worst offending eigenvalue instead of the first one encountered. --- src/implementations/matrixfunctions.jl | 62 +++++++++++--------------- 1 file changed, 27 insertions(+), 35 deletions(-) diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index ffe15fa2b..298da7261 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -5,20 +5,17 @@ # and throw a `DomainError` for eigenvalues that are genuinely negative, since then the # result cannot be expressed with the same (real) scalar type. function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) - for i in eachindex(λ) - x = λ[i] - if x < -atol - throw( - DomainError( - x, - "The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " * - "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." - ) + λmin = minimum(λ; init = zero(eltype(λ))) + if λmin < -atol + throw( + DomainError( + λmin, + "The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " * + "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." ) - elseif x < 0 - λ[i] = zero(x) - end + ) end + λ .= max.(λ, zero(eltype(λ))) return λ end @@ -33,37 +30,32 @@ end # Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative # real axis obstruct a real result; complex-conjugate pairs do not. function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real) - for i in eachindex(λ) - x = λ[i] - if abs(imag(x)) <= atol && real(x) < 0 - if real(x) < -atol - throw( - DomainError( - x, - "The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " * - "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." - ) - ) - else - λ[i] = zero(x) - end - end + onaxis = x -> abs(imag(x)) <= atol && real(x) < 0 + λmin = mapreduce(x -> onaxis(x) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ)))) + if λmin < -atol + throw( + DomainError( + λmin, + "The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " * + "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." + ) + ) end + λ .= ifelse.(onaxis.(λ), zero(eltype(λ)), λ) return λ end # Reject (numerically) zero eigenvalues for functions that are undefined there, # e.g. `logarithm` and `power` with a negative fractional power. function _check_nonzero_eigenvalues(λ, atol::Real) - for x in λ - if abs(x) <= atol - throw( - DomainError( - x, - "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined." - ) + amin = minimum(abs, λ; init = typemax(real(eltype(λ)))) + if amin <= atol + throw( + DomainError( + amin, + "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined." ) - end + ) end return λ end From 759fe50f3ef168eeb3c88092a9d8ce79fa0b06ee Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 27 Jul 2026 17:16:30 -0400 Subject: [PATCH 08/22] fix: correct the out-of-domain `DomainError` messages The `lazy"..."` literal introduced for the negative-eigenvalue error does not process the `\` line continuation, so the rendered message contained a literal backslash, newline and indentation. Use ordinary string concatenation, as elsewhere in this file, and report the tolerance under its user-facing name `domain_atol` rather than `atol`. Restore the distinction between the real ("a negative real eigenvalue") and complex ("an eigenvalue on the negative real axis") variants by passing the description to the shared helper, and give `_check_nonzero_eigenvalues` the same `@noinline` throw helper so both domain checks stay free of error-path code on the GPU. Thread `domain_atol` back into the diagonal kernels invoked by the `squareroot` eigendecomposition paths, and let the `DiagonalAlgorithm` kernel reuse `_clamp_domain_eigenvalues!` again so it matches `power!`. Co-Authored-By: Claude Opus 5 (1M context) --- src/implementations/matrixfunctions.jl | 50 +++++++++++++------------- src/implementations/squareroot.jl | 6 ++-- src/interface/logarithm.jl | 12 +++---- 3 files changed, 30 insertions(+), 38 deletions(-) diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index 298da7261..4e72dacf7 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -1,20 +1,33 @@ # Shared helpers for matrix functions with a restricted domain # ------------------------------------------------------------- +# The throwing branches live in `@noinline` helpers so that the reductions and broadcasts +# below stay free of error-path code, which keeps them GPU friendly. +@noinline function throw_negative_eigenvalue(λmin, atol, what) + return throw( + DomainError( + λmin, + "The matrix has $what beyond `domain_atol = $atol` and the result of this matrix function is complex. " * + "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." + ) + ) +end + +@noinline function throw_zero_eigenvalue(amin, atol) + return throw( + DomainError( + amin, + "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined." + ) + ) +end + # Clamp real eigenvalues that are negative within `atol` (rounding artifacts) to zero, # and throw a `DomainError` for eigenvalues that are genuinely negative, since then the # result cannot be expressed with the same (real) scalar type. function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) λmin = minimum(λ; init = zero(eltype(λ))) - if λmin < -atol - throw( - DomainError( - λmin, - "The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " * - "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." - ) - ) - end + λmin < -atol && throw_negative_eigenvalue(λmin, atol, "a negative real eigenvalue") λ .= max.(λ, zero(eltype(λ))) return λ end @@ -32,15 +45,7 @@ end function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real) onaxis = x -> abs(imag(x)) <= atol && real(x) < 0 λmin = mapreduce(x -> onaxis(x) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ)))) - if λmin < -atol - throw( - DomainError( - λmin, - "The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " * - "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." - ) - ) - end + λmin < -atol && throw_negative_eigenvalue(λmin, atol, "an eigenvalue on the negative real axis") λ .= ifelse.(onaxis.(λ), zero(eltype(λ)), λ) return λ end @@ -49,14 +54,7 @@ end # e.g. `logarithm` and `power` with a negative fractional power. function _check_nonzero_eigenvalues(λ, atol::Real) amin = minimum(abs, λ; init = typemax(real(eltype(λ)))) - if amin <= atol - throw( - DomainError( - amin, - "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined." - ) - ) - end + amin <= atol && throw_zero_eigenvalue(amin, atol) return λ end diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index aa5661b81..2a5221c89 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -38,9 +38,8 @@ function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaLA) # `LinearAlgebra.sqrt` of a real matrix is real whenever the principal square root is, # so a complex result with a real output signals a genuine domain violation sqrtAc = LinearAlgebra.sqrt(A) - if eltype(sqrtAc) <: Complex && !(eltype(sqrtA) <: Complex) + eltype(sqrtAc) <: Complex && !(eltype(sqrtA) <: Complex) && throw(_realness_domainerror(squareroot!)) - end copy!(sqrtA, sqrtAc) return sqrtA end @@ -50,8 +49,7 @@ function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEigh) D, V = eigh_full!(A, alg.eigh_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) # `sqrt(A) = (V * D^(1/4)) * (V * D^(1/4))'` is hermitian by construction - sqrtD = squareroot!(D, D, diag_alg) - Vs = rmul!(V, squareroot!(sqrtD, sqrtD, diag_alg)) + Vs = rmul!(V, power!(D, 1 // 4, D, diag_alg)) return _mul_herm!(sqrtA, Vs) end diff --git a/src/interface/logarithm.jl b/src/interface/logarithm.jl index 39acb133c..b84499c14 100644 --- a/src/interface/logarithm.jl +++ b/src/interface/logarithm.jl @@ -11,18 +11,14 @@ Compute the principal logarithm `logA` of the square matrix `A`, i.e. the logari whose eigenvalues have imaginary part in `(-π, π]`. The scalar type of the output matches that of the input. -As a consequence, a real matrix with eigenvalues on the negative real axis, for which -the principal logarithm is complex, leads to a `DomainError`; pass a complex matrix -to obtain the principal value. -A matrix with (numerically) zero eigenvalues has no logarithm and also leads to a -`DomainError`. +As a consequence, a real matrix with eigenvalues on the negative real axis, for which the principal logarithm is complex, leads to a `DomainError`; pass a complex matrix to obtain the principal value. +A matrix with (numerically) zero eigenvalues has no logarithm and also leads to a `DomainError`. Both checks use a tolerance `domain_atol`, which can be specified for the algorithms that support it and defaults to [`default_domain_atol`](@ref). !!! note - The bang method `logarithm!` optionally accepts the output structure and - possibly destroys the input matrix `A`. Always use the return value of the function - as it may not always be possible to use the provided `logA` as output. + The bang method `logarithm!` optionally accepts the output structure and possibly destroys the input matrix `A`. + Always use the return value of the function as it may not always be possible to use the provided `logA` as output. """ @functiondef logarithm From f2110d441f38f1b556436eccc2736249f847f1b0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 27 Jul 2026 19:05:40 -0400 Subject: [PATCH 09/22] test: cover `squareroot`, `logarithm` and `power` on GPU The domain helpers were rewritten as reductions and broadcasts so they would run on GPU arrays, but nothing exercised that. Add `test_{squareroot,logarithm,power}_gpu`, modelled on the existing `test_exponential_gpu`, comparing device results against the CPU reference for the `DiagonalAlgorithm` fast path and the `MatrixFunctionViaEigh` path. Both domain outcomes are covered, since the clamp and the throw are precisely what the scalar loops were replaced by: roundoff-scale negative eigenvalues must clamp, genuinely negative ones must throw a `DomainError`, and for `logarithm` and negative fractional `power` a zero eigenvalue must throw as well. Scalar indexing anywhere along these paths errors under the GPUArrays guard, so a regression to the previous loops would fail these tests. Co-Authored-By: Claude Opus 5 (1M context) --- test/logarithm.jl | 43 ++++++++++++++++++++++++++++++++++++++ test/power.jl | 51 ++++++++++++++++++++++++++++++++++++++++++++++ test/squareroot.jl | 48 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 142 insertions(+) diff --git a/test/logarithm.jl b/test/logarithm.jl index 2307141a5..5a1e7f6d5 100644 --- a/test/logarithm.jl +++ b/test/logarithm.jl @@ -4,6 +4,7 @@ using TestExtras using StableRNGs using MatrixAlgebraKit: diagview using LinearAlgebra +using CUDA, AMDGPU using LinearAlgebra: exp BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) @@ -100,3 +101,45 @@ end end @test_throws DomainError logarithm(Diagonal(zeros(T, m))) end + +# GPU tests +# --------- +# As for `squareroot`, the `DiagonalAlgorithm` kernel and the shared domain helpers are +# backend-generic, and `MatrixFunctionViaEigh` goes through the device `eigh_full!`. Both +# domain rejections (negative real axis, and the zero eigenvalue where no logarithm exists) +# are exercised, since `_check_nonzero_eigenvalues` is a reduction too. If any step fell back +# to scalar indexing these would error under GPUArrays' scalar-indexing guard. +function test_logarithm_gpu(ArrayT, T) + rng = StableRNG(123) + m = 54 + + # Diagonal fast path + data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) + @test Array(diagview(logarithm(Diagonal(ArrayT(data))))) ≈ + diagview(logarithm(Diagonal(data))) + + # hermitian positive definite, via the eigenvalue decomposition + X = randn(rng, T, m, m) + A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) + A_gpu = ArrayT(A) + alg_gpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A_gpu)) + alg_cpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A)) + @test Array(logarithm(A_gpu, alg_gpu)) ≈ logarithm(A, alg_cpu) + + # domain handling: negative real axis, and no logarithm at a zero eigenvalue + T <: Real && @test_throws DomainError logarithm(Diagonal(ArrayT(-data))) + @test_throws DomainError logarithm(Diagonal(ArrayT(zeros(T, m)))) + return nothing +end + +if CUDA.functional() + @testset "logarithm on CUDA for T = $T" for T in BLASFloats + test_logarithm_gpu(CuArray, T) + end +end + +if AMDGPU.functional() + @testset "logarithm on AMDGPU for T = $T" for T in BLASFloats + test_logarithm_gpu(ROCArray, T) + end +end diff --git a/test/power.jl b/test/power.jl index 99cf659b1..6322c784d 100644 --- a/test/power.jl +++ b/test/power.jl @@ -4,6 +4,7 @@ using TestExtras using StableRNGs using MatrixAlgebraKit: diagview using LinearAlgebra +using CUDA, AMDGPU BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) @@ -124,3 +125,53 @@ end @test_throws SingularException power(Diagonal(zeros(T, m)), -1) @test_throws DomainError power(Diagonal(zeros(T, m)), -0.5) end + +# GPU tests +# --------- +# As for `squareroot`, the `DiagonalAlgorithm` kernel and the shared domain helpers are +# backend-generic, and `MatrixFunctionViaEigh` goes through the device `eigh_full!`. Both the +# integer branch (which only checks for singularity) and the fractional branch (which clamps +# and rejects out-of-domain eigenvalues) are exercised. If any step fell back to scalar +# indexing these would error under GPUArrays' scalar-indexing guard. +function test_power_gpu(ArrayT, T) + rng = StableRNG(123) + m = 54 + + # Diagonal fast path, integer and fractional exponents + data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) + for p in (2, -1, 0.5, -0.25) + @test Array(diagview(power(Diagonal(ArrayT(data)), p))) ≈ + diagview(power(Diagonal(data), p)) + end + + # hermitian positive definite, via the eigenvalue decomposition + X = randn(rng, T, m, m) + A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) + A_gpu = ArrayT(A) + alg_gpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A_gpu)) + alg_cpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A)) + for p in (2, -1, 0.5, -0.5) + @test Array(power(A_gpu, p, alg_gpu)) ≈ power(A, p, alg_cpu) + end + + # domain handling: integer powers ignore the domain, fractional ones do not + if T <: Real + @test Array(diagview(power(Diagonal(ArrayT(-data)), 2))) ≈ data .^ 2 + @test_throws DomainError power(Diagonal(ArrayT(-data)), 0.5) + end + @test_throws SingularException power(Diagonal(ArrayT(zeros(T, m))), -1) + @test_throws DomainError power(Diagonal(ArrayT(zeros(T, m))), -0.5) + return nothing +end + +if CUDA.functional() + @testset "power on CUDA for T = $T" for T in BLASFloats + test_power_gpu(CuArray, T) + end +end + +if AMDGPU.functional() + @testset "power on AMDGPU for T = $T" for T in BLASFloats + test_power_gpu(ROCArray, T) + end +end diff --git a/test/squareroot.jl b/test/squareroot.jl index 0d02a2f5a..478548606 100644 --- a/test/squareroot.jl +++ b/test/squareroot.jl @@ -4,6 +4,7 @@ using TestExtras using StableRNGs using MatrixAlgebraKit: diagview using LinearAlgebra +using CUDA, AMDGPU BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) @@ -104,3 +105,50 @@ end @test squareroot(Diagonal(T[-eps(T), 1])) ≈ Diagonal(T[0, 1]) end end + +# GPU tests +# --------- +# The `DiagonalAlgorithm` kernel and the shared domain helpers are written as reductions and +# broadcasts rather than scalar loops, so the same code runs on GPU; the `MatrixFunctionViaEigh` +# path additionally goes through the device `eigh_full!`. Compare against the CPU reference and +# exercise both domain outcomes, since the clamp and the throw are exactly what the reductions +# replaced. If any step fell back to scalar indexing these would error under GPUArrays' +# scalar-indexing guard. +function test_squareroot_gpu(ArrayT, T) + rng = StableRNG(123) + m = 54 + + # Diagonal fast path + data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : randn(rng, T, m) + @test Array(diagview(squareroot(Diagonal(ArrayT(data))))) ≈ + diagview(squareroot(Diagonal(data))) + + # hermitian positive definite, via the eigenvalue decomposition + X = randn(rng, T, m, m) + A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) + A_gpu = ArrayT(A) + alg_gpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A_gpu)) + alg_cpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A)) + sqrtA = squareroot(A_gpu, alg_gpu) + @test Array(sqrtA) ≈ squareroot(A, alg_cpu) + @test Array(sqrtA) * Array(sqrtA) ≈ A + + # domain handling: roundoff-scale negatives are clamped, genuine ones throw + if T <: Real + @test Array(diagview(squareroot(Diagonal(ArrayT(T[-eps(T), 1]))))) ≈ T[0, 1] + @test_throws DomainError squareroot(Diagonal(ArrayT(-data))) + end + return nothing +end + +if CUDA.functional() + @testset "squareroot on CUDA for T = $T" for T in BLASFloats + test_squareroot_gpu(CuArray, T) + end +end + +if AMDGPU.functional() + @testset "squareroot on AMDGPU for T = $T" for T in BLASFloats + test_squareroot_gpu(ROCArray, T) + end +end From 2db280b3231569a3a84d66ecef73d6125758e74a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 08:58:29 -0400 Subject: [PATCH 10/22] refactor: build `Diagonal` copies via `map_diagonal` `map_diagonal(f, src)` is exactly `diagonal(f.(diagview(src)))`, so the three `copy_input` methods can use it instead of constructing a `Diagonal` directly. This keeps them consistent with the rest of the package, which routes diagonal construction through `diagonal` rather than the `LinearAlgebra` type. Co-Authored-By: Claude Opus 5 (1M context) --- src/implementations/logarithm.jl | 2 +- src/implementations/power.jl | 2 +- src/implementations/squareroot.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index f86b11d58..e59140362 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -3,7 +3,7 @@ function copy_input(::typeof(logarithm), A::AbstractMatrix) return copy!(similar(A, float(eltype(A))), A) end -copy_input(::typeof(logarithm), A::Diagonal) = Diagonal(float.(diagview(A))) +copy_input(::typeof(logarithm), A::Diagonal) = map_diagonal(float, A) function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, alg::AbstractAlgorithm) m = LinearAlgebra.checksquare(A) diff --git a/src/implementations/power.jl b/src/implementations/power.jl index c0e0b7cb5..fff0ab7be 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -3,7 +3,7 @@ function copy_input(::typeof(power), A::AbstractMatrix, p::Real) return copy!(similar(A, float(eltype(A))), A), p end -copy_input(::typeof(power), A::Diagonal, p::Real) = Diagonal(float.(diagview(A))), p +copy_input(::typeof(power), A::Diagonal, p::Real) = map_diagonal(float, A), p function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, alg::AbstractAlgorithm) m = LinearAlgebra.checksquare(A) diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index 2a5221c89..67b08c3b3 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -3,7 +3,7 @@ function copy_input(::typeof(squareroot), A::AbstractMatrix) return copy!(similar(A, float(eltype(A))), A) end -copy_input(::typeof(squareroot), A::Diagonal) = Diagonal(float.(diagview(A))) +copy_input(::typeof(squareroot), A::Diagonal) = map_diagonal(float, A) function check_input(::typeof(squareroot!), A::AbstractMatrix, sqrtA, alg::AbstractAlgorithm) m = LinearAlgebra.checksquare(A) From 647c8dd6f1fed8f92c137e86ee0a03b7b1b197ff Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 11:29:54 -0400 Subject: [PATCH 11/22] test: move the matrix-function tests into the TestSuite MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The four matrix functions were tested by hand-written, CPU-only top-level files plus per-backend duplicates under `test/genericschur/` and `test/genericlinearalgebra/`, while every decomposition lives in `test/testsuite/` as array-type-parameterized functions driven by a thin `test//.jl`. Close that gap: `test/testsuite/matrixfunctions/` now holds the bodies and `test/matrixfunctions/` the drivers. Assertions become invariants rather than comparisons against `LinearAlgebra`, so one body serves host, CUDA, ROCm and downstream array types: `sqrtA * sqrtA ≈ A`, `exponential(logA) ≈ A`, `exp(A) * exp(-A) ≈ I`, and elementwise spectrum mapping via `eigh_vals`. Reference cross-checks against `LinearAlgebra` are kept, but confined to `test_*_reference`, which the drivers only call for host arrays. GPU coverage now falls out of the array-type parameterization, so the bespoke `test_*_gpu` functions are gone. Three new instantiators (`instantiate_offaxis_matrix`, `instantiate_posdef_matrix`, `instantiate_hermitian_spectrum`) give the controlled spectra these functions need, reusing `instantiate_unitary` and the `Diagonal` handling already in the suite. The domain rules differ per function and are now encoded rather than approximated: a roundoff-negative eigenvalue clamps for `squareroot` and positive fractional `power`, but still throws for `logarithm`, where clamping onto zero leaves the matrix singular. `MatrixFunctionViaLA` inspects only the realness of its result, so it accepts singular input that the spectrum-based algorithms reject; the `eigh`-based algorithms reject a negative eigenvalue whatever the scalar type, since the root of a hermitian matrix with a negative eigenvalue is not hermitian. Also drops the deprecated `LAPACK_Simple`/`LAPACK_QRIteration`/ `GS_QRIteration`/`GLA_QRIteration` spellings, which accounted for the deprecation warnings these tests emitted. Generic eltypes name their driver explicitly (`QRIteration(; driver = GS())`), since with both generic packages loaded `default_driver` resolves `QRIteration` to GLA, which provides no `geev!`. Full suite: 93826 passing, up from 90917. Co-Authored-By: Claude Opus 5 (1M context) --- test/common/exponential.jl | 145 -------------- test/genericlinearalgebra/exponential.jl | 46 ----- test/genericlinearalgebra/logarithm.jl | 28 --- test/genericlinearalgebra/power.jl | 29 --- test/genericlinearalgebra/squareroot.jl | 30 --- test/genericschur/exponential.jl | 47 ----- test/genericschur/logarithm.jl | 32 --- test/genericschur/power.jl | 31 --- test/genericschur/squareroot.jl | 32 --- test/logarithm.jl | 145 -------------- test/matrixfunctions/exponential.jl | 114 +++++++++++ test/matrixfunctions/logarithm.jl | 113 +++++++++++ test/matrixfunctions/power.jl | 113 +++++++++++ test/matrixfunctions/squareroot.jl | 112 +++++++++++ test/power.jl | 177 ----------------- test/squareroot.jl | 154 --------------- test/testsuite/TestSuite.jl | 41 ++++ test/testsuite/matrixfunctions/exponential.jl | 172 ++++++++++++++++ test/testsuite/matrixfunctions/logarithm.jl | 121 ++++++++++++ test/testsuite/matrixfunctions/power.jl | 184 ++++++++++++++++++ test/testsuite/matrixfunctions/squareroot.jl | 119 +++++++++++ 21 files changed, 1089 insertions(+), 896 deletions(-) delete mode 100644 test/common/exponential.jl delete mode 100644 test/genericlinearalgebra/exponential.jl delete mode 100644 test/genericlinearalgebra/logarithm.jl delete mode 100644 test/genericlinearalgebra/power.jl delete mode 100644 test/genericlinearalgebra/squareroot.jl delete mode 100644 test/genericschur/exponential.jl delete mode 100644 test/genericschur/logarithm.jl delete mode 100644 test/genericschur/power.jl delete mode 100644 test/genericschur/squareroot.jl delete mode 100644 test/logarithm.jl create mode 100644 test/matrixfunctions/exponential.jl create mode 100644 test/matrixfunctions/logarithm.jl create mode 100644 test/matrixfunctions/power.jl create mode 100644 test/matrixfunctions/squareroot.jl delete mode 100644 test/power.jl delete mode 100644 test/squareroot.jl create mode 100644 test/testsuite/matrixfunctions/exponential.jl create mode 100644 test/testsuite/matrixfunctions/logarithm.jl create mode 100644 test/testsuite/matrixfunctions/power.jl create mode 100644 test/testsuite/matrixfunctions/squareroot.jl diff --git a/test/common/exponential.jl b/test/common/exponential.jl deleted file mode 100644 index 981a0a265..000000000 --- a/test/common/exponential.jl +++ /dev/null @@ -1,145 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using LinearAlgebra: exp -using CUDA, AMDGPU - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) -GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) - -@testset "exponential! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - A = LinearAlgebra.normalize!(randn(rng, T, m, m)) - Ac = copy(A) - expA = LinearAlgebra.exp(A) - - expA2 = @constinferred exponential(A) - @test expA ≈ expA2 - @test A == Ac - - algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple()), MatrixFunctionViaTaylor()) - @testset "algorithm $alg" for alg in algs - expA2 = @constinferred exponential(A, alg) - @test expA ≈ expA2 - @test A == Ac - end - - @test_throws DomainError exponential(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) -end - -@testset "exponential! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - A = randn(rng, T, m, m) - τ = randn(rng, T) - Ac = copy(A) - - Aτ = A * τ - expAτ = LinearAlgebra.exp(Aτ) - - expAτ2 = @constinferred exponential((τ, A)) - @test expAτ ≈ expAτ2 - @test A == Ac - - algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple()), MatrixFunctionViaTaylor()) - @testset "algorithm $alg" for alg in algs - expAτ2 = @constinferred exponential((τ, A), alg) - @test expAτ ≈ expAτ2 - @test A == Ac - end - - @test_throws DomainError exponential((τ, A); alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) -end - -@testset "exponential! for non-Matrix input $T" for T in BLASFloats - rng = StableRNG(123) - m = 12 - A = LinearAlgebra.normalize!(randn(rng, T, m, m)) - expA = LinearAlgebra.exp(A) - - wrappers = ( - ("view", B -> view(B, :, :)), - ("PermutedDimsArray", B -> PermutedDimsArray(permutedims(B), (2, 1))), - ("ReshapedArray", B -> reshape(view(vec(B), 1:(m * m)), m, m)), - ) - @testset "$name" for (name, wrap) in wrappers - W = wrap(copy(A)) - @test !(W isa Matrix) - @test exponential!(W) ≈ expA - end -end - -@testset "exponential! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - m = 54 - - A = Diagonal(randn(rng, T, m)) - τ = randn(rng, T) - Ac = copy(A) - - expA = LinearAlgebra.exp(A) - - expA2 = @constinferred exponential(A) - @test expA ≈ expA2 - @test A == Ac -end - -@testset "exponential! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - m = 1 - - A = Diagonal(randn(rng, T, m)) - τ = randn(rng, T) - Ac = copy(A) - - Aτ = A * τ - expAτ = LinearAlgebra.exp(Aτ) - - expAτ2 = @constinferred exponential((τ, A)) - @test expAτ ≈ expAτ2 - @test A == Ac -end - -# GPU tests -# --------- -# The Taylor exponential is backend-generic, so the same code runs on GPU. Compare device -# results against the CPU reference, exercising both `balance` settings (fix 1 & 3) and the -# scaled `(τ, A)` entrypoint. A badly-scaled matrix exercises the balancing path. If any step -# fell back to scalar indexing these would error under GPUArrays' scalar-indexing guard. -function test_exponential_gpu(ArrayT, T) - rng = StableRNG(123) - m = 54 - A = randn(rng, T, m, m) ./ (2 * m) - τ = randn(rng, T) - - # badly-scaled similarity transform Aᵢⱼ ← Aᵢⱼ sᵢ / sⱼ, to give balancing work to do - s = exp10.(range(-real(T)(3), real(T)(3), length = m)) - Abad = A .* s ./ transpose(s) - - for M in (A, Abad) - M_gpu = ArrayT(M) - for alg in (MatrixFunctionViaTaylor(), MatrixFunctionViaTaylor(; balance = false)) - @test Array(exponential(M_gpu, alg)) ≈ exponential(M, alg) - @test Array(exponential((τ, M_gpu), alg)) ≈ exponential((τ, M), alg) - end - end - return nothing -end - -if CUDA.functional() - @testset "exponential on CUDA for T = $T" for T in BLASFloats - test_exponential_gpu(CuArray, T) - end -end - -if AMDGPU.functional() - @testset "exponential on AMDGPU for T = $T" for T in BLASFloats - test_exponential_gpu(ROCArray, T) - end -end diff --git a/test/genericlinearalgebra/exponential.jl b/test/genericlinearalgebra/exponential.jl deleted file mode 100644 index 5743d3b4d..000000000 --- a/test/genericlinearalgebra/exponential.jl +++ /dev/null @@ -1,46 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericLinearAlgebra - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "exponential! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 54 - - A = project_hermitian!(randn(rng, T, m, m)) - D, V = @constinferred eigh_full(A) - algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) - @testset "algorithm $alg" for alg in algs - expA = @constinferred exponential!(copy(A); alg) - expA2 = @constinferred exponential(A; alg) - @test expA2 ≈ expA - - Dexp, Vexp = @constinferred eigh_full(expA) - @test diagview(Dexp) ≈ LinearAlgebra.exp.(diagview(D)) - end -end - -using GenericSchur -@testset "exponential! for T1 = $T1, T2 = $T2" for T1 in GenericFloats, T2 in GenericFloats - rng = StableRNG(123) - m = 54 - A = project_hermitian!(randn(rng, T1, m, m)) - τ = randn(rng, T2) - - D, V = @constinferred eigh_full(A) - algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) - @testset "algorithm $alg" for alg in algs - expτA = @constinferred exponential!((τ, copy(A)); alg) - expτA2 = @constinferred exponential((τ, A); alg) - @test expτA2 ≈ expτA - - Dexp, Vexp = @constinferred eig_full(expτA) - - @test sort(diagview(Dexp); by = real) ≈ sort(LinearAlgebra.exp.(diagview(D) .* τ); by = real) - end -end diff --git a/test/genericlinearalgebra/logarithm.jl b/test/genericlinearalgebra/logarithm.jl deleted file mode 100644 index 4cd69d342..000000000 --- a/test/genericlinearalgebra/logarithm.jl +++ /dev/null @@ -1,28 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericLinearAlgebra - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "logarithm! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 24 - - X = randn(rng, T, m, m) - A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I - D, V = @constinferred eigh_full(A) - - algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) - @testset "algorithm $alg" for alg in algs - logA = @constinferred logarithm!(copy(A); alg) - logA2 = @constinferred logarithm(A; alg) - @test logA2 ≈ logA - - Dlog, Vlog = @constinferred eigh_full(logA) - @test diagview(Dlog) ≈ log.(diagview(D)) - end -end diff --git a/test/genericlinearalgebra/power.jl b/test/genericlinearalgebra/power.jl deleted file mode 100644 index 1170be498..000000000 --- a/test/genericlinearalgebra/power.jl +++ /dev/null @@ -1,29 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericLinearAlgebra - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "power! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 24 - - X = randn(rng, T, m, m) - A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I - D, V = @constinferred eigh_full(A) - - algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) - @testset "algorithm $alg" for alg in algs - @test power(A, 3, alg) ≈ A * A * A - @test power(A, -1, alg) * A ≈ LinearAlgebra.I - @test power(A, big"0.5", alg) ≈ squareroot(A; alg) - - powA = @constinferred power(A, big"0.5", alg) - Dpow, Vpow = @constinferred eigh_full(powA) - @test diagview(Dpow) ≈ sqrt.(diagview(D)) - end -end diff --git a/test/genericlinearalgebra/squareroot.jl b/test/genericlinearalgebra/squareroot.jl deleted file mode 100644 index 1f9ee93af..000000000 --- a/test/genericlinearalgebra/squareroot.jl +++ /dev/null @@ -1,30 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericLinearAlgebra - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "squareroot! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 24 - - X = randn(rng, T, m, m) - A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I - D, V = @constinferred eigh_full(A) - - algs = (MatrixFunctionViaEigh(GLA_QRIteration()),) - @testset "algorithm $alg" for alg in algs - sqrtA = @constinferred squareroot!(copy(A); alg) - sqrtA2 = @constinferred squareroot(A; alg) - @test sqrtA2 ≈ sqrtA - @test sqrtA * sqrtA ≈ A - @test LinearAlgebra.ishermitian(sqrtA) - - Dsqrt, Vsqrt = @constinferred eigh_full(sqrtA) - @test diagview(Dsqrt) ≈ sqrt.(diagview(D)) - end -end diff --git a/test/genericschur/exponential.jl b/test/genericschur/exponential.jl deleted file mode 100644 index e2d26e2bd..000000000 --- a/test/genericschur/exponential.jl +++ /dev/null @@ -1,47 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericSchur - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "exponential! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 54 - - A = randn(rng, T, m, m) - D, V = @constinferred eig_full(A) - algs = (MatrixFunctionViaEig(GS_QRIteration()),) - expA_LA = @constinferred exponential(A) - @testset "algorithm $alg" for alg in algs - expA = @constinferred exponential!(copy(A)) - expA2 = @constinferred exponential(A; alg = alg) - @test expA ≈ expA_LA - @test expA2 ≈ expA - - Dexp, Vexp = @constinferred eig_full(expA) - @test sort(diagview(Dexp); by = imag) ≈ sort(LinearAlgebra.exp.(diagview(D)); by = imag) - end -end - -@testset "exponential! for T1 = $T1, T2 = $T2" for T1 in GenericFloats, T2 in GenericFloats - rng = StableRNG(123) - m = 54 - - A = randn(rng, T1, m, m) - τ = randn(rng, T2) - - D, V = @constinferred eig_full(A) - algs = (MatrixFunctionViaEig(GS_QRIteration()),) - @testset "algorithm $alg" for alg in algs - expτA = @constinferred exponential!((τ, copy(A))) - expτA2 = @constinferred exponential((τ, A); alg) - @test expτA2 ≈ expτA - - Dexp, Vexp = @constinferred eig_full(expτA) - @test sort(diagview(Dexp); by = x -> (imag(x), real(x))) ≈ sort(LinearAlgebra.exp.(diagview(D) .* τ); by = x -> (imag(x), real(x))) - end -end diff --git a/test/genericschur/logarithm.jl b/test/genericschur/logarithm.jl deleted file mode 100644 index aeaf2c7d4..000000000 --- a/test/genericschur/logarithm.jl +++ /dev/null @@ -1,32 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericSchur - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "logarithm! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 24 - - # spectrum inside a disk around 1, away from the negative real axis and zero - A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) - D, V = @constinferred eig_full(A) - - logA = @constinferred logarithm(A) - @test eltype(logA) == T - @test exponential(logA) ≈ A - - algs = (MatrixFunctionViaEig(GS_QRIteration()),) - @testset "algorithm $alg" for alg in algs - logA2 = @constinferred logarithm(A; alg) - @test logA2 ≈ logA - - Dlog, Vlog = @constinferred eig_full(logA2) - by = x -> (real(x), imag(x)) - @test sort(diagview(Dlog); by) ≈ sort(log.(diagview(D)); by) - end -end diff --git a/test/genericschur/power.jl b/test/genericschur/power.jl deleted file mode 100644 index f4b4a907b..000000000 --- a/test/genericschur/power.jl +++ /dev/null @@ -1,31 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericSchur - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "power! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 24 - - # spectrum inside a disk around 1, away from the negative real axis and zero - A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) - - powA = @constinferred power(A, 3) - @test powA ≈ A * A * A - @test eltype(powA) == T - - powA = @constinferred power(A, big"0.5") - @test powA ≈ squareroot(A) - - algs = (MatrixFunctionViaEig(GS_QRIteration()),) - @testset "algorithm $alg" for alg in algs - @test power(A, 3, alg) ≈ A * A * A - @test power(A, -1, alg) * A ≈ LinearAlgebra.I - @test power(A, big"0.5", alg) ≈ squareroot(A) - end -end diff --git a/test/genericschur/squareroot.jl b/test/genericschur/squareroot.jl deleted file mode 100644 index 33a2140e4..000000000 --- a/test/genericschur/squareroot.jl +++ /dev/null @@ -1,32 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using GenericSchur - -GenericFloats = (BigFloat, Complex{BigFloat}) - -@testset "squareroot! for T = $T" for T in GenericFloats - rng = StableRNG(123) - m = 24 - - # spectrum inside a disk around 1, away from the negative real axis - A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) - D, V = @constinferred eig_full(A) - - sqrtA = @constinferred squareroot(A) - @test sqrtA * sqrtA ≈ A - @test eltype(sqrtA) == T - - algs = (MatrixFunctionViaEig(GS_QRIteration()),) - @testset "algorithm $alg" for alg in algs - sqrtA2 = @constinferred squareroot(A; alg) - @test sqrtA2 ≈ sqrtA - - Dsqrt, Vsqrt = @constinferred eig_full(sqrtA2) - by = x -> (real(x), imag(x)) - @test sort(diagview(Dsqrt); by) ≈ sort(sqrt.(diagview(D)); by) - end -end diff --git a/test/logarithm.jl b/test/logarithm.jl deleted file mode 100644 index 5a1e7f6d5..000000000 --- a/test/logarithm.jl +++ /dev/null @@ -1,145 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using CUDA, AMDGPU -using LinearAlgebra: exp - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) -GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) - -@testset "logarithm! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - # spectrum inside a disk around 1, away from the negative real axis and zero - A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) - Ac = copy(A) - logA = LinearAlgebra.log(A) - - logA2 = @constinferred logarithm(A) - @test logA ≈ logA2 - @test exp(logA2) ≈ A - @test A == Ac - - algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) - @testset "algorithm $alg" for alg in algs - logA2 = @constinferred logarithm(A, alg) - @test logA ≈ logA2 - @test A == Ac - end - - @test_throws DomainError logarithm(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) - - # roundtrip with exponential - @test logarithm(exponential(logA)) ≈ logA -end - -@testset "logarithm! for hermitian T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - X = randn(rng, T, m, m) - A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) - Ac = copy(A) - logA = LinearAlgebra.log(LinearAlgebra.Hermitian(A)) - - algs = ( - MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()), - MatrixFunctionViaEig(LAPACK_Simple()), - ) - @testset "algorithm $alg" for alg in algs - logA2 = @constinferred logarithm(A, alg) - @test logA ≈ logA2 - @test A == Ac - end -end - -@testset "logarithm! domain handling for T = $T" for T in (Float32, Float64) - rng = StableRNG(123) - m = 4 - - X = randn(rng, T, m, m) - V = Matrix(LinearAlgebra.qr(X).Q) - A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V' - A = (A + A') / 2 - - # negative eigenvalue: DomainError for real input, complex principal value otherwise - @test_throws DomainError logarithm(A) - @test_throws DomainError logarithm(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) - @test_throws DomainError logarithm(A; alg = MatrixFunctionViaEig(LAPACK_Simple())) - - logA = @constinferred logarithm(complex.(A)) - @test exp(logA) ≈ A - - # (numerically) zero eigenvalue: no logarithm exists - Asing = V * LinearAlgebra.Diagonal(T[0, 1, 2, 3]) * V' - Asing = (Asing + Asing') / 2 - @test_throws DomainError logarithm(Asing; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) - @test_throws DomainError logarithm(Asing; alg = MatrixFunctionViaEig(LAPACK_Simple())) - @test_throws DomainError logarithm(complex.(Asing); alg = MatrixFunctionViaEig(LAPACK_Simple())) -end - -@testset "logarithm! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - m = 54 - - data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) - A = Diagonal(data) - Ac = copy(A) - logA = LinearAlgebra.log(A) - - logA2 = @constinferred logarithm(A) - @test logA2 isa Diagonal - @test logA ≈ logA2 - @test A == Ac - - if T <: Real - @test_throws DomainError logarithm(Diagonal(-data)) - end - @test_throws DomainError logarithm(Diagonal(zeros(T, m))) -end - -# GPU tests -# --------- -# As for `squareroot`, the `DiagonalAlgorithm` kernel and the shared domain helpers are -# backend-generic, and `MatrixFunctionViaEigh` goes through the device `eigh_full!`. Both -# domain rejections (negative real axis, and the zero eigenvalue where no logarithm exists) -# are exercised, since `_check_nonzero_eigenvalues` is a reduction too. If any step fell back -# to scalar indexing these would error under GPUArrays' scalar-indexing guard. -function test_logarithm_gpu(ArrayT, T) - rng = StableRNG(123) - m = 54 - - # Diagonal fast path - data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) - @test Array(diagview(logarithm(Diagonal(ArrayT(data))))) ≈ - diagview(logarithm(Diagonal(data))) - - # hermitian positive definite, via the eigenvalue decomposition - X = randn(rng, T, m, m) - A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) - A_gpu = ArrayT(A) - alg_gpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A_gpu)) - alg_cpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A)) - @test Array(logarithm(A_gpu, alg_gpu)) ≈ logarithm(A, alg_cpu) - - # domain handling: negative real axis, and no logarithm at a zero eigenvalue - T <: Real && @test_throws DomainError logarithm(Diagonal(ArrayT(-data))) - @test_throws DomainError logarithm(Diagonal(ArrayT(zeros(T, m)))) - return nothing -end - -if CUDA.functional() - @testset "logarithm on CUDA for T = $T" for T in BLASFloats - test_logarithm_gpu(CuArray, T) - end -end - -if AMDGPU.functional() - @testset "logarithm on AMDGPU for T = $T" for T in BLASFloats - test_logarithm_gpu(ROCArray, T) - end -end diff --git a/test/matrixfunctions/exponential.jl b/test/matrixfunctions/exponential.jl new file mode 100644 index 000000000..a16031ffa --- /dev/null +++ b/test/matrixfunctions/exponential.jl @@ -0,0 +1,114 @@ +using MatrixAlgebraKit +using LinearAlgebra: Diagonal +using MatrixAlgebraKit: GLA, GS +using CUDA, AMDGPU +using GenericSchur, GenericLinearAlgebra + +if @isdefined(fast_tests) && fast_tests + BLASFloats = (Float64, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +else + BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +end +# only the `Diagonal` fast path applies to these, as they have no `eig`/`eigh` support +DiagonalOnlyFloats = (Float16, ComplexF16) + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +m = 54 + +# CPU tests +# --------- +if !is_buildkite + # LAPACK algorithms: + for T in BLASFloats + TestSuite.seed_rng!(123) + LAPACK_EIG_ALGS = ( + MatrixFunctionViaLA(), + MatrixFunctionViaEig(QRIteration()), + MatrixFunctionViaTaylor(), + ) + LAPACK_EIGH_ALGS = (MatrixFunctionViaEigh(QRIteration()), MatrixFunctionViaEigh(DivideAndConquer())) + TestSuite.test_exponential(T, (m, m)) + TestSuite.test_exponential_algs(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_exponential_scaled(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_exponential_hermitian(T, (m, m), LAPACK_EIGH_ALGS) + TestSuite.test_exponential_taylor(T, (m, m)) + TestSuite.test_exponential_reference(T, (m, m)) + TestSuite.test_exponential_wrappers(T, (12, 12)) + # `eigh` rejects non-hermitian input rather than projecting it + TestSuite.test_exponential_domain(T, (m, m), LAPACK_EIGH_ALGS) + end + + # Generic floats: `eig` comes from GenericSchur, `eigh` from GenericLinearAlgebra. Both are + # loaded here, so name the driver explicitly instead of relying on `default_driver`. The + # native Taylor algorithm needs no LAPACK support and applies at arbitrary precision. + for T in GenericFloats + TestSuite.seed_rng!(123) + GS_ALGS = (MatrixFunctionViaEig(QRIteration(; driver = GS())), MatrixFunctionViaTaylor()) + GLA_ALGS = (MatrixFunctionViaEigh(QRIteration(; driver = GLA())),) + TestSuite.test_exponential_algs(T, (24, 24), GS_ALGS) + TestSuite.test_exponential_scaled(T, (24, 24), GS_ALGS) + TestSuite.test_exponential_hermitian(T, (24, 24), GLA_ALGS) + TestSuite.test_exponential_taylor(T, (24, 24)) + end + + # Diagonal: + for T in (BLASFloats..., GenericFloats..., DiagonalOnlyFloats...) + TestSuite.seed_rng!(123) + AT = Diagonal{T, Vector{T}} + test_spectrum = !(T in DiagonalOnlyFloats) + TestSuite.test_exponential(AT, m) + TestSuite.test_exponential_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_exponential_scaled(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_exponential_hermitian(AT, m, (DiagonalAlgorithm(),); test_spectrum) + TestSuite.test_exponential_reference(AT, m; test_hermitian = !(T in GenericFloats)) + end +end + +# CUDA tests +# ---------- +# Unlike the other matrix functions, a general dense matrix *is* supported on device, through the +# native `MatrixFunctionViaTaylor`. The `logarithm` roundtrip is unavailable there, so those calls +# rely on the `exp(A) * exp(-A) ≈ I` invariant instead. +if CUDA.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + CUDA_EIGH_ALGS = (MatrixFunctionViaEigh(Jacobi()), MatrixFunctionViaEigh(DivideAndConquer())) + TestSuite.test_exponential_algs(CuMatrix{T}, (m, m), (MatrixFunctionViaTaylor(),); test_roundtrip = false) + TestSuite.test_exponential_scaled(CuMatrix{T}, (m, m), (MatrixFunctionViaTaylor(),)) + TestSuite.test_exponential_taylor(CuMatrix{T}, (m, m)) + TestSuite.test_exponential_hermitian(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS) + TestSuite.test_exponential_domain(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS) + + AT = Diagonal{T, CuVector{T}} + TestSuite.test_exponential(AT, m) + TestSuite.test_exponential_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_exponential_scaled(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_exponential_hermitian(AT, m, (DiagonalAlgorithm(),)) + end +end + +# AMDGPU tests +# ------------ +if AMDGPU.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + ROC_EIGH_ALGS = (MatrixFunctionViaEigh(Jacobi()), MatrixFunctionViaEigh(DivideAndConquer())) + TestSuite.test_exponential_algs(ROCMatrix{T}, (m, m), (MatrixFunctionViaTaylor(),); test_roundtrip = false) + TestSuite.test_exponential_scaled(ROCMatrix{T}, (m, m), (MatrixFunctionViaTaylor(),)) + TestSuite.test_exponential_taylor(ROCMatrix{T}, (m, m)) + TestSuite.test_exponential_hermitian(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS) + TestSuite.test_exponential_domain(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS) + + AT = Diagonal{T, ROCVector{T}} + TestSuite.test_exponential(AT, m) + TestSuite.test_exponential_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_exponential_scaled(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_exponential_hermitian(AT, m, (DiagonalAlgorithm(),)) + end +end diff --git a/test/matrixfunctions/logarithm.jl b/test/matrixfunctions/logarithm.jl new file mode 100644 index 000000000..46a1a5b42 --- /dev/null +++ b/test/matrixfunctions/logarithm.jl @@ -0,0 +1,113 @@ +using MatrixAlgebraKit +using LinearAlgebra: Diagonal +using MatrixAlgebraKit: GLA, GS +using CUDA, AMDGPU +using GenericSchur, GenericLinearAlgebra + +if @isdefined(fast_tests) && fast_tests + BLASFloats = (Float64, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +else + BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +end +# only the `Diagonal` fast path applies to these, as they have no `eig`/`eigh` support +DiagonalOnlyFloats = (Float16, ComplexF16) + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +m = 54 +md = 4 # domain tests prescribe the full spectrum, so keep them small + +# CPU tests +# --------- +if !is_buildkite + # LAPACK algorithms: + for T in BLASFloats + TestSuite.seed_rng!(123) + LAPACK_EIG_ALGS = (MatrixFunctionViaLA(), MatrixFunctionViaEig(QRIteration())) + LAPACK_EIGH_ALGS = ( + MatrixFunctionViaEigh(QRIteration()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_logarithm(T, (m, m)) + TestSuite.test_logarithm_algs(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_logarithm_hermitian(T, (m, m), LAPACK_EIG_ALGS; exact_hermiticity = false) + TestSuite.test_logarithm_hermitian(T, (m, m), LAPACK_EIGH_ALGS) + TestSuite.test_logarithm_reference(T, (m, m)) + TestSuite.test_logarithm_domain(T, (md, md), (MatrixFunctionViaLA(),); supports_domain_atol = false) + TestSuite.test_logarithm_domain(T, (md, md), (MatrixFunctionViaEig(QRIteration()),)) + TestSuite.test_logarithm_domain(T, (md, md), LAPACK_EIGH_ALGS; hermitian_output = true) + end + + # Generic floats: `eig` comes from GenericSchur, `eigh` from GenericLinearAlgebra. Both are + # loaded here, so name the driver explicitly instead of relying on `default_driver`. + for T in GenericFloats + TestSuite.seed_rng!(123) + GS_ALGS = (MatrixFunctionViaEig(QRIteration(; driver = GS())),) + GLA_ALGS = (MatrixFunctionViaEigh(QRIteration(; driver = GLA())),) + TestSuite.test_logarithm_algs(T, (24, 24), GS_ALGS) + TestSuite.test_logarithm_hermitian(T, (24, 24), GS_ALGS; exact_hermiticity = false) + TestSuite.test_logarithm_hermitian(T, (24, 24), GLA_ALGS) + TestSuite.test_logarithm_domain(T, (md, md), GS_ALGS) + TestSuite.test_logarithm_domain(T, (md, md), GLA_ALGS; hermitian_output = true) + end + + # Diagonal: + for T in (BLASFloats..., GenericFloats..., DiagonalOnlyFloats...) + TestSuite.seed_rng!(123) + AT = Diagonal{T, Vector{T}} + test_spectrum = !(T in DiagonalOnlyFloats) + TestSuite.test_logarithm(AT, m) + TestSuite.test_logarithm_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_logarithm_hermitian(AT, m, (DiagonalAlgorithm(),); test_spectrum) + TestSuite.test_logarithm_reference(AT, m; test_hermitian = !(T in GenericFloats)) + TestSuite.test_logarithm_domain(AT, md, (DiagonalAlgorithm(),)) + end +end + +# CUDA tests +# ---------- +# General dense matrices are not supported on device: `MatrixFunctionViaLA` would call LAPACK on +# device memory, and `MatrixFunctionViaEig` scalar-indexes in its `lu!`-based solve. Hermitian +# input via `eigh` and the `Diagonal` fast path are both backend-generic. +if CUDA.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + CUDA_EIGH_ALGS = ( + MatrixFunctionViaEigh(Jacobi()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_logarithm_hermitian(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS) + TestSuite.test_logarithm_domain(CuMatrix{T}, (md, md), CUDA_EIGH_ALGS; hermitian_output = true) + + AT = Diagonal{T, CuVector{T}} + TestSuite.test_logarithm(AT, m) + TestSuite.test_logarithm_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_logarithm_hermitian(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_logarithm_domain(AT, md, (DiagonalAlgorithm(),)) + end +end + +# AMDGPU tests +# ------------ +if AMDGPU.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + ROC_EIGH_ALGS = ( + MatrixFunctionViaEigh(Jacobi()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_logarithm_hermitian(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS) + TestSuite.test_logarithm_domain(ROCMatrix{T}, (md, md), ROC_EIGH_ALGS; hermitian_output = true) + + AT = Diagonal{T, ROCVector{T}} + TestSuite.test_logarithm(AT, m) + TestSuite.test_logarithm_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_logarithm_hermitian(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_logarithm_domain(AT, md, (DiagonalAlgorithm(),)) + end +end diff --git a/test/matrixfunctions/power.jl b/test/matrixfunctions/power.jl new file mode 100644 index 000000000..208a48607 --- /dev/null +++ b/test/matrixfunctions/power.jl @@ -0,0 +1,113 @@ +using MatrixAlgebraKit +using LinearAlgebra: Diagonal +using MatrixAlgebraKit: GLA, GS +using CUDA, AMDGPU +using GenericSchur, GenericLinearAlgebra + +if @isdefined(fast_tests) && fast_tests + BLASFloats = (Float64, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +else + BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +end +# only the `Diagonal` fast path applies to these, as they have no `eig`/`eigh` support +DiagonalOnlyFloats = (Float16, ComplexF16) + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +m = 54 +md = 4 # domain tests prescribe the full spectrum, so keep them small + +# CPU tests +# --------- +if !is_buildkite + # LAPACK algorithms: + for T in BLASFloats + TestSuite.seed_rng!(123) + LAPACK_EIG_ALGS = (MatrixFunctionViaLA(), MatrixFunctionViaEig(QRIteration())) + LAPACK_EIGH_ALGS = ( + MatrixFunctionViaEigh(QRIteration()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_power(T, (m, m)) + TestSuite.test_power_algs(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_power_hermitian(T, (m, m), LAPACK_EIG_ALGS; exact_hermiticity = false) + TestSuite.test_power_hermitian(T, (m, m), LAPACK_EIGH_ALGS) + TestSuite.test_power_reference(T, (m, m)) + TestSuite.test_power_domain(T, (md, md), (MatrixFunctionViaLA(),); supports_domain_atol = false) + TestSuite.test_power_domain(T, (md, md), (MatrixFunctionViaEig(QRIteration()),)) + TestSuite.test_power_domain(T, (md, md), LAPACK_EIGH_ALGS; hermitian_output = true) + end + + # Generic floats: `eig` comes from GenericSchur, `eigh` from GenericLinearAlgebra. Both are + # loaded here, so name the driver explicitly instead of relying on `default_driver`. + for T in GenericFloats + TestSuite.seed_rng!(123) + GS_ALGS = (MatrixFunctionViaEig(QRIteration(; driver = GS())),) + GLA_ALGS = (MatrixFunctionViaEigh(QRIteration(; driver = GLA())),) + TestSuite.test_power_algs(T, (24, 24), GS_ALGS) + TestSuite.test_power_hermitian(T, (24, 24), GS_ALGS; exact_hermiticity = false) + TestSuite.test_power_hermitian(T, (24, 24), GLA_ALGS) + TestSuite.test_power_domain(T, (md, md), GS_ALGS) + TestSuite.test_power_domain(T, (md, md), GLA_ALGS; hermitian_output = true) + end + + # Diagonal: + for T in (BLASFloats..., GenericFloats..., DiagonalOnlyFloats...) + TestSuite.seed_rng!(123) + AT = Diagonal{T, Vector{T}} + test_spectrum = !(T in DiagonalOnlyFloats) + TestSuite.test_power(AT, m) + TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),); test_spectrum) + TestSuite.test_power_reference(AT, m; test_hermitian = !(T in GenericFloats)) + TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) + end +end + +# CUDA tests +# ---------- +# General dense matrices are not supported on device: `MatrixFunctionViaLA` would call LAPACK on +# device memory, and `MatrixFunctionViaEig` scalar-indexes in its `lu!`-based solve. Hermitian +# input via `eigh` and the `Diagonal` fast path are both backend-generic. +if CUDA.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + CUDA_EIGH_ALGS = ( + MatrixFunctionViaEigh(Jacobi()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_power_hermitian(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS) + TestSuite.test_power_domain(CuMatrix{T}, (md, md), CUDA_EIGH_ALGS; hermitian_output = true) + + AT = Diagonal{T, CuVector{T}} + TestSuite.test_power(AT, m) + TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) + end +end + +# AMDGPU tests +# ------------ +if AMDGPU.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + ROC_EIGH_ALGS = ( + MatrixFunctionViaEigh(Jacobi()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_power_hermitian(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS) + TestSuite.test_power_domain(ROCMatrix{T}, (md, md), ROC_EIGH_ALGS; hermitian_output = true) + + AT = Diagonal{T, ROCVector{T}} + TestSuite.test_power(AT, m) + TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) + end +end diff --git a/test/matrixfunctions/squareroot.jl b/test/matrixfunctions/squareroot.jl new file mode 100644 index 000000000..a800a751d --- /dev/null +++ b/test/matrixfunctions/squareroot.jl @@ -0,0 +1,112 @@ +using MatrixAlgebraKit +using LinearAlgebra: Diagonal +using MatrixAlgebraKit: GLA, GS +using CUDA, AMDGPU +using GenericSchur, GenericLinearAlgebra + +if @isdefined(fast_tests) && fast_tests + BLASFloats = (Float64, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +else + BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) + GenericFloats = (BigFloat, Complex{BigFloat}) +end +# only the `Diagonal` fast path applies to these, as they have no `eig`/`eigh` support +DiagonalOnlyFloats = (Float16, ComplexF16) + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +m = 54 +md = 4 # domain tests prescribe the full spectrum, so keep them small + +# CPU tests +# --------- +if !is_buildkite + # LAPACK algorithms: + for T in BLASFloats + TestSuite.seed_rng!(123) + LAPACK_EIG_ALGS = (MatrixFunctionViaLA(), MatrixFunctionViaEig(QRIteration())) + LAPACK_EIGH_ALGS = ( + MatrixFunctionViaEigh(QRIteration()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_squareroot(T, (m, m)) + TestSuite.test_squareroot_algs(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_squareroot_hermitian(T, (m, m), LAPACK_EIG_ALGS; exact_hermiticity = false) + TestSuite.test_squareroot_hermitian(T, (m, m), LAPACK_EIGH_ALGS) + TestSuite.test_squareroot_reference(T, (m, m)) + TestSuite.test_squareroot_domain(T, (md, md), LAPACK_EIG_ALGS) + TestSuite.test_squareroot_domain(T, (md, md), LAPACK_EIGH_ALGS; hermitian_output = true) + end + + # Generic floats: `eig` comes from GenericSchur, `eigh` from GenericLinearAlgebra. Both are + # loaded here, so name the driver explicitly instead of relying on `default_driver`. + for T in GenericFloats + TestSuite.seed_rng!(123) + GS_ALGS = (MatrixFunctionViaEig(QRIteration(; driver = GS())),) + GLA_ALGS = (MatrixFunctionViaEigh(QRIteration(; driver = GLA())),) + TestSuite.test_squareroot_algs(T, (24, 24), GS_ALGS) + TestSuite.test_squareroot_hermitian(T, (24, 24), GS_ALGS; exact_hermiticity = false) + TestSuite.test_squareroot_hermitian(T, (24, 24), GLA_ALGS) + TestSuite.test_squareroot_domain(T, (md, md), GS_ALGS) + TestSuite.test_squareroot_domain(T, (md, md), GLA_ALGS; hermitian_output = true) + end + + # Diagonal: + for T in (BLASFloats..., GenericFloats..., DiagonalOnlyFloats...) + TestSuite.seed_rng!(123) + AT = Diagonal{T, Vector{T}} + test_spectrum = !(T in DiagonalOnlyFloats) + TestSuite.test_squareroot(AT, m) + TestSuite.test_squareroot_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_squareroot_hermitian(AT, m, (DiagonalAlgorithm(),); test_spectrum) + TestSuite.test_squareroot_reference(AT, m; test_hermitian = !(T in GenericFloats)) + TestSuite.test_squareroot_domain(AT, md, (DiagonalAlgorithm(),)) + end +end + +# CUDA tests +# ---------- +# General dense matrices are not supported on device: `MatrixFunctionViaLA` would call LAPACK on +# device memory, and `MatrixFunctionViaEig` scalar-indexes in its `lu!`-based solve. Hermitian +# input via `eigh` and the `Diagonal` fast path are both backend-generic. +if CUDA.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + CUDA_EIGH_ALGS = ( + MatrixFunctionViaEigh(Jacobi()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_squareroot_hermitian(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS) + TestSuite.test_squareroot_domain(CuMatrix{T}, (md, md), CUDA_EIGH_ALGS; hermitian_output = true) + + AT = Diagonal{T, CuVector{T}} + TestSuite.test_squareroot(AT, m) + TestSuite.test_squareroot_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_squareroot_hermitian(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_squareroot_domain(AT, md, (DiagonalAlgorithm(),)) + end +end + +# AMDGPU tests +# ------------ +if AMDGPU.functional() + for T in BLASFloats + TestSuite.seed_rng!(123) + ROC_EIGH_ALGS = ( + MatrixFunctionViaEigh(Jacobi()), + MatrixFunctionViaEigh(DivideAndConquer()), + ) + TestSuite.test_squareroot_hermitian(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS) + TestSuite.test_squareroot_domain(ROCMatrix{T}, (md, md), ROC_EIGH_ALGS; hermitian_output = true) + + AT = Diagonal{T, ROCVector{T}} + TestSuite.test_squareroot(AT, m) + TestSuite.test_squareroot_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_squareroot_hermitian(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_squareroot_domain(AT, md, (DiagonalAlgorithm(),)) + end +end diff --git a/test/power.jl b/test/power.jl deleted file mode 100644 index 6322c784d..000000000 --- a/test/power.jl +++ /dev/null @@ -1,177 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using CUDA, AMDGPU - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) -GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) - -@testset "power! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - # spectrum inside a disk around 1, away from the negative real axis and zero - A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) - Ac = copy(A) - - # integer powers - @testset "integer p = $p" for p in (0, 1, 3, -2) - powA = A^p - powA2 = @constinferred power(A, p) - @test powA ≈ powA2 - @test A == Ac - - algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) - @testset "algorithm $alg" for alg in algs - powA2 = @constinferred power(A, p, alg) - @test powA ≈ powA2 - @test A == Ac - end - end - - # fractional powers - @testset "fractional p = $p" for p in (0.5, 0.75, -0.25, 3.5) - powA = A^p - powA2 = @constinferred power(A, p) - @test powA ≈ powA2 - @test A == Ac - - algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) - @testset "algorithm $alg" for alg in algs - powA2 = @constinferred power(A, p, alg) - @test powA ≈ powA2 - @test A == Ac - end - end - - @test power(A, 0.5) ≈ squareroot(A) - @test power(A, 3.0) ≈ power(A, 3) - @test_throws DomainError power(A, 0.5; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) -end - -@testset "power! for hermitian T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - X = randn(rng, T, m, m) - A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) - Ac = copy(A) - - algs = ( - MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()), - MatrixFunctionViaEig(LAPACK_Simple()), - ) - @testset "p = $p, algorithm $alg" for p in (2, -1, 0.5, -0.5), alg in algs - powA = LinearAlgebra.Hermitian(A)^p - powA2 = @constinferred power(A, p, alg) - @test powA ≈ powA2 - @test A == Ac - end - - @test LinearAlgebra.ishermitian(power(A, 0.5, MatrixFunctionViaEigh(LAPACK_QRIteration()))) -end - -@testset "power! domain handling for T = $T" for T in (Float32, Float64) - rng = StableRNG(123) - m = 4 - - X = randn(rng, T, m, m) - V = Matrix(LinearAlgebra.qr(X).Q) - A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V' - A = (A + A') / 2 - - # integer powers of matrices with negative eigenvalues are fine - @test power(A, 3) ≈ A * A * A - @test power(A, 3, MatrixFunctionViaEigh(LAPACK_QRIteration())) ≈ A * A * A - @test power(A, 3, MatrixFunctionViaEig(LAPACK_Simple())) ≈ A * A * A - - # fractional powers require staying off the negative real axis for real input - @test_throws DomainError power(A, 0.5) - @test_throws DomainError power(A, 0.5, MatrixFunctionViaEigh(LAPACK_QRIteration())) - @test_throws DomainError power(A, 0.5, MatrixFunctionViaEig(LAPACK_Simple())) - powA = @constinferred power(complex.(A), 0.5) - @test powA * powA ≈ A - - # (numerically) singular matrices with negative fractional powers - Asing = V * LinearAlgebra.Diagonal(T[0, 1, 2, 3]) * V' - Asing = (Asing + Asing') / 2 - @test_throws DomainError power(Asing, -0.5, MatrixFunctionViaEigh(LAPACK_QRIteration())) - @test_throws DomainError power(Asing, -0.5, MatrixFunctionViaEig(LAPACK_Simple())) -end - -@testset "power! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - m = 54 - - data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) - A = Diagonal(data) - Ac = copy(A) - - @testset "p = $p" for p in (2, -1, 0.5) - powA = A^p - powA2 = @constinferred power(A, p) - @test powA2 isa Diagonal - @test powA ≈ powA2 - @test A == Ac - end - - if T <: Real - @test_throws DomainError power(Diagonal(-data), 0.5) - @test power(Diagonal(-data), 2) ≈ Diagonal(data .^ 2) - end - @test_throws SingularException power(Diagonal(zeros(T, m)), -1) - @test_throws DomainError power(Diagonal(zeros(T, m)), -0.5) -end - -# GPU tests -# --------- -# As for `squareroot`, the `DiagonalAlgorithm` kernel and the shared domain helpers are -# backend-generic, and `MatrixFunctionViaEigh` goes through the device `eigh_full!`. Both the -# integer branch (which only checks for singularity) and the fractional branch (which clamps -# and rejects out-of-domain eigenvalues) are exercised. If any step fell back to scalar -# indexing these would error under GPUArrays' scalar-indexing guard. -function test_power_gpu(ArrayT, T) - rng = StableRNG(123) - m = 54 - - # Diagonal fast path, integer and fractional exponents - data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T)) - for p in (2, -1, 0.5, -0.25) - @test Array(diagview(power(Diagonal(ArrayT(data)), p))) ≈ - diagview(power(Diagonal(data), p)) - end - - # hermitian positive definite, via the eigenvalue decomposition - X = randn(rng, T, m, m) - A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) - A_gpu = ArrayT(A) - alg_gpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A_gpu)) - alg_cpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A)) - for p in (2, -1, 0.5, -0.5) - @test Array(power(A_gpu, p, alg_gpu)) ≈ power(A, p, alg_cpu) - end - - # domain handling: integer powers ignore the domain, fractional ones do not - if T <: Real - @test Array(diagview(power(Diagonal(ArrayT(-data)), 2))) ≈ data .^ 2 - @test_throws DomainError power(Diagonal(ArrayT(-data)), 0.5) - end - @test_throws SingularException power(Diagonal(ArrayT(zeros(T, m))), -1) - @test_throws DomainError power(Diagonal(ArrayT(zeros(T, m))), -0.5) - return nothing -end - -if CUDA.functional() - @testset "power on CUDA for T = $T" for T in BLASFloats - test_power_gpu(CuArray, T) - end -end - -if AMDGPU.functional() - @testset "power on AMDGPU for T = $T" for T in BLASFloats - test_power_gpu(ROCArray, T) - end -end diff --git a/test/squareroot.jl b/test/squareroot.jl deleted file mode 100644 index 478548606..000000000 --- a/test/squareroot.jl +++ /dev/null @@ -1,154 +0,0 @@ -using MatrixAlgebraKit -using Test -using TestExtras -using StableRNGs -using MatrixAlgebraKit: diagview -using LinearAlgebra -using CUDA, AMDGPU - -BLASFloats = (Float32, Float64, ComplexF32, ComplexF64) -GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat}) - -@testset "squareroot! for T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - # spectrum inside a disk around 1, away from the negative real axis - A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m)) - Ac = copy(A) - sqrtA = LinearAlgebra.sqrt(A) - - sqrtA2 = @constinferred squareroot(A) - @test sqrtA ≈ sqrtA2 - @test sqrtA2 * sqrtA2 ≈ A - @test A == Ac - - algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple())) - @testset "algorithm $alg" for alg in algs - sqrtA2 = @constinferred squareroot(A, alg) - @test sqrtA ≈ sqrtA2 - @test A == Ac - end - - @test_throws DomainError squareroot(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) -end - -@testset "squareroot! for hermitian T = $T" for T in BLASFloats - rng = StableRNG(123) - m = 54 - - X = randn(rng, T, m, m) - A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) - Ac = copy(A) - sqrtA = LinearAlgebra.sqrt(LinearAlgebra.Hermitian(A)) - - algs = ( - MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()), - MatrixFunctionViaEig(LAPACK_Simple()), - ) - @testset "algorithm $alg" for alg in algs - sqrtA2 = @constinferred squareroot(A, alg) - @test sqrtA ≈ sqrtA2 - @test A == Ac - end - - @test LinearAlgebra.ishermitian(squareroot(A, MatrixFunctionViaEigh(LAPACK_QRIteration()))) -end - -@testset "squareroot! domain handling for T = $T" for T in (Float32, Float64) - rng = StableRNG(123) - m = 4 - - # genuinely negative eigenvalue: DomainError for real input, complex principal value otherwise - X = randn(rng, T, m, m) - V = Matrix(LinearAlgebra.qr(X).Q) - A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V' - A = (A + A') / 2 - - @test_throws DomainError squareroot(A) - @test_throws DomainError squareroot(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration())) - @test_throws DomainError squareroot(A; alg = MatrixFunctionViaEig(LAPACK_Simple())) - - sqrtA = @constinferred squareroot(complex.(A)) - @test sqrtA * sqrtA ≈ A - - # roundoff-scale negative eigenvalue: clamped to zero instead of throwing - Aclamp = V * LinearAlgebra.Diagonal(T[-10 * eps(T), 1, 2, 3]) * V' - Aclamp = (Aclamp + Aclamp') / 2 - @testset "algorithm $alg" for alg in ( - MatrixFunctionViaEigh(LAPACK_QRIteration()), - MatrixFunctionViaEig(LAPACK_Simple()), - ) - sqrtA2 = @constinferred squareroot(Aclamp, alg) - @test eltype(sqrtA2) == T - @test sqrtA2 * sqrtA2 ≈ Aclamp atol = sqrt(eps(T)) - end -end - -@testset "squareroot! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...) - rng = StableRNG(123) - m = 54 - - data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : randn(rng, T, m) - A = Diagonal(data) - Ac = copy(A) - sqrtA = LinearAlgebra.sqrt(A) - - sqrtA2 = @constinferred squareroot(A) - @test sqrtA2 isa Diagonal - @test sqrtA ≈ sqrtA2 - @test A == Ac - - if T <: Real - @test_throws DomainError squareroot(Diagonal(-data)) - # roundoff-scale negative entries are clamped - @test squareroot(Diagonal(T[-eps(T), 1])) ≈ Diagonal(T[0, 1]) - end -end - -# GPU tests -# --------- -# The `DiagonalAlgorithm` kernel and the shared domain helpers are written as reductions and -# broadcasts rather than scalar loops, so the same code runs on GPU; the `MatrixFunctionViaEigh` -# path additionally goes through the device `eigh_full!`. Compare against the CPU reference and -# exercise both domain outcomes, since the clamp and the throw are exactly what the reductions -# replaced. If any step fell back to scalar indexing these would error under GPUArrays' -# scalar-indexing guard. -function test_squareroot_gpu(ArrayT, T) - rng = StableRNG(123) - m = 54 - - # Diagonal fast path - data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : randn(rng, T, m) - @test Array(diagview(squareroot(Diagonal(ArrayT(data))))) ≈ - diagview(squareroot(Diagonal(data))) - - # hermitian positive definite, via the eigenvalue decomposition - X = randn(rng, T, m, m) - A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I)) - A_gpu = ArrayT(A) - alg_gpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A_gpu)) - alg_cpu = MatrixFunctionViaEigh(MatrixAlgebraKit.select_algorithm(eigh_full, A)) - sqrtA = squareroot(A_gpu, alg_gpu) - @test Array(sqrtA) ≈ squareroot(A, alg_cpu) - @test Array(sqrtA) * Array(sqrtA) ≈ A - - # domain handling: roundoff-scale negatives are clamped, genuine ones throw - if T <: Real - @test Array(diagview(squareroot(Diagonal(ArrayT(T[-eps(T), 1]))))) ≈ T[0, 1] - @test_throws DomainError squareroot(Diagonal(ArrayT(-data))) - end - return nothing -end - -if CUDA.functional() - @testset "squareroot on CUDA for T = $T" for T in BLASFloats - test_squareroot_gpu(CuArray, T) - end -end - -if AMDGPU.functional() - @testset "squareroot on AMDGPU for T = $T" for T in BLASFloats - test_squareroot_gpu(ROCArray, T) - end -end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 36ae68304..5dc29ec9a 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -102,6 +102,40 @@ function instantiate_rank_deficient_matrix(::Type{T}, sz; trunc = truncrank(div( return Diagonal(diag(mul!(A, V, C))) end +# Matrix functions such as `squareroot`, `logarithm` and fractional `power` are only defined +# for matrices whose spectrum avoids (part of) the negative real axis, so their tests need +# matrices with a controlled spectrum rather than the plain `randn` of `instantiate_matrix`. + +# Spectral radius at most one, so that `exponential` stays well inside the principal branch and +# `logarithm` inverts it. +function instantiate_smallnorm_matrix(T, sz) + A = instantiate_matrix(T, sz) + return A / norm(A) +end + +# Spectrum inside the unit disk around 1, i.e. clear of the negative real axis and of zero. +instantiate_offaxis_matrix(T, sz) = instantiate_smallnorm_matrix(T, sz) + I + +# Hermitian positive definite, so that the `eigh`-based algorithms apply. +function instantiate_posdef_matrix(T, sz) + A = instantiate_matrix(T, sz) + return project_hermitian!(A * A') + I +end + +# Hermitian with the prescribed (real) spectrum `λ`, for probing the domain boundary. +function instantiate_hermitian_spectrum(T, sz, λ) + A = instantiate_matrix(T, sz) + n = size(A, 1) + @assert length(λ) == n + dv = A isa Diagonal ? diagview(A) : A + Ddiag = similar(dv, eltype(A), n) + # convert on the host first: `copyto!` onto a device array does not convert eltypes + copyto!(Ddiag, convert(Vector{eltype(A)}, collect(λ))) + A isa Diagonal && return Diagonal(Ddiag) + V = instantiate_unitary(T, A, n) + return project_hermitian!(V * Diagonal(Ddiag) * V') +end + include("ad_utils.jl") include("projections.jl") @@ -117,6 +151,13 @@ include("decompositions/eigh.jl") include("decompositions/orthnull.jl") include("decompositions/svd.jl") +# Matrix functions +# ---------------- +include("matrixfunctions/exponential.jl") +include("matrixfunctions/squareroot.jl") +include("matrixfunctions/logarithm.jl") +include("matrixfunctions/power.jl") + # Mooncake # -------- include("mooncake/mooncake.jl") diff --git a/test/testsuite/matrixfunctions/exponential.jl b/test/testsuite/matrixfunctions/exponential.jl new file mode 100644 index 000000000..1c162094d --- /dev/null +++ b/test/testsuite/matrixfunctions/exponential.jl @@ -0,0 +1,172 @@ +using TestExtras +using LinearAlgebra: LinearAlgebra, I +using MatrixAlgebraKit: ishermitian + +# `exponential` has no restricted domain. Two invariants are used: +# +# * `exp(A) * exp(-A) ≈ I`, which holds for every algorithm and every backend, since `A` and +# `-A` commute. This is the primary check. +# * `logarithm(exp(A)) ≈ A`, which additionally needs `A` inside the principal branch — +# guaranteed by `instantiate_smallnorm_matrix`, as a spectral radius of at most one keeps +# every eigenvalue's imaginary part well inside `(-π, π]`. Gated behind `test_roundtrip`, +# since `logarithm` of a general matrix is not available on device. + +function test_exponential(T::Type, sz; test_roundtrip = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "exponential $summary_str" begin + A = instantiate_smallnorm_matrix(T, sz) + Ac = deepcopy(A) + + expA = @testinferred exponential(A) + @test eltype(expA) == eltype(A) + @test expA * exponential(-A) ≈ I + @test A == Ac + test_roundtrip && @test logarithm(expA) ≈ A + + # the in-place method may not be able to reuse the provided output + expA2 = @testinferred exponential!(deepcopy(A), deepcopy(expA)) + @test expA2 ≈ expA + end +end + +function test_exponential_algs(T::Type, sz, algs; test_roundtrip = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "exponential algorithm $alg $summary_str" for alg in algs + A = instantiate_smallnorm_matrix(T, sz) + Ac = deepcopy(A) + + expA = @testinferred exponential(A, alg) + @test eltype(expA) == eltype(A) + @test expA * exponential(-A, alg) ≈ I + @test A == Ac + test_roundtrip && @test logarithm(expA) ≈ A + end +end + +# The scaled entrypoint `exponential((τ, A))` computes `exp(τ * A)`. +function test_exponential_scaled(T::Type, sz, algs; kwargs...) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "exponential scaled algorithm $alg $summary_str" for alg in algs + A = instantiate_smallnorm_matrix(T, sz) + Ac = deepcopy(A) + + # both a scalar of the matrix eltype and a real one, to exercise promotion + @testset "τ::$(typeof(τ))" for τ in (randn(rng, eltype(T)), randn(rng, R)) + expτA = @testinferred exponential((τ, A), alg) + @test eltype(expτA) == eltype(A) + @test expτA ≈ exponential(τ * A, alg) + @test expτA * exponential((-τ, A), alg) ≈ I + @test A == Ac + end + end +end + +# See the corresponding comment in `squareroot.jl` for `exact_hermiticity`. `alg` is reused for +# the `logarithm` roundtrip, so this must not be called with `MatrixFunctionViaTaylor`, which +# only implements `exponential`. +function test_exponential_hermitian( + T::Type, sz, algs; + exact_hermiticity = true, test_spectrum = true, kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "exponential hermitian algorithm $alg $summary_str" for alg in algs + A = project_hermitian!(instantiate_smallnorm_matrix(T, sz)) + Ac = deepcopy(A) + + expA = @testinferred exponential(A, alg) + @test eltype(expA) == eltype(A) + @test expA * exponential(-A, alg) ≈ I + @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 + @test ishermitian(expA) + test_spectrum && @test eigh_vals(expA) ≈ exp.(eigh_vals(A)) + else + @test ishermitian(expA; rtol = precision(T)) + end + end +end + +# `MatrixFunctionViaTaylor` is a native scaling-and-squaring evaluation, so it is the only +# algorithm that applies to a general matrix on device and at arbitrary precision. Its balancing +# step is exercised with a badly-scaled similarity transform `Aᵢⱼ ← Aᵢⱼ sᵢ / sⱼ`. Dense input +# only, as the scaling would densify a `Diagonal`. +function test_exponential_taylor(T::Type, sz; kwargs...) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "exponential Taylor $summary_str" begin + A = instantiate_smallnorm_matrix(T, sz) + n = size(A, 1) + s = similar(A, R, n) + copyto!(s, exp10.(range(-R(3), R(3), length = n))) + Abad = A .* s ./ transpose(s) + + @testset "balance = $balance" for balance in (true, false) + alg = MatrixFunctionViaTaylor(; balance) + # `Abad` is deliberately ill-conditioned, so `exp(A) * exp(-A) ≈ I` loses too many + # digits to assert there; the balancing invariance below is the check that matters. + expA = @testinferred exponential(A, alg) + @test eltype(expA) == eltype(A) + @test expA * exponential(-A, alg) ≈ I + + for M in (A, Abad) + expM = @testinferred exponential(M, alg) + @test eltype(expM) == eltype(M) + # the balancing similarity must not change the result + @test expM ≈ exponential(M, MatrixFunctionViaTaylor(; balance = !balance)) + end + end + end +end + +# `MatrixFunctionViaEigh` requires hermitian input, and says so rather than silently projecting. +function test_exponential_domain(T::Type, sz, algs; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "exponential domain algorithm $alg $summary_str" for alg in algs + A = instantiate_smallnorm_matrix(T, sz) + @test_throws DomainError exponential(A, alg) + end +end + +# Input that is not a `Matrix`: the kernels must not assume a strided layout. +function test_exponential_wrappers(T::Type, sz; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "exponential non-Matrix input $summary_str" begin + A = instantiate_smallnorm_matrix(T, sz) + m = size(A, 1) + expA = exponential(A) + + wrappers = ( + ("view", B -> view(B, :, :)), + ("PermutedDimsArray", B -> PermutedDimsArray(permutedims(B), (2, 1))), + ("ReshapedArray", B -> reshape(view(vec(B), 1:(m * m)), m, m)), + ) + @testset "$name" for (name, wrap) in wrappers + W = wrap(deepcopy(A)) + @test !(W isa Matrix) + @test exponential!(W) ≈ expA + end + end +end + +# Cross-check against `LinearAlgebra`, which only applies to host arrays. Pass +# `test_hermitian = false` for generic eltypes, as `LinearAlgebra` has no matrix functions for a +# `Hermitian` wrapper outside the BLAS floats. +function test_exponential_reference(T::Type, sz; test_hermitian = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "exponential vs LinearAlgebra $summary_str" begin + A = instantiate_smallnorm_matrix(T, sz) + @test exponential(A) ≈ LinearAlgebra.exp(A) + + if test_hermitian + H = project_hermitian!(instantiate_smallnorm_matrix(T, sz)) + @test exponential(H) ≈ LinearAlgebra.exp(LinearAlgebra.Hermitian(H)) + end + end +end diff --git a/test/testsuite/matrixfunctions/logarithm.jl b/test/testsuite/matrixfunctions/logarithm.jl new file mode 100644 index 000000000..e94bc4968 --- /dev/null +++ b/test/testsuite/matrixfunctions/logarithm.jl @@ -0,0 +1,121 @@ +using TestExtras +using LinearAlgebra: LinearAlgebra, I +using MatrixAlgebraKit: ishermitian + +# `exponential` is the inverse of `logarithm` on the principal branch, and its default algorithm +# is backend-generic, so it serves as the reference invariant here. + +function test_logarithm(T::Type, sz; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "logarithm $summary_str" begin + A = instantiate_offaxis_matrix(T, sz) + Ac = deepcopy(A) + + logA = @testinferred logarithm(A) + @test eltype(logA) == eltype(A) + @test exponential(logA) ≈ A + @test A == Ac + + # the in-place method may not be able to reuse the provided output + logA2 = @testinferred logarithm!(deepcopy(A), deepcopy(logA)) + @test logA2 ≈ logA + end +end + +function test_logarithm_algs(T::Type, sz, algs; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "logarithm algorithm $alg $summary_str" for alg in algs + A = instantiate_offaxis_matrix(T, sz) + Ac = deepcopy(A) + + logA = @testinferred logarithm(A, alg) + @test eltype(logA) == eltype(A) + @test exponential(logA) ≈ A + @test A == Ac + end +end + +# See the corresponding comment in `squareroot.jl` for `exact_hermiticity`. +function test_logarithm_hermitian( + T::Type, sz, algs; + exact_hermiticity = true, test_spectrum = true, kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "logarithm hermitian algorithm $alg $summary_str" for alg in algs + A = instantiate_posdef_matrix(T, sz) + Ac = deepcopy(A) + + logA = @testinferred logarithm(A, alg) + @test eltype(logA) == eltype(A) + @test exponential(logA) ≈ A + @test A == Ac + + if exact_hermiticity + @test ishermitian(logA) + test_spectrum && @test eigh_vals(logA) ≈ log.(eigh_vals(A)) + else + @test ishermitian(logA; rtol = precision(T)) + end + end +end + +# Domain handling. `logarithm` is undefined both on the negative real axis and at zero, and the +# two interact: an eigenvalue that is negative only by roundoff gets clamped onto the boundary, +# i.e. to zero, where there is still no logarithm. So unlike `squareroot`, the clamp does not +# rescue such a matrix and a `DomainError` is still the correct outcome. +# +# `hermitian_output = true`: see `squareroot.jl`. +# `supports_domain_atol = false`: for `MatrixFunctionViaLA`, which inspects only the realness of +# the result rather than the spectrum, so it silently accepts singular input. +function test_logarithm_domain( + T::Type, sz, algs; + hermitian_output = false, supports_domain_atol = true, kwargs... + ) + R = real(eltype(T)) + n = sz isa Tuple ? first(sz) : sz + summary_str = testargs_summary(T, sz) + return @testset "logarithm domain algorithm $alg $summary_str" for alg in algs + # eigenvalue on the negative real axis + λ = collect(R, 1:n) + λ[1] = -one(R) + A = instantiate_hermitian_spectrum(T, sz, λ) + + if eltype(T) <: Real || hermitian_output + @test_throws DomainError logarithm(A, alg) + else + logA = @testinferred logarithm(A, alg) + @test exponential(logA) ≈ A + end + + supports_domain_atol || continue + + # (numerically) zero eigenvalue: no logarithm exists + λzero = collect(R, 1:n) + λzero[1] = zero(R) + @test_throws DomainError logarithm(instantiate_hermitian_spectrum(T, sz, λzero), alg) + + # roundoff-scale negative eigenvalue: clamped onto the boundary, which is still singular + λtiny = collect(R, 1:n) + λtiny[1] = -10 * eps(R) + @test_throws DomainError logarithm(instantiate_hermitian_spectrum(T, sz, λtiny), alg) + end +end + +# Cross-check against `LinearAlgebra`, which only applies to host arrays. Pass +# `test_hermitian = false` for generic eltypes, as `LinearAlgebra` has no matrix functions for a +# `Hermitian` wrapper outside the BLAS floats. +# +# Note this must not be called for *dense* generic eltypes at all: +# `LinearAlgebra.log(::UpperTriangular{BigFloat})` does not terminate. +function test_logarithm_reference(T::Type, sz; test_hermitian = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "logarithm vs LinearAlgebra $summary_str" begin + A = instantiate_offaxis_matrix(T, sz) + @test logarithm(A) ≈ LinearAlgebra.log(A) + + if test_hermitian + H = instantiate_posdef_matrix(T, sz) + @test logarithm(H) ≈ LinearAlgebra.log(LinearAlgebra.Hermitian(H)) + end + end +end diff --git a/test/testsuite/matrixfunctions/power.jl b/test/testsuite/matrixfunctions/power.jl new file mode 100644 index 000000000..065e7159e --- /dev/null +++ b/test/testsuite/matrixfunctions/power.jl @@ -0,0 +1,184 @@ +using TestExtras +using LinearAlgebra: LinearAlgebra, I, SingularException +using MatrixAlgebraKit: ishermitian + +# `power` takes the exponent as a second positional argument, and splits into an integer branch +# (repeated multiplication, defined for any square matrix) and a fractional branch (the principal +# power, restricted to the domain). Both are exercised throughout. +# +# Integer and fractional exponents are always iterated in separate, uniformly typed testsets: +# a tuple mixing `Int` and floats would infer as a `Union` and defeat `@testinferred`. + +function test_power(T::Type, sz; ps = (0, 1, 2, -1, 3), qs = (1 // 2, 3 // 4, -1 // 4), kwargs...) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "power $summary_str" begin + A = instantiate_offaxis_matrix(T, sz) + Ac = deepcopy(A) + + @testset "integer p = $p" for p in ps + powA = @testinferred power(A, p) + @test eltype(powA) == eltype(A) + @test A == Ac + # an integral exponent is recognised at runtime, as `Base.^` does + @test powA ≈ power(A, convert(R, p)) + end + + @testset "fractional p = $q" for q in qs + powA = @testinferred power(A, convert(R, q)) + @test eltype(powA) == eltype(A) + @test A == Ac + end + + @test power(A, 0) ≈ I + @test power(A, 1) ≈ A + @test power(A, 2) ≈ A * A + @test power(A, 3) ≈ A * A * A + @test power(A, -1) * A ≈ I + + # `squareroot` is the `p = 1/2` case, and exponents add + @test power(A, one(R) / 2) ≈ squareroot(A) + @test power(A, one(R) / 4) * power(A, one(R) / 4) ≈ power(A, one(R) / 2) + + # the in-place method may not be able to reuse the provided output + powA = power(A, 2) + @test @testinferred(power!(deepcopy(A), 2, deepcopy(powA))) ≈ powA + end +end + +function test_power_algs(T::Type, sz, algs; ps = (2, -1), qs = (1 // 2, -1 // 4), kwargs...) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "power algorithm $alg $summary_str" for alg in algs + A = instantiate_offaxis_matrix(T, sz) + Ac = deepcopy(A) + + @testset "integer p = $p" for p in ps + powA = @testinferred power(A, p, alg) + @test eltype(powA) == eltype(A) + @test A == Ac + end + + @testset "fractional p = $q" for q in qs + powA = @testinferred power(A, convert(R, q), alg) + @test eltype(powA) == eltype(A) + @test A == Ac + end + + @test power(A, 2, alg) ≈ A * A + @test power(A, -1, alg) * A ≈ I + @test power(A, one(R) / 2, alg) ≈ squareroot(A, alg) + end +end + +# See the corresponding comment in `squareroot.jl` for `exact_hermiticity`. +function test_power_hermitian( + T::Type, sz, algs; + exact_hermiticity = true, test_spectrum = true, kwargs... + ) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "power hermitian algorithm $alg $summary_str" for alg in algs + A = instantiate_posdef_matrix(T, sz) + Ac = deepcopy(A) + + check = function (powA, p) + @test eltype(powA) == eltype(A) + @test A == Ac + if exact_hermiticity + @test ishermitian(powA) + # `power` maps the spectrum elementwise, up to the reordering a negative + # exponent induces + test_spectrum && @test eigh_vals(powA) ≈ sort(eigh_vals(A) .^ p) + else + @test ishermitian(powA; rtol = precision(T)) + end + return nothing + end + + @testset "integer p = $p" for p in (2, -1) + check(@testinferred(power(A, p, alg)), p) + end + @testset "fractional p = $q" for q in (1 // 2, -1 // 2) + p = convert(R, q) + check(@testinferred(power(A, p, alg)), p) + end + end +end + +# Domain handling. Integer exponents are defined for any square matrix and must ignore the +# domain entirely; fractional exponents are principal powers and must respect it. Unlike +# `logarithm`, clamping a roundoff-negative eigenvalue onto zero is harmless for a positive +# fractional exponent, since `0^p` is well defined there. +# +# `hermitian_output = true` and `supports_domain_atol = false`: see `squareroot.jl` and +# `logarithm.jl`. `test_singular = true` only where the zero eigenvalue is exactly representable +# (the `Diagonal` path); the `eig`-based kernels test `any(iszero, λ)` on *computed* eigenvalues, +# which a numerically singular matrix does not satisfy. +function test_power_domain( + T::Type, sz, algs; + hermitian_output = false, supports_domain_atol = true, test_singular = false, kwargs... + ) + R = real(eltype(T)) + n = sz isa Tuple ? first(sz) : sz + half = one(R) / 2 + summary_str = testargs_summary(T, sz) + return @testset "power domain algorithm $alg $summary_str" for alg in algs + # eigenvalue on the negative real axis + λ = collect(R, 1:n) + λ[1] = -one(R) + A = instantiate_hermitian_spectrum(T, sz, λ) + + # integer exponents are unaffected by the domain + @test power(A, 2, alg) ≈ A * A + @test power(A, 3, alg) ≈ A * A * A + + if eltype(T) <: Real || hermitian_output + @test_throws DomainError power(A, half, alg) + else + powA = @testinferred power(A, half, alg) + @test powA * powA ≈ A + end + + supports_domain_atol || continue + + # roundoff-scale negative eigenvalue: clamped onto the boundary, which is fine for p > 0 + λtiny = collect(R, 1:n) + λtiny[1] = -10 * eps(R) + Atiny = instantiate_hermitian_spectrum(T, sz, λtiny) + powAtiny = @testinferred power(Atiny, half, alg) + @test eltype(powAtiny) == eltype(Atiny) + @test powAtiny * powAtiny ≈ Atiny atol = sqrt(eps(R)) + + # a negative fractional exponent additionally requires a nonzero spectrum + λzero = collect(R, 1:n) + λzero[1] = zero(R) + Azero = instantiate_hermitian_spectrum(T, sz, λzero) + @test_throws DomainError power(Azero, -half, alg) + test_singular && @test_throws SingularException power(Azero, -1, alg) + end +end + +# Cross-check against `LinearAlgebra`, which only applies to host arrays. Pass +# `test_hermitian = false` for generic eltypes, as `LinearAlgebra` has no matrix functions for a +# `Hermitian` wrapper outside the BLAS floats. +function test_power_reference(T::Type, sz; test_hermitian = true, kwargs...) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "power vs LinearAlgebra $summary_str" begin + A = instantiate_offaxis_matrix(T, sz) + @testset "integer p = $p" for p in (2, -1) + @test power(A, p) ≈ A^p + end + @testset "fractional p = $q" for q in (1 // 2, -1 // 4) + p = convert(R, q) + @test power(A, p) ≈ A^p + end + + if test_hermitian + H = instantiate_posdef_matrix(T, sz) + @test power(H, 2) ≈ LinearAlgebra.Hermitian(H)^2 + @test power(H, one(R) / 2) ≈ LinearAlgebra.Hermitian(H)^(one(R) / 2) + end + end +end diff --git a/test/testsuite/matrixfunctions/squareroot.jl b/test/testsuite/matrixfunctions/squareroot.jl new file mode 100644 index 000000000..1e2b2b453 --- /dev/null +++ b/test/testsuite/matrixfunctions/squareroot.jl @@ -0,0 +1,119 @@ +using TestExtras +using LinearAlgebra: LinearAlgebra, I +using MatrixAlgebraKit: ishermitian + +# The assertions below are invariants rather than comparisons against a reference +# implementation, so that the same bodies apply on GPU and to downstream array types. +# `test_squareroot_reference` is the exception and is host-only by construction. + +function test_squareroot(T::Type, sz; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "squareroot $summary_str" begin + A = instantiate_offaxis_matrix(T, sz) + Ac = deepcopy(A) + + sqrtA = @testinferred squareroot(A) + @test eltype(sqrtA) == eltype(A) + @test sqrtA * sqrtA ≈ A + @test A == Ac + + # the in-place method may not be able to reuse the provided output + sqrtA2 = @testinferred squareroot!(deepcopy(A), deepcopy(sqrtA)) + @test sqrtA2 ≈ sqrtA + + # `squareroot` is the `p = 1/2` case of `power` + @test sqrtA ≈ power(A, one(real(eltype(A))) / 2) + end +end + +function test_squareroot_algs(T::Type, sz, algs; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "squareroot algorithm $alg $summary_str" for alg in algs + A = instantiate_offaxis_matrix(T, sz) + Ac = deepcopy(A) + + sqrtA = @testinferred squareroot(A, alg) + @test eltype(sqrtA) == eltype(A) + @test sqrtA * sqrtA ≈ A + @test A == Ac + end +end + +# Hermitian positive definite input. +# +# The `eigh`-based kernels build the result as a symmetric product (`_mul_herm!`), so it is +# hermitian to the last bit; pass `exact_hermiticity = false` for algorithms that route through +# the general `eig` path instead, whose output is only approximately hermitian. The elementwise +# spectrum check needs exact hermiticity, since `eigh_vals` rejects anything else. +function test_squareroot_hermitian( + T::Type, sz, algs; + exact_hermiticity = true, test_spectrum = true, kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "squareroot hermitian algorithm $alg $summary_str" for alg in algs + A = instantiate_posdef_matrix(T, sz) + Ac = deepcopy(A) + + sqrtA = @testinferred squareroot(A, alg) + @test eltype(sqrtA) == eltype(A) + @test sqrtA * sqrtA ≈ A + @test A == Ac + + if exact_hermiticity + @test ishermitian(sqrtA) + # the square root maps the spectrum elementwise + test_spectrum && @test eigh_vals(sqrtA) ≈ sqrt.(eigh_vals(A)) + else + @test ishermitian(sqrtA; rtol = precision(T)) + end + end +end + +# Domain handling: a matrix whose spectrum reaches the negative real axis has a complex +# principal square root, which a type-stable real output cannot represent. +# +# Pass `hermitian_output = true` for algorithms that promise a hermitian result (i.e. the +# `eigh`-based ones). Those must reject a negative eigenvalue whatever the scalar type, since +# the square root of a hermitian matrix with a negative eigenvalue is not hermitian. +function test_squareroot_domain(T::Type, sz, algs; hermitian_output = false, kwargs...) + R = real(eltype(T)) + n = sz isa Tuple ? first(sz) : sz + summary_str = testargs_summary(T, sz) + return @testset "squareroot domain algorithm $alg $summary_str" for alg in algs + # genuinely negative eigenvalue + λ = collect(R, 1:n) + λ[1] = -one(R) + A = instantiate_hermitian_spectrum(T, sz, λ) + + if eltype(T) <: Real || hermitian_output + @test_throws DomainError squareroot(A, alg) + else + sqrtA = @testinferred squareroot(A, alg) + @test sqrtA * sqrtA ≈ A + end + + # roundoff-scale negative eigenvalue: clamped onto the boundary rather than rejected + λclamp = collect(R, 1:n) + λclamp[1] = -10 * eps(R) + Aclamp = instantiate_hermitian_spectrum(T, sz, λclamp) + sqrtAclamp = @testinferred squareroot(Aclamp, alg) + @test eltype(sqrtAclamp) == eltype(Aclamp) + @test sqrtAclamp * sqrtAclamp ≈ Aclamp atol = sqrt(eps(R)) + end +end + +# Cross-check against `LinearAlgebra`, which only applies to host arrays. Pass +# `test_hermitian = false` for generic eltypes, as `LinearAlgebra` has no matrix functions for a +# `Hermitian` wrapper outside the BLAS floats. +function test_squareroot_reference(T::Type, sz; test_hermitian = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "squareroot vs LinearAlgebra $summary_str" begin + A = instantiate_offaxis_matrix(T, sz) + @test squareroot(A) ≈ LinearAlgebra.sqrt(A) + + if test_hermitian + H = instantiate_posdef_matrix(T, sz) + @test squareroot(H) ≈ LinearAlgebra.sqrt(LinearAlgebra.Hermitian(H)) + end + end +end From 6427635fec8bba5ab317cf30d89be81867a665af Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 12:44:49 -0400 Subject: [PATCH 12/22] review comment --- docs/src/user_interface/matrix_functions.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/user_interface/matrix_functions.md b/docs/src/user_interface/matrix_functions.md index 71570ba3c..0888861cd 100644 --- a/docs/src/user_interface/matrix_functions.md +++ b/docs/src/user_interface/matrix_functions.md @@ -42,7 +42,8 @@ MatrixAlgebraKit.MatrixFunctionViaEigh The functions below ([`squareroot`](@ref), [`logarithm`](@ref) and [`power`](@ref) with fractional powers) are only defined for matrices whose eigenvalues avoid (part of) the negative real axis, and their principal values are complex whenever eigenvalues on that axis are present. In MatrixAlgebraKit, we aim to keep type stability, and thus the scalar type of the output always matches that of the input. -As such, a real matrix with eigenvalues on the negative real axis leads to a `DomainError`, you should pass a complex matrix instead to obtain the complex principal value. +As such, a real matrix with eigenvalues on the negative real axis leads to a `DomainError`. +You should pass a complex matrix instead to obtain the complex principal value. To avoid spurious errors for eigenvalues that lie on the negative real axis only because of rounding errors (e.g. a positive semidefinite matrix with a tiny negative eigenvalue), eigenvalues within an absolute tolerance `domain_atol` of the domain boundary are clamped onto it. This tolerance defaults to [`default_domain_atol`](@ref) and can be specified explicitly for the algorithms that support it, e.g. `MatrixFunctionViaEigh(eigh_alg; domain_atol=...)`. From dd70ffb612533e563d721ddcaf4a0b7ea87bbd31 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 13:11:10 -0400 Subject: [PATCH 13/22] special case power!(A, 0/1) --- src/implementations/power.jl | 9 ++++++ src/interface/power.jl | 3 ++ test/matrixfunctions/power.jl | 9 ++++++ test/testsuite/matrixfunctions/power.jl | 38 ++++++++++++++++++++++++- 4 files changed, 58 insertions(+), 1 deletion(-) diff --git a/src/implementations/power.jl b/src/implementations/power.jl index fff0ab7be..232fd2201 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -30,11 +30,14 @@ power!(A::AbstractMatrix, p::Real, out, alg::DefaultAlgorithm) = power!(A, p, ou # ------- initialize_output(::typeof(power!), A::AbstractMatrix, p::Real, ::AbstractAlgorithm) = A + # Implementation # -------------- function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaLA) check_input(power!, A, p, powA, alg) isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `power`")) + iszero(p) && return one!(powA) + isone(p) && ((powA === A || copy!(powA, A)); return powA) powAc = A^p if eltype(powAc) <: Complex && !(eltype(powA) <: Complex) # `LinearAlgebra` computes fractional powers of real matrices in complex @@ -53,6 +56,8 @@ end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) check_input(power!, A, p, powA, alg) + iszero(p) && return one!(powA) + isone(p) && ((powA === A || copy!(powA, A)); return powA) D, V = eigh_full!(A, alg.eigh_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if isinteger(p) @@ -68,6 +73,8 @@ end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) check_input(power!, A, p, powA, alg) + iszero(p) && return one!(powA) + isone(p) && ((powA === A || copy!(powA, A)); return powA) D, V = eig_full!(A, alg.eig_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) if eltype(A) <: Real @@ -85,6 +92,8 @@ end # -------------- function power!(A::AbstractMatrix, p::Real, powA, alg::DiagonalAlgorithm) check_input(power!, A, p, powA, alg) + iszero(p) && return one!(powA) + isone(p) && ((powA === A || copy!(powA, A)); return powA) λ = diagview(powA) copyto!(λ, diagview(A)) if isinteger(p) diff --git a/src/interface/power.jl b/src/interface/power.jl index 6da9b42e8..6232b6e22 100644 --- a/src/interface/power.jl +++ b/src/interface/power.jl @@ -11,6 +11,9 @@ Compute the matrix power `powA = A^p` of the square matrix `A`. For integer `p` this is defined for any square matrix (invertible if `p < 0`); for fractional `p` the principal power `exp(p * log(A))` is computed. +The exponents `p = 0` and `p = 1` are resolved directly as `I` and `A`, without computing any +decomposition, so they apply to every square matrix regardless of the algorithm selected. + The scalar type of the output matches that of the input. As a consequence, for fractional `p`, a real matrix with eigenvalues on the negative real axis, for which the principal power is complex, leads to a `DomainError`; pass a diff --git a/test/matrixfunctions/power.jl b/test/matrixfunctions/power.jl index 208a48607..c738ff593 100644 --- a/test/matrixfunctions/power.jl +++ b/test/matrixfunctions/power.jl @@ -35,6 +35,8 @@ if !is_buildkite ) TestSuite.test_power(T, (m, m)) TestSuite.test_power_algs(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_power_trivial(T, (m, m), LAPACK_EIG_ALGS) + TestSuite.test_power_trivial(T, (m, m), LAPACK_EIGH_ALGS; test_rejected_input = true) TestSuite.test_power_hermitian(T, (m, m), LAPACK_EIG_ALGS; exact_hermiticity = false) TestSuite.test_power_hermitian(T, (m, m), LAPACK_EIGH_ALGS) TestSuite.test_power_reference(T, (m, m)) @@ -50,6 +52,8 @@ if !is_buildkite GS_ALGS = (MatrixFunctionViaEig(QRIteration(; driver = GS())),) GLA_ALGS = (MatrixFunctionViaEigh(QRIteration(; driver = GLA())),) TestSuite.test_power_algs(T, (24, 24), GS_ALGS) + TestSuite.test_power_trivial(T, (24, 24), GS_ALGS) + TestSuite.test_power_trivial(T, (24, 24), GLA_ALGS; test_rejected_input = true) TestSuite.test_power_hermitian(T, (24, 24), GS_ALGS; exact_hermiticity = false) TestSuite.test_power_hermitian(T, (24, 24), GLA_ALGS) TestSuite.test_power_domain(T, (md, md), GS_ALGS) @@ -63,6 +67,7 @@ if !is_buildkite test_spectrum = !(T in DiagonalOnlyFloats) TestSuite.test_power(AT, m) TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_trivial(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),); test_spectrum) TestSuite.test_power_reference(AT, m; test_hermitian = !(T in GenericFloats)) TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) @@ -82,11 +87,13 @@ if CUDA.functional() MatrixFunctionViaEigh(DivideAndConquer()), ) TestSuite.test_power_hermitian(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS) + TestSuite.test_power_trivial(CuMatrix{T}, (m, m), CUDA_EIGH_ALGS; test_rejected_input = true) TestSuite.test_power_domain(CuMatrix{T}, (md, md), CUDA_EIGH_ALGS; hermitian_output = true) AT = Diagonal{T, CuVector{T}} TestSuite.test_power(AT, m) TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_trivial(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) end @@ -102,11 +109,13 @@ if AMDGPU.functional() MatrixFunctionViaEigh(DivideAndConquer()), ) TestSuite.test_power_hermitian(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS) + TestSuite.test_power_trivial(ROCMatrix{T}, (m, m), ROC_EIGH_ALGS; test_rejected_input = true) TestSuite.test_power_domain(ROCMatrix{T}, (md, md), ROC_EIGH_ALGS; hermitian_output = true) AT = Diagonal{T, ROCVector{T}} TestSuite.test_power(AT, m) TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) + TestSuite.test_power_trivial(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) end diff --git a/test/testsuite/matrixfunctions/power.jl b/test/testsuite/matrixfunctions/power.jl index 065e7159e..6276f23d2 100644 --- a/test/testsuite/matrixfunctions/power.jl +++ b/test/testsuite/matrixfunctions/power.jl @@ -1,6 +1,6 @@ using TestExtras using LinearAlgebra: LinearAlgebra, I, SingularException -using MatrixAlgebraKit: ishermitian +using MatrixAlgebraKit: ishermitian, one! # `power` takes the exponent as a second positional argument, and splits into an integer branch # (repeated multiplication, defined for any square matrix) and a fractional branch (the principal @@ -71,6 +71,42 @@ function test_power_algs(T::Type, sz, algs; ps = (2, -1), qs = (1 // 2, -1 // 4) end end +# `A^0 = I` and `A^1 = A` hold for every square matrix, and both are short-circuited before any +# decomposition is computed, so the results are exact rather than merely approximate. +# +# `test_rejected_input = true` for algorithms that cannot handle a general matrix at all (the +# `eigh`-based ones): feeding them input they would reject proves the decomposition really is +# skipped, rather than being computed and happening to give the right answer. +function test_power_trivial(T::Type, sz, algs; test_rejected_input = false, kwargs...) + R = real(eltype(T)) + summary_str = testargs_summary(T, sz) + return @testset "power trivial exponents algorithm $alg $summary_str" for alg in algs + A = instantiate_offaxis_matrix(T, sz) + Id = one!(deepcopy(A)) + + @testset "integer p = $p" for p in (0, 1) + powA = @testinferred power(A, p, alg) + @test powA == (iszero(p) ? Id : A) + # in-place, where the output aliases the input + @test power!(deepcopy(A), p, alg) == (iszero(p) ? Id : A) + end + @testset "float p = $p" for p in (zero(R), one(R)) + powA = @testinferred power(A, p, alg) + @test powA == (iszero(p) ? Id : A) + end + + if test_rejected_input + B = instantiate_smallnorm_matrix(T, sz) + IdB = one!(deepcopy(B)) + # confirm the premise: this algorithm cannot decompose `B` + @test_throws DomainError power(B, 2, alg) + # yet the trivial exponents never reach the decomposition + @test power(B, 0, alg) == IdB + @test power(B, 1, alg) == B + end + end +end + # See the corresponding comment in `squareroot.jl` for `exact_hermiticity`. function test_power_hermitian( T::Type, sz, algs; From 3c65bc7b5fd9bc6b0275b6512ef71f6a96878efd Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 13:29:42 -0400 Subject: [PATCH 14/22] refactor: share input handling and reconstruction between matrix functions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `squareroot`, `logarithm` and `power` all map a square matrix to a matrix of the same size and scalar type, so their `copy_input`, `check_input` and `initialize_output` contracts were identical modulo the function name -- six methods repeated three times over. Fold the bodies into `_matrixfunction_copy_input` and `_matrixfunction_check_input`, whose own dispatch absorbs the `Diagonal` and `DiagonalAlgorithm` specializations, so each function keeps just the one line that names it. The three `MatrixFunctionViaEig` kernels shared an identical skeleton, differing only in which scalar function is applied to the eigenvalues: solve `V f(D) V⁻¹` and, for real input, discard the imaginary part left by the complex decomposition. That becomes `_apply_eig!`, taking the already transformed eigenvalues. Similarly `_apply_eigh!` captures the `V f(D) V'`-then-project shape shared by `logarithm` and integer `power`, and pairs with the existing `_mul_herm!` for the cases where a symmetric product makes the result hermitian by construction instead. Deliberately left alone: the `MatrixFunctionViaLA` kernels, whose bodies look alike but must not be merged -- `_copy_result!` did exactly that and was removed in 2698e26, because `squareroot`/`logarithm` need a strict complex-result check while `power` has to tolerate rounding-level imaginary components. `exponential` also stays out, as its real/complex branch keys on the scalar `τ` as well as the matrix eltype. Behaviour is unchanged: the matrix-function tests report the same 4701 passing as before the refactor. Co-Authored-By: Claude Opus 5 (1M context) --- src/implementations/logarithm.jl | 36 ++++--------------- src/implementations/matrixfunctions.jl | 49 ++++++++++++++++++++++++++ src/implementations/power.jl | 43 +++++----------------- src/implementations/squareroot.jl | 31 +++------------- 4 files changed, 69 insertions(+), 90 deletions(-) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index e59140362..315543b81 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -1,24 +1,9 @@ # Inputs # ------ -function copy_input(::typeof(logarithm), A::AbstractMatrix) - return copy!(similar(A, float(eltype(A))), A) -end -copy_input(::typeof(logarithm), A::Diagonal) = map_diagonal(float, A) +copy_input(::typeof(logarithm), A::AbstractMatrix) = _matrixfunction_copy_input(A) function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, alg::AbstractAlgorithm) - m = LinearAlgebra.checksquare(A) - @check_size(logA, (m, m)) - @check_scalar(logA, A) - return nothing -end - -function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, ::DiagonalAlgorithm) - m = LinearAlgebra.checksquare(A) - @assert isdiag(A) - @assert logA isa Diagonal - @check_size(logA, (m, m)) - @check_scalar(logA, A) - return nothing + return _matrixfunction_check_input(A, logA, alg) end # Algorithm selection @@ -48,24 +33,17 @@ end function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh) check_input(logarithm!, A, logA, alg) D, V = eigh_full!(A, alg.eigh_alg) - VD = V * logarithm!(D, D, DiagonalAlgorithm(; domain_atol = alg.domain_atol)) - mul!(logA, VD, V') - return project_hermitian!(logA) + diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) + return _apply_eigh!(logA, V, logarithm!(D, D, diag_alg)) end function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig) check_input(logarithm!, A, logA, alg) D, V = eig_full!(A, alg.eig_alg) + # a real result requires the spectrum to stay off the negative real axis + eltype(A) <: Real && _clamp_domain_eigenvalues!(D, alg.domain_atol) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) - if eltype(A) <: Real - _clamp_domain_eigenvalues!(D, alg.domain_atol) - VlogD = V * logarithm!(D, D, diag_alg) - logAc = rdiv!(VlogD, LinearAlgebra.lu!(V)) - return logA .= real.(logAc) - else - logA .= V .* transpose(diagview(logarithm!(D, D, diag_alg))) - return rdiv!(logA, LinearAlgebra.lu!(V)) - end + return _apply_eig!(logA, V, logarithm!(D, D, diag_alg)) end # Diagonal logic diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index 4e72dacf7..9e6a966d9 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -1,3 +1,52 @@ +# Shared input and output handling +# -------------------------------- +# `squareroot`, `logarithm` and `power` all map a square matrix to a matrix of the same size and +# scalar type, so they agree on `copy_input`, `check_input` and `initialize_output`; only their +# kernels differ. Each implementation file forwards to the helpers below, which keeps the +# per-function definitions to the one line that names the function. + +_matrixfunction_copy_input(A::AbstractMatrix) = copy!(similar(A, float(eltype(A))), A) +_matrixfunction_copy_input(A::Diagonal) = map_diagonal(float, A) + +function _matrixfunction_check_input(A::AbstractMatrix, out, ::AbstractAlgorithm) + m = LinearAlgebra.checksquare(A) + @check_size(out, (m, m)) + @check_scalar(out, A) + return nothing +end + +function _matrixfunction_check_input(A::AbstractMatrix, out, ::DiagonalAlgorithm) + m = LinearAlgebra.checksquare(A) + @assert isdiag(A) + @assert out isa Diagonal + @check_size(out, (m, m)) + @check_scalar(out, A) + return nothing +end + +# Shared reconstruction from an eigenvalue decomposition +# ------------------------------------------------------ +# Both take the already-transformed eigenvalues `fD = f(D)` and rebuild `f(A)`. + +# `f(A) = V f(D) V⁻¹` for a general `A`. A real matrix has a complex decomposition but a real +# `f(A)` whenever `f(D)` closes under conjugation, so the imaginary part is dropped after the +# solve; the callers are responsible for rejecting the inputs where it would not be. +function _apply_eig!(fA, V, fD) + if eltype(fA) <: Real + VfD = V * fD + fAc = rdiv!(VfD, LinearAlgebra.lu!(V)) + return fA .= real.(fAc) + else + fA .= V .* transpose(diagview(fD)) + return rdiv!(fA, LinearAlgebra.lu!(V)) + end +end + +# `f(A) = V f(D) V'` for a hermitian `A`. The product is hermitian only up to roundoff, so it is +# projected afterwards; where `f` admits a square root, prefer building `f(A)` as a symmetric +# product via `_mul_herm!`, which is hermitian by construction. +_apply_eigh!(fA, V, fD) = project_hermitian!(mul!(fA, V * fD, V')) + # Shared helpers for matrix functions with a restricted domain # ------------------------------------------------------------- diff --git a/src/implementations/power.jl b/src/implementations/power.jl index 232fd2201..ce92df8ed 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -1,24 +1,9 @@ # Inputs # ------ -function copy_input(::typeof(power), A::AbstractMatrix, p::Real) - return copy!(similar(A, float(eltype(A))), A), p -end -copy_input(::typeof(power), A::Diagonal, p::Real) = map_diagonal(float, A), p +copy_input(::typeof(power), A::AbstractMatrix, p::Real) = _matrixfunction_copy_input(A), p function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, alg::AbstractAlgorithm) - m = LinearAlgebra.checksquare(A) - @check_size(powA, (m, m)) - @check_scalar(powA, A) - return nothing -end - -function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, ::DiagonalAlgorithm) - m = LinearAlgebra.checksquare(A) - @assert isdiag(A) - @assert powA isa Diagonal - @check_size(powA, (m, m)) - @check_scalar(powA, A) - return nothing + return _matrixfunction_check_input(A, powA, alg) end # Algorithm selection @@ -60,15 +45,9 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) isone(p) && ((powA === A || copy!(powA, A)); return powA) D, V = eigh_full!(A, alg.eigh_alg) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) - if isinteger(p) - VD = V * power!(D, p, D, diag_alg) - mul!(powA, VD, V') - return project_hermitian!(powA) - else - # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction - Vs = rmul!(V, power!(D, p / 2, D, diag_alg)) - return _mul_herm!(powA, Vs) - end + isinteger(p) && return _apply_eigh!(powA, V, power!(D, p, D, diag_alg)) + # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction + return _mul_herm!(powA, rmul!(V, power!(D, p / 2, D, diag_alg))) end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) @@ -76,16 +55,10 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) iszero(p) && return one!(powA) isone(p) && ((powA === A || copy!(powA, A)); return powA) D, V = eig_full!(A, alg.eig_alg) + # only a fractional power of a real matrix needs the spectrum off the negative real axis + eltype(A) <: Real && !isinteger(p) && _clamp_domain_eigenvalues!(D, alg.domain_atol) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) - if eltype(A) <: Real - isinteger(p) || _clamp_domain_eigenvalues!(D, alg.domain_atol) - VpD = V * power!(D, p, D, diag_alg) - powAc = rdiv!(VpD, LinearAlgebra.lu!(V)) - return powA .= real.(powAc) - else - powA .= V .* transpose(diagview(power!(D, p, D, diag_alg))) - return rdiv!(powA, LinearAlgebra.lu!(V)) - end + return _apply_eig!(powA, V, power!(D, p, D, diag_alg)) end # Diagonal logic diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index 67b08c3b3..bcd940251 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -1,24 +1,9 @@ # Inputs # ------ -function copy_input(::typeof(squareroot), A::AbstractMatrix) - return copy!(similar(A, float(eltype(A))), A) -end -copy_input(::typeof(squareroot), A::Diagonal) = map_diagonal(float, A) +copy_input(::typeof(squareroot), A::AbstractMatrix) = _matrixfunction_copy_input(A) function check_input(::typeof(squareroot!), A::AbstractMatrix, sqrtA, alg::AbstractAlgorithm) - m = LinearAlgebra.checksquare(A) - @check_size(sqrtA, (m, m)) - @check_scalar(sqrtA, A) - return nothing -end - -function check_input(::typeof(squareroot!), A::AbstractMatrix, sqrtA, ::DiagonalAlgorithm) - m = LinearAlgebra.checksquare(A) - @assert isdiag(A) - @assert sqrtA isa Diagonal - @check_size(sqrtA, (m, m)) - @check_scalar(sqrtA, A) - return nothing + return _matrixfunction_check_input(A, sqrtA, alg) end # Algorithm selection @@ -56,16 +41,10 @@ end function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig) check_input(squareroot!, A, sqrtA, alg) D, V = eig_full!(A, alg.eig_alg) + # a real result requires the spectrum to stay off the negative real axis + eltype(A) <: Real && _clamp_domain_eigenvalues!(D, alg.domain_atol) diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) - if eltype(A) <: Real - _clamp_domain_eigenvalues!(D, alg.domain_atol) - VsqrtD = V * squareroot!(D, D, diag_alg) - sqrtAc = rdiv!(VsqrtD, LinearAlgebra.lu!(V)) - return sqrtA .= real.(sqrtAc) - else - sqrtA .= V .* transpose(diagview(squareroot!(D, D, diag_alg))) - return rdiv!(sqrtA, LinearAlgebra.lu!(V)) - end + return _apply_eig!(sqrtA, V, squareroot!(D, D, diag_alg)) end # Diagonal logic From 077cc7ba1938b90593772119a13600ec73eddee1 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 15:27:43 -0400 Subject: [PATCH 15/22] fix: make `exponential` hermitian to the last bit for complex input MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `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) --- src/implementations/exponential.jl | 19 ++++++++++++------- test/testsuite/matrixfunctions/exponential.jl | 16 +++++++++++----- 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/implementations/exponential.jl b/src/implementations/exponential.jl index 3e9c92412..4490db40e 100644 --- a/src/implementations/exponential.jl +++ b/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 diff --git a/test/testsuite/matrixfunctions/exponential.jl b/test/testsuite/matrixfunctions/exponential.jl index 1c162094d..93f0ce7d6 100644 --- a/test/testsuite/matrixfunctions/exponential.jl +++ b/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 From 4e3794428c79ebcf877180ec3b0eab4642516e0a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 16:02:40 -0400 Subject: [PATCH 16/22] docs: reference #261 in the changelog and record the behaviour changes Adds the PR link to the existing entries, which every other entry in the file carries, plus the two user-visible changes made since: the trivial exponents of `power` short-circuiting the decomposition, and `exponential` of a complex hermitian matrix now being exactly hermitian. Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/changelog.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 18e579c0b..60bf7f4cf 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -22,9 +22,9 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added -- New matrix functions `squareroot`, `logarithm` and `power` (with both integer and fractional exponents), supporting the `MatrixFunctionViaLA`, `MatrixFunctionViaEig`, `MatrixFunctionViaEigh` and `DiagonalAlgorithm` algorithms. -- The scalar type of the output of these matrix functions matches that of the input, and out-of-domain eigenvalues (e.g. on the negative real axis for a real input) throw a `DomainError`; eigenvalues that violate the domain within a tolerance `domain_atol` (defaulting to `default_domain_atol`) are treated as rounding artifacts and clamped to the domain boundary. -- `MatrixFunctionViaEig` and `MatrixFunctionViaEigh` accept a `domain_atol` keyword argument to control this tolerance. +- New matrix functions `squareroot`, `logarithm` and `power` (with both integer and fractional exponents), supporting the `MatrixFunctionViaLA`, `MatrixFunctionViaEig`, `MatrixFunctionViaEigh` and `DiagonalAlgorithm` algorithms ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). +- The scalar type of the output of these matrix functions matches that of the input, and out-of-domain eigenvalues (e.g. on the negative real axis for a real input) throw a `DomainError`; eigenvalues that violate the domain within a tolerance `domain_atol` (defaulting to `default_domain_atol`) are treated as rounding artifacts and clamped to the domain boundary ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). +- `MatrixFunctionViaEig` and `MatrixFunctionViaEigh` accept a `domain_atol` keyword argument to control this tolerance ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). ### Changed @@ -34,6 +34,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Fixed +- `exponential` of a complex hermitian matrix via `MatrixFunctionViaEigh` is now hermitian to the last bit, rather than only up to rounding errors, matching the guarantee the real case already provided ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). + ### Performance ## [0.6.9](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/compare/v0.6.8...v0.6.9) - 2026-07-10 From a82a76f9c972fef72d820c0baf2fe572e729f5a0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 28 Jul 2026 16:12:09 -0400 Subject: [PATCH 17/22] docs: update the test README for the matrixfunctions group The README added in #262 describes the layout this branch changes: the matrix exponential has moved out of `common`, the `genericschur` and `genericlinearalgebra` groups are gone now that each driver covers the generic element types alongside the BLAS floats, and `matrixfunctions` is a new group. Also records why the `eig`-based algorithms have to name their driver explicitly for those element types, and lists the spectrum-controlled instantiators the matrix functions need. Co-Authored-By: Claude Opus 5 (1M context) --- test/README.md | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/test/README.md b/test/README.md index 18eb02a1d..96cd0c7a2 100644 --- a/test/README.md +++ b/test/README.md @@ -73,13 +73,18 @@ file placed directly in `test/` would belong to no group and never run in CI. | Group | Contents | |-------|----------| -| `common` | Algorithm selection and defaults, truncation strategies, projections, matrix exponential, Aqua code-quality checks | +| `common` | Algorithm selection and defaults, truncation strategies, projections, Aqua code-quality checks | | `decompositions` | `qr`, `lq`, `svd`, `eig`, `eigh`, `gen_eig`, `schur`, `polar`, `orthnull` on CPU and GPU array types | +| `matrixfunctions` | `exponential`, `squareroot`, `logarithm`, `power` on CPU and GPU array types | | `chainrules` | ChainRulesCore rules, exercised through ChainRulesTestUtils and Zygote | | `mooncake` | Mooncake AD rules | | `enzyme` | Enzyme AD rules, exercised through EnzymeTestUtils | -| `genericlinearalgebra` | `MatrixAlgebraKitGenericLinearAlgebraExt` | -| `genericschur` | `MatrixAlgebraKitGenericSchurExt` | + +The generic element types (`BigFloat` and friends) are not a group of their own: each driver covers +them alongside the BLAS floats, naming the `MatrixAlgebraKitGenericSchurExt` and +`MatrixAlgebraKitGenericLinearAlgebraExt` algorithms explicitly. Both extensions are loaded at once, +so a bare `QRIteration()` resolves its driver to GenericLinearAlgebra, which provides no `geev!`; +`eig`-based algorithms therefore have to spell out `QRIteration(; driver = GS())`. Two directories are *not* groups: `testsuite/` holds the shared implementation described below and is excluded both from discovery (in `runtests.jl`) and from CI (`exclude: '["testsuite"]'`). @@ -125,6 +130,11 @@ Supporting infrastructure in the module: - Predicates used throughout the assertions: `isleftnull`, `isrightnull`, `isleftcomplete`, `isrightcomplete`, `has_positive_diagonal`. - `instantiate_unitary`, `instantiate_rank_deficient_matrix` — inputs with prescribed structure. +- `instantiate_smallnorm_matrix`, `instantiate_offaxis_matrix`, `instantiate_posdef_matrix`, + `instantiate_hermitian_spectrum` — inputs with a prescribed *spectrum*, which the matrix functions + need: `squareroot`, `logarithm` and fractional `power` are only defined away from the negative + real axis, and `logarithm` also requires a nonzero spectrum, so the plain `randn` of + `instantiate_matrix` will not do. `test/linearmap.jl` is a further helper, defining a `LinearMap` wrapper that is deliberately *not* an `AbstractMatrix`, used to check the generic code paths. It is excluded from discovery and included From 234d4920c4ef463cc3e5d876d0a3dedd634da0f7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 29 Jul 2026 08:09:31 -0400 Subject: [PATCH 18/22] docs: fix the Pkg workspace link in the test README `pkgdocs.julialang.org/dev/workspaces/` is a 404 -- workspaces have no page of their own, they are a section of the `Project.toml` reference. Point at `v1/toml-files/#Workspaces` instead: the anchor resolves to "The `[workspace]` section", and the versioned channel will not move the way the `dev` docs did. Co-Authored-By: Claude Opus 5 (1M context) --- test/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/README.md b/test/README.md index 96cd0c7a2..e01ea9915 100644 --- a/test/README.md +++ b/test/README.md @@ -4,7 +4,7 @@ Tests are driven by [ParallelTestRunner.jl](https://github.com/JuliaTesting/Para Every `.jl` file under `test/` is auto-discovered and executed in its own worker process, so test files share no state and must be self-contained. -The test environment is a [Pkg workspace](https://pkgdocs.julialang.org/dev/workspaces/) member: +The test environment is a [Pkg workspace](https://pkgdocs.julialang.org/v1/toml-files/#Workspaces) member: `test/Project.toml` declares the test dependencies and resolves the parent package through `[sources]`, while the whole workspace shares a single `Manifest.toml` at the repository root. From c2750a7e25ffff66f5c1a47e19dde1b725025868 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 29 Jul 2026 10:18:40 -0400 Subject: [PATCH 19/22] copyto! -> copy! --- src/implementations/exponential.jl | 4 ++-- src/implementations/logarithm.jl | 2 +- src/implementations/power.jl | 2 +- src/implementations/squareroot.jl | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/implementations/exponential.jl b/src/implementations/exponential.jl index 4490db40e..bd19824cf 100644 --- a/src/implementations/exponential.jl +++ b/src/implementations/exponential.jl @@ -140,7 +140,7 @@ function exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaTaylor) iszero(d[1]) && return one!(expA) powers = Vector{Base.promote_op(similar, typeof(A))}(undef, p₀) - powers[1] = eltype(powers) === typeof(A) ? A : copyto!(similar(A), A) + powers[1] = eltype(powers) === typeof(A) ? A : copy!(similar(A), A) for p in 2:p₀ powers[p] = similar(powers[1]) @@ -179,7 +179,7 @@ function exponential!(A::AbstractMatrix, expA, alg::MatrixFunctionViaTaylor) if dobalance expA .= scale .* X ./ transpose(scale) else - X === expA || copyto!(expA, X) + X === expA || copy!(expA, X) end return expA end diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index 315543b81..44984f323 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -51,7 +51,7 @@ end function logarithm!(A::AbstractMatrix, logA, alg::DiagonalAlgorithm) check_input(logarithm!, A, logA, alg) λ = diagview(logA) - copyto!(λ, diagview(A)) + copy!(λ, diagview(A)) atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) _check_nonzero_eigenvalues(λ, atol) if eltype(λ) <: Real diff --git a/src/implementations/power.jl b/src/implementations/power.jl index ce92df8ed..5b2d9c20c 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -68,7 +68,7 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::DiagonalAlgorithm) iszero(p) && return one!(powA) isone(p) && ((powA === A || copy!(powA, A)); return powA) λ = diagview(powA) - copyto!(λ, diagview(A)) + copy!(λ, diagview(A)) if isinteger(p) p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) else diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index bcd940251..d2e515546 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -52,7 +52,7 @@ end function squareroot!(A::AbstractMatrix, sqrtA, alg::DiagonalAlgorithm) check_input(squareroot!, A, sqrtA, alg) λ = diagview(sqrtA) - copyto!(λ, diagview(A)) + copy!(λ, diagview(A)) if eltype(λ) <: Real atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) _clamp_domain_eigenvalues!(λ, atol) From 4849584018a99e29090f39a4563a83930ed275f0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 29 Jul 2026 10:19:54 -0400 Subject: [PATCH 20/22] code style suggestion --- src/implementations/power.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/implementations/power.jl b/src/implementations/power.jl index 5b2d9c20c..1476e5c2d 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -33,9 +33,9 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaLA) atol = defaulttol(powA) * norm(powAc, Inf) all(x -> abs(imag(x)) <= atol, powAc) || throw(_realness_domainerror(power!)) powA .= real.(powAc) - return powA + else + copy!(powA, powAc) end - copy!(powA, powAc) return powA end From 9f527afce17d1b7143038958856a8f76a8d8a17e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 29 Jul 2026 10:51:12 -0400 Subject: [PATCH 21/22] tolerance in Float64 --- src/implementations/matrixfunctions.jl | 13 +++++-------- src/interface/matrixfunctions.jl | 22 ++++++++++------------ 2 files changed, 15 insertions(+), 20 deletions(-) diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index 9e6a966d9..0d5808aed 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -71,29 +71,26 @@ end ) end +_clamp_domain_eigenvalues!(D::Diagonal, atol::Real) = + _clamp_domain_eigenvalues!(diagview(D), atol) + # Clamp real eigenvalues that are negative within `atol` (rounding artifacts) to zero, # and throw a `DomainError` for eigenvalues that are genuinely negative, since then the # result cannot be expressed with the same (real) scalar type. function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) λmin = minimum(λ; init = zero(eltype(λ))) + atol = atol < 0 ? default_domain_atol(λ) : oftype(λmin, atol) λmin < -atol && throw_negative_eigenvalue(λmin, atol, "a negative real eigenvalue") λ .= max.(λ, zero(eltype(λ))) return λ end -# Convenience method for the eigenvalues of a decomposition, deriving the default -# tolerance from the eigenvalues themselves when `domain_atol` is `nothing`. -function _clamp_domain_eigenvalues!(D::Diagonal, domain_atol::Union{Nothing, Real}) - λ = diagview(D) - atol = something(domain_atol, default_domain_atol(λ)) - return _clamp_domain_eigenvalues!(λ, atol) -end - # Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative # real axis obstruct a real result; complex-conjugate pairs do not. function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real) onaxis = x -> abs(imag(x)) <= atol && real(x) < 0 λmin = mapreduce(x -> onaxis(x) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ)))) + atol = atol < 0 ? default_domain_atol(λ) : oftype(λmin, atol) λmin < -atol && throw_negative_eigenvalue(λmin, atol, "an eigenvalue on the negative real axis") λ .= ifelse.(onaxis.(λ), zero(eltype(λ)), λ) return λ diff --git a/src/interface/matrixfunctions.jl b/src/interface/matrixfunctions.jl index a0e0a7304..39527ba42 100644 --- a/src/interface/matrixfunctions.jl +++ b/src/interface/matrixfunctions.jl @@ -38,17 +38,16 @@ For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`l as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default tolerance [`default_domain_atol`](@ref). """ -struct MatrixFunctionViaEigh{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm +struct MatrixFunctionViaEigh{A <: AbstractAlgorithm} <: AbstractAlgorithm eigh_alg::A - domain_atol::T -end -function MatrixFunctionViaEigh(eigh_alg::AbstractAlgorithm; domain_atol::Union{Nothing, Real} = nothing) - return MatrixFunctionViaEigh(eigh_alg, domain_atol) + domain_atol::Float64 # negative value for runtime defaults end +MatrixFunctionViaEigh(eigh_alg::AbstractAlgorithm; domain_atol::Real = -1.0) = + MatrixFunctionViaEigh(eigh_alg, Float64(domain_atol)) function Base.show(io::IO, alg::MatrixFunctionViaEigh) print(io, "MatrixFunctionViaEigh(") _show_alg(io, alg.eigh_alg) - isnothing(alg.domain_atol) || print(io, "; domain_atol=", alg.domain_atol) + alg.domain_atol < 0 || print(io, "; domain_atol=", alg.domain_atol) return print(io, ")") end @@ -62,16 +61,15 @@ For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`l as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default tolerance [`default_domain_atol`](@ref). """ -struct MatrixFunctionViaEig{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm +struct MatrixFunctionViaEig{A <: AbstractAlgorithm} <: AbstractAlgorithm eig_alg::A - domain_atol::T -end -function MatrixFunctionViaEig(eig_alg::AbstractAlgorithm; domain_atol::Union{Nothing, Real} = nothing) - return MatrixFunctionViaEig(eig_alg, domain_atol) + domain_atol::Float64 # negative value for runtime defaults end +MatrixFunctionViaEig(eig_alg::AbstractAlgorithm; domain_atol::Real = -1.0) = + MatrixFunctionViaEig(eig_alg, Float64(domain_atol)) function Base.show(io::IO, alg::MatrixFunctionViaEig) print(io, "MatrixFunctionViaEig(") _show_alg(io, alg.eig_alg) - isnothing(alg.domain_atol) || print(io, "; domain_atol=", alg.domain_atol) + alg.domain_atol < 0 || print(io, "; domain_atol=", alg.domain_atol) return print(io, ")") end From f94f6bc55e517a8a61f67e23386985f54c6e272b Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 29 Jul 2026 17:56:20 -0400 Subject: [PATCH 22/22] matrix function tolerances are hard --- docs/src/changelog.md | 2 - docs/src/user_interface/algorithms.md | 2 +- docs/src/user_interface/matrix_functions.md | 48 ++++++- ext/MatrixAlgebraKitGenericSchurExt.jl | 2 +- src/common/defaults.jl | 13 +- src/implementations/logarithm.jl | 27 ++-- src/implementations/matrixfunctions.jl | 131 +++++++++++++++---- src/implementations/power.jl | 45 ++++--- src/implementations/squareroot.jl | 27 ++-- src/interface/logarithm.jl | 5 +- src/interface/matrixfunctions.jl | 26 ++-- src/interface/power.jl | 10 +- src/interface/squareroot.jl | 6 +- test/matrixfunctions/power.jl | 6 +- test/testsuite/TestSuite.jl | 12 ++ test/testsuite/matrixfunctions/logarithm.jl | 30 ++++- test/testsuite/matrixfunctions/power.jl | 31 +++-- test/testsuite/matrixfunctions/squareroot.jl | 16 ++- 18 files changed, 323 insertions(+), 116 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 60bf7f4cf..e7c3ea011 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -23,8 +23,6 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Added - New matrix functions `squareroot`, `logarithm` and `power` (with both integer and fractional exponents), supporting the `MatrixFunctionViaLA`, `MatrixFunctionViaEig`, `MatrixFunctionViaEigh` and `DiagonalAlgorithm` algorithms ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). -- The scalar type of the output of these matrix functions matches that of the input, and out-of-domain eigenvalues (e.g. on the negative real axis for a real input) throw a `DomainError`; eigenvalues that violate the domain within a tolerance `domain_atol` (defaulting to `default_domain_atol`) are treated as rounding artifacts and clamped to the domain boundary ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). -- `MatrixFunctionViaEig` and `MatrixFunctionViaEigh` accept a `domain_atol` keyword argument to control this tolerance ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)). ### Changed diff --git a/docs/src/user_interface/algorithms.md b/docs/src/user_interface/algorithms.md index dcb93f9ad..a9677d48f 100644 --- a/docs/src/user_interface/algorithms.md +++ b/docs/src/user_interface/algorithms.md @@ -103,7 +103,7 @@ The following algorithms for matrix functions are available. | Algorithm | Applicable matrix functions | Key keyword arguments | |:----------|:--------------------------|:----------------------| | [`MatrixFunctionViaTaylor`](@ref) | exponential | `tol`, `balance` | -| [`MatrixFunctionViaLA`](@ref) | exponential, squareroot, logarithm, power | | +| [`MatrixFunctionViaLA`](@ref) | exponential, squareroot, logarithm, power | `domain_atol` | | [`MatrixFunctionViaEig`](@ref) | exponential, squareroot, logarithm, power | `eig_alg`, `domain_atol` | | [`MatrixFunctionViaEigh`](@ref) | exponential, squareroot, logarithm, power | `eigh_alg`, `domain_atol` | diff --git a/docs/src/user_interface/matrix_functions.md b/docs/src/user_interface/matrix_functions.md index 0888861cd..439ec54d5 100644 --- a/docs/src/user_interface/matrix_functions.md +++ b/docs/src/user_interface/matrix_functions.md @@ -38,14 +38,56 @@ MatrixAlgebraKit.MatrixFunctionViaEig MatrixAlgebraKit.MatrixFunctionViaEigh ``` -## Domain considerations +## [Domain considerations](@id sec_matrixfunction_domain) The functions below ([`squareroot`](@ref), [`logarithm`](@ref) and [`power`](@ref) with fractional powers) are only defined for matrices whose eigenvalues avoid (part of) the negative real axis, and their principal values are complex whenever eigenvalues on that axis are present. In MatrixAlgebraKit, we aim to keep type stability, and thus the scalar type of the output always matches that of the input. As such, a real matrix with eigenvalues on the negative real axis leads to a `DomainError`. You should pass a complex matrix instead to obtain the complex principal value. -To avoid spurious errors for eigenvalues that lie on the negative real axis only because of rounding errors (e.g. a positive semidefinite matrix with a tiny negative eigenvalue), eigenvalues within an absolute tolerance `domain_atol` of the domain boundary are clamped onto it. -This tolerance defaults to [`default_domain_atol`](@ref) and can be specified explicitly for the algorithms that support it, e.g. `MatrixFunctionViaEigh(eigh_alg; domain_atol=...)`. +To avoid spurious errors for eigenvalues that lie on the negative real axis only because of rounding errors (e.g. a positive semidefinite matrix with a tiny negative eigenvalue), an absolute tolerance `domain_atol` decides which eigenvalues are treated as rounding artifacts. +It can be specified for every algorithm, e.g. `MatrixFunctionViaEigh(eigh_alg; domain_atol=...)` or `squareroot(A; domain_atol=...)`, and defaults to [`default_domain_atol`](@ref). + +### The tolerance has two opposite roles + +Whether the boundary point `λ = 0` belongs to the domain differs per function, and this flips what raising `domain_atol` does. + +| Function | Domain for a real result | `domain_atol` is | Raising it | +|:---------|:-------------------------|:-----------------|:-----------| +| [`squareroot`](@ref) | `λ ∉ ℝ₋`, boundary included | a clamping radius: eigenvalues within `domain_atol` of the negative real axis are moved onto it | accepts more matrices | +| [`power`](@ref), fractional `p > 0` | as `squareroot`, since `0^p = 0` | as `squareroot` | accepts more matrices | +| [`logarithm`](@ref) | `λ ∉ ℝ₋ ∪ {0}`, boundary excluded | a rejection radius: eigenvalues within `domain_atol` of the origin are `DomainError`s, since `log(0)` does not exist | rejects more matrices | +| [`power`](@ref), fractional `p < 0` | as `logarithm`, since `0^p` diverges | as `logarithm` | rejects more matrices | +| [`power`](@ref), integer `p < 0` | invertible matrices | a rejection radius around the origin | rejects more matrices | +| [`power`](@ref), integer `p ≥ 0` | unrestricted | unused | — | +| [`exponential`](@ref) | unrestricted | unused | — | + +So for `logarithm` and negative `power`, raising `domain_atol` never rescues a matrix that was rejected: there is no boundary point to clamp onto, and the tolerance only widens the neighbourhood of the origin that is rejected. +Raising it past a small negative eigenvalue merely turns a "negative eigenvalue" `DomainError` into a "numerically zero eigenvalue" one; *lowering* it is what admits an eigenvalue that is small but genuinely nonzero. + +Clamping is backward stable, but only to the size of the eigenvalue that was discarded, and the forward error it incurs is *not* of that size. +For `squareroot`, clamping an eigenvalue at `-δ` perturbs the result by `O(√δ)`, so an accepted result computed at the default tolerance can differ from the exact principal value by considerably more than the tolerance itself. + +### What the tolerance is measured on + +For the eigenvalue-decomposition-based algorithms, `domain_atol` is the distance of an eigenvalue to the domain boundary, and its default reflects how accurately that algorithm obtains the spectrum. + +| Algorithm | `domain_atol` measures | Default | +|:----------|:-----------------------|:--------| +| [`DiagonalAlgorithm`](@ref) | the eigenvalue itself, which is exact input here | `n * eps * maximum(abs, λ)` | +| [`MatrixFunctionViaEigh`](@ref) | the eigenvalue, accurate to the backward error of a stable hermitian eigensolver | `n * eps * maximum(abs, λ)` | +| [`MatrixFunctionViaEig`](@ref) | the eigenvalue, whose accuracy is additionally limited by the conditioning of the eigenvectors | `defaulttol(λ) * maximum(abs, λ)` | +| [`MatrixFunctionViaLA`](@ref) | the imaginary part of the *result* | `defaulttol * norm(f(A), Inf)` | + +The first two use the same rule as `LinearAlgebra.sqrt(::Hermitian; rtol = eps(T) * size(A, 1))`, so for hermitian input MatrixAlgebraKit and `LinearAlgebra` accept and reject the same matrices. +`MatrixFunctionViaEig` is deliberately looser, since a defective eigenvalue of multiplicity `k` is only resolved to `eps^(1/k)`. + +`MatrixFunctionViaLA` is the exception, because `LinearAlgebra` never exposes the spectrum: it decides internally whether a real result exists and returns a complex matrix when it does not. +The tolerance therefore bounds the imaginary part of the result instead, which is a different quantity on a different scale — for `squareroot`, an eigenvalue at `-δ` shows up as an imaginary part of order `√δ`. +Such a tolerance is not optional there: `LinearAlgebra` casts a fractional power back to a real matrix only when its imaginary part vanishes identically, so `A^p` of a real matrix with complex-conjugate eigenvalues is complex even when the spectrum stays well clear of the negative real axis. + +!!! warning "`MatrixFunctionViaLA` cannot detect a singular matrix" + Because it inspects only the result, `MatrixFunctionViaLA` does not enforce the nonzero-eigenvalue condition of `logarithm` and of `power` with a negative exponent: `LinearAlgebra` treats a numerically zero eigenvalue as in-domain and returns a finite result whose entries are merely large. + Use [`MatrixFunctionViaEig`](@ref) or [`MatrixFunctionViaEigh`](@ref) if you need such input rejected. ```@docs; canonical=false MatrixAlgebraKit.default_domain_atol diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index 93c8f5823..d04ca1a2b 100644 --- a/ext/MatrixAlgebraKitGenericSchurExt.jl +++ b/ext/MatrixAlgebraKitGenericSchurExt.jl @@ -26,7 +26,7 @@ for default_f_algorithm in ( :default_power_algorithm, ) @eval function MatrixAlgebraKit.$default_f_algorithm( - type::Type{T}; domain_atol::Union{Nothing, Real} = nothing, kwargs... + type::Type{T}; domain_atol::Real = -1.0, kwargs... ) where {T <: StridedMatrix{<:GSFloat}} eig_alg = MatrixAlgebraKit.default_eig_algorithm(type; kwargs...) return MatrixFunctionViaEig(eig_alg; domain_atol) diff --git a/src/common/defaults.jl b/src/common/defaults.jl index 008d2ca6b..541b3f0f9 100644 --- a/src/common/defaults.jl +++ b/src/common/defaults.jl @@ -44,13 +44,22 @@ Default tolerance for deciding to warn if the provided `A` is not hermitian. default_hermitian_tol(A) = eps(norm(A, Inf))^(3 / 4) """ - default_domain_atol(λ) + default_domain_atol(λ, alg) Default absolute tolerance for deciding when the eigenvalues `λ` should be considered to lie outside of the domain of a matrix function, e.g. on the negative real axis for [`squareroot`](@ref) and [`logarithm`](@ref) of a real matrix. + +The tolerance has to absorb the error with which `alg` obtained `λ`, and that error differs by +orders of magnitude between the algorithms, so the default is algorithm-dependent. +See [Domain considerations](@ref sec_matrixfunction_domain) for the resulting values. """ -default_domain_atol(λ) = defaulttol(λ) * maximum(abs, λ; init = abs(zero(eltype(λ)))) +function default_domain_atol end + +# roundoff scales only with number of elements and spectrum - e.g. hermitian precision +# conditioning is less strict +_roundoff_domain_atol(λ) = length(λ) * eps(real(float(one(eltype(λ))))) * maximum(abs, λ; init = abs(zero(eltype(λ)))) +_conditioning_domain_atol(λ) = defaulttol(λ) * maximum(abs, λ; init = abs(zero(eltype(λ)))) const DEFAULT_FIXGAUGE = Ref(true) diff --git a/src/implementations/logarithm.jl b/src/implementations/logarithm.jl index 44984f323..be035c924 100644 --- a/src/implementations/logarithm.jl +++ b/src/implementations/logarithm.jl @@ -19,30 +19,34 @@ initialize_output(::typeof(logarithm!), A::AbstractMatrix, ::AbstractAlgorithm) # -------------- function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaLA) check_input(logarithm!, A, logA, alg) - isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `logarithm`")) - # `LinearAlgebra.log` of a real matrix is real whenever the principal logarithm is, - # so a complex result with a real output signals a genuine domain violation + domain_atol = _la_domain_atol(alg, logarithm!) + # `LinearAlgebra.log` of a real matrix is real whenever the principal logarithm is. Note that a + # (numerically) zero eigenvalue goes undetected here, as documented in the manual. logAc = LinearAlgebra.log(A) if eltype(logAc) <: Complex && !(eltype(logA) <: Complex) - throw(_realness_domainerror(logarithm!)) + _la_project_real!(logA, logAc, domain_atol, logarithm!) + else + copy!(logA, logAc) end - copy!(logA, logAc) return logA end function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh) check_input(logarithm!, A, logA, alg) D, V = eigh_full!(A, alg.eigh_alg) - diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) + diag_alg = DiagonalAlgorithm(; domain_atol = _resolve_domain_atol(diagview(D), alg)) return _apply_eigh!(logA, V, logarithm!(D, D, diag_alg)) end function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig) check_input(logarithm!, A, logA, alg) D, V = eig_full!(A, alg.eig_alg) + λ = diagview(D) + atol = _resolve_domain_atol(λ, alg) + _check_nonzero_eigenvalues(λ, atol) # a real result requires the spectrum to stay off the negative real axis - eltype(A) <: Real && _clamp_domain_eigenvalues!(D, alg.domain_atol) - diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) + eltype(A) <: Real && _check_domain_eigenvalues(λ, atol, false) + diag_alg = DiagonalAlgorithm(; domain_atol = atol) return _apply_eig!(logA, V, logarithm!(D, D, diag_alg)) end @@ -52,11 +56,10 @@ function logarithm!(A::AbstractMatrix, logA, alg::DiagonalAlgorithm) check_input(logarithm!, A, logA, alg) λ = diagview(logA) copy!(λ, diagview(A)) - atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) + atol = _resolve_domain_atol(λ, alg) + # `log(0)` does not exist, so the origin is excluded from the domain and nothing is clamped _check_nonzero_eigenvalues(λ, atol) - if eltype(λ) <: Real - _clamp_domain_eigenvalues!(λ, atol) - end + eltype(λ) <: Real && _check_domain_eigenvalues(λ, atol, false) λ .= log.(λ) return logA end diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl index 0d5808aed..631ee3873 100644 --- a/src/implementations/matrixfunctions.jl +++ b/src/implementations/matrixfunctions.jl @@ -49,15 +49,40 @@ _apply_eigh!(fA, V, fD) = project_hermitian!(mul!(fA, V * fD, V')) # Shared helpers for matrix functions with a restricted domain # ------------------------------------------------------------- +# `_clamp_domain_eigenvalues!` is for the functions whose domain includes its boundary and +# `_check_domain_eigenvalues` for the ones that exclude it; see the manual section on domain +# considerations for what that means for `domain_atol`. + +# Each algorithm obtains its eigenvalues with a different accuracy, so each brings its own default. +default_domain_atol(λ, ::DiagonalAlgorithm) = _roundoff_domain_atol(λ) +default_domain_atol(λ, ::MatrixFunctionViaEigh) = _roundoff_domain_atol(λ) +default_domain_atol(λ, ::MatrixFunctionViaEig) = _conditioning_domain_atol(λ) + +# A negative tolerance denotes the runtime default, which keeps the algorithm types concrete. +_domain_atol(alg::Union{MatrixFunctionViaEig, MatrixFunctionViaEigh}) = alg.domain_atol +_domain_atol(alg::DiagonalAlgorithm) = get(alg.kwargs, :domain_atol, -1.0) + +# Callers resolve the default before delegating to an inner `DiagonalAlgorithm`, so that the +# tolerance follows the algorithm that computed the eigenvalues. +function _resolve_domain_atol(λ, alg) + atol = _domain_atol(alg) + R = real(float(eltype(λ))) + return atol < 0 ? convert(R, default_domain_atol(λ, alg)) : convert(R, atol) +end # The throwing branches live in `@noinline` helpers so that the reductions and broadcasts # below stay free of error-path code, which keeps them GPU friendly. -@noinline function throw_negative_eigenvalue(λmin, atol, what) +@noinline function throw_negative_eigenvalue(λmin, atol, what, clampable) + advice = if clampable + "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." + else + "Pass a complex matrix to obtain the principal value; increasing `domain_atol` cannot help, as this matrix function is " * + "undefined on the domain boundary itself." + end return throw( DomainError( λmin, - "The matrix has $what beyond `domain_atol = $atol` and the result of this matrix function is complex. " * - "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact." + "The matrix has $what beyond `domain_atol = $atol` and the result of this matrix function is complex. " * advice ) ) end @@ -66,51 +91,101 @@ end return throw( DomainError( amin, - "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined." + "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined. " * + "Decrease `domain_atol` if the eigenvalue is genuine and should not be treated as zero." ) ) end -_clamp_domain_eigenvalues!(D::Diagonal, atol::Real) = - _clamp_domain_eigenvalues!(diagview(D), atol) - -# Clamp real eigenvalues that are negative within `atol` (rounding artifacts) to zero, -# and throw a `DomainError` for eigenvalues that are genuinely negative, since then the -# result cannot be expressed with the same (real) scalar type. -function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) +# Real eigenvalues that are negative beyond `atol` cannot be expressed with the same (real) scalar +# type. `clampable` only selects the advice in the error message. +function _check_domain_eigenvalues(λ::AbstractVector{<:Real}, atol::Real, clampable::Bool = true) λmin = minimum(λ; init = zero(eltype(λ))) - atol = atol < 0 ? default_domain_atol(λ) : oftype(λmin, atol) - λmin < -atol && throw_negative_eigenvalue(λmin, atol, "a negative real eigenvalue") - λ .= max.(λ, zero(eltype(λ))) + λmin < -atol && throw_negative_eigenvalue(λmin, atol, "a negative real eigenvalue", clampable) return λ end # Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative # real axis obstruct a real result; complex-conjugate pairs do not. +_onaxis(x, atol) = abs(imag(x)) <= atol && real(x) < 0 + +function _check_domain_eigenvalues(λ::AbstractVector{<:Complex}, atol::Real, clampable::Bool = true) + λmin = mapreduce(x -> _onaxis(x, atol) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ)))) + λmin < -atol && throw_negative_eigenvalue(λmin, atol, "an eigenvalue on the negative real axis", clampable) + return λ +end + +# Move the eigenvalues that violate the domain within `atol` onto the boundary. +function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real) + _check_domain_eigenvalues(λ, atol) + λ .= max.(λ, zero(eltype(λ))) + return λ +end + function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real) - onaxis = x -> abs(imag(x)) <= atol && real(x) < 0 - λmin = mapreduce(x -> onaxis(x) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ)))) - atol = atol < 0 ? default_domain_atol(λ) : oftype(λmin, atol) - λmin < -atol && throw_negative_eigenvalue(λmin, atol, "an eigenvalue on the negative real axis") - λ .= ifelse.(onaxis.(λ), zero(eltype(λ)), λ) + _check_domain_eigenvalues(λ, atol) + λ .= ifelse.(_onaxis.(λ, atol), zero(eltype(λ)), λ) return λ end # Reject (numerically) zero eigenvalues for functions that are undefined there, -# e.g. `logarithm` and `power` with a negative fractional power. +# e.g. `logarithm` and `power` with a negative exponent. function _check_nonzero_eigenvalues(λ, atol::Real) amin = minimum(abs, λ; init = typemax(real(eltype(λ)))) amin <= atol && throw_zero_eigenvalue(amin, atol) return λ end -# For `MatrixFunctionViaLA`, domain violations surface as a complex result from -# `LinearAlgebra` while the output should remain real. -function _realness_domainerror(f) - return DomainError( - f, - "The result of this matrix function applied to the given real matrix is complex (eigenvalues on the negative real axis). " * - "Pass a complex matrix to obtain the principal value, or use `MatrixFunctionViaEigh`/`MatrixFunctionViaEig` with a suitable " * - "`domain_atol` if the offending eigenvalues are rounding artifacts." +# Shared helpers for `MatrixFunctionViaLA` +# --------------------------------------- +# `LinearAlgebra` never exposes the spectrum, so the domain check happens in result space and +# `domain_atol` bounds the imaginary part of `f(A)` instead; see the manual for the consequences. + +@noinline function throw_la_kwargs(f, ks) + return throw( + ArgumentError("`MatrixFunctionViaLA` only accepts the `domain_atol` keyword argument for `$f`, got $ks") + ) +end + +# `MatrixFunctionViaLA` accepts generic keywords, so the kernels validate the ones they support. +function _la_domain_atol(alg::MatrixFunctionViaLA, f) + ks = keys(alg.kwargs) + (isempty(ks) || ks == (:domain_atol,)) || throw_la_kwargs(f, ks) + return get(alg.kwargs, :domain_atol, -1.0) +end + +@noinline function throw_complex_result(f, atol, imagmax) + return throw( + DomainError( + f, + "The result of this matrix function applied to the given real matrix is complex (eigenvalues on the negative real axis): " * + "its imaginary part reaches $imagmax, beyond `domain_atol = $atol`. Pass a complex matrix to obtain the principal " * + "value, or increase `domain_atol` if the imaginary part is a rounding artifact." + ) ) end + +@noinline function throw_nonfinite_result(f) + return throw( + DomainError( + f, + "The result of this matrix function is not finite, which signals a (numerically) singular input for which it is undefined. " * + "Use `MatrixFunctionViaEig`/`MatrixFunctionViaEigh` to have the spectrum itself checked against `domain_atol`." + ) + ) +end + +# Project a complex `LinearAlgebra` result onto the real output. A rounding-level imaginary part is +# not a domain violation: `LinearAlgebra` casts back to real only when the imaginary part vanishes +# identically, so `schurpow` yields a complex matrix even for an in-domain real one. +function _la_project_real!(fA, fAc, domain_atol::Real, f) + all(isfinite, fAc) || throw_nonfinite_result(f) + R = real(eltype(fA)) + # the working precision is that of the output: `LinearAlgebra` computes in complex arithmetic + # throughout, so e.g. a `Float32` input promotes all the way to `ComplexF64` + atol = domain_atol < 0 ? defaulttol(fA) * convert(R, norm(fAc, Inf)) : convert(R, domain_atol) + imagmax = convert(R, maximum(abs ∘ imag, fAc; init = zero(real(eltype(fAc))))) + imagmax <= atol || throw_complex_result(f, atol, imagmax) + fA .= real.(fAc) + return fA +end diff --git a/src/implementations/power.jl b/src/implementations/power.jl index 1476e5c2d..a6583a789 100644 --- a/src/implementations/power.jl +++ b/src/implementations/power.jl @@ -20,19 +20,12 @@ initialize_output(::typeof(power!), A::AbstractMatrix, p::Real, ::AbstractAlgori # -------------- function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaLA) check_input(power!, A, p, powA, alg) - isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `power`")) + domain_atol = _la_domain_atol(alg, power!) iszero(p) && return one!(powA) isone(p) && ((powA === A || copy!(powA, A)); return powA) powAc = A^p if eltype(powAc) <: Complex && !(eltype(powA) <: Complex) - # `LinearAlgebra` computes fractional powers of real matrices in complex - # arithmetic and only casts back to real when the result is exactly real, - # so rounding-level imaginary components do not signal a domain violation. - # The tolerance is based on the working precision, which may be lower than - # the result eltype suggests (e.g. `Float32` input promotes to `ComplexF64`). - atol = defaulttol(powA) * norm(powAc, Inf) - all(x -> abs(imag(x)) <= atol, powAc) || throw(_realness_domainerror(power!)) - powA .= real.(powAc) + _la_project_real!(powA, powAc, domain_atol, power!) else copy!(powA, powAc) end @@ -44,10 +37,14 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh) iszero(p) && return one!(powA) isone(p) && ((powA === A || copy!(powA, A)); return powA) D, V = eigh_full!(A, alg.eigh_alg) - diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) - isinteger(p) && return _apply_eigh!(powA, V, power!(D, p, D, diag_alg)) - # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction - return _mul_herm!(powA, rmul!(V, power!(D, p / 2, D, diag_alg))) + diag_alg = DiagonalAlgorithm(; domain_atol = _resolve_domain_atol(diagview(D), alg)) + if isinteger(p) + _apply_eigh!(powA, V, power!(D, p, D, diag_alg)) + else + # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction + _mul_herm!(powA, rmul!(V, power!(D, p / 2, D, diag_alg))) + end + return powA end function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) @@ -55,9 +52,15 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig) iszero(p) && return one!(powA) isone(p) && ((powA === A || copy!(powA, A)); return powA) D, V = eig_full!(A, alg.eig_alg) + λ = diagview(D) + atol = _resolve_domain_atol(λ, alg) + # a negative exponent excludes the origin from the domain, whether or not it is an integer + p < 0 && _check_nonzero_eigenvalues(λ, atol) # only a fractional power of a real matrix needs the spectrum off the negative real axis - eltype(A) <: Real && !isinteger(p) && _clamp_domain_eigenvalues!(D, alg.domain_atol) - diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) + if eltype(A) <: Real && !isinteger(p) + p < 0 ? _check_domain_eigenvalues(λ, atol, false) : _clamp_domain_eigenvalues!(λ, atol) + end + diag_alg = DiagonalAlgorithm(; domain_atol = atol) return _apply_eig!(powA, V, power!(D, p, D, diag_alg)) end @@ -69,12 +72,14 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::DiagonalAlgorithm) isone(p) && ((powA === A || copy!(powA, A)); return powA) λ = diagview(powA) copy!(λ, diagview(A)) - if isinteger(p) - p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0)) - else - atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) + # a nonnegative integer exponent is defined for every square matrix and needs no tolerance + if p < 0 || !isinteger(p) + atol = _resolve_domain_atol(λ, alg) + # a negative exponent excludes the origin from the domain, whether or not it is an integer p < 0 && _check_nonzero_eigenvalues(λ, atol) - eltype(λ) <: Real && _clamp_domain_eigenvalues!(λ, atol) + if eltype(λ) <: Real && !isinteger(p) + p < 0 ? _check_domain_eigenvalues(λ, atol, false) : _clamp_domain_eigenvalues!(λ, atol) + end end λ .= λ .^ p return powA diff --git a/src/implementations/squareroot.jl b/src/implementations/squareroot.jl index d2e515546..490970fc2 100644 --- a/src/implementations/squareroot.jl +++ b/src/implementations/squareroot.jl @@ -19,20 +19,21 @@ initialize_output(::typeof(squareroot!), A::AbstractMatrix, ::AbstractAlgorithm) # -------------- function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaLA) check_input(squareroot!, A, sqrtA, alg) - isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `squareroot`")) - # `LinearAlgebra.sqrt` of a real matrix is real whenever the principal square root is, - # so a complex result with a real output signals a genuine domain violation + domain_atol = _la_domain_atol(alg, squareroot!) + # `LinearAlgebra.sqrt` of a real matrix is real whenever the principal square root is sqrtAc = LinearAlgebra.sqrt(A) - eltype(sqrtAc) <: Complex && !(eltype(sqrtA) <: Complex) && - throw(_realness_domainerror(squareroot!)) - copy!(sqrtA, sqrtAc) + if eltype(sqrtAc) <: Complex && !(eltype(sqrtA) <: Complex) + _la_project_real!(sqrtA, sqrtAc, domain_atol, squareroot!) + else + copy!(sqrtA, sqrtAc) + end return sqrtA end function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEigh) check_input(squareroot!, A, sqrtA, alg) D, V = eigh_full!(A, alg.eigh_alg) - diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) + diag_alg = DiagonalAlgorithm(; domain_atol = _resolve_domain_atol(diagview(D), alg)) # `sqrt(A) = (V * D^(1/4)) * (V * D^(1/4))'` is hermitian by construction Vs = rmul!(V, power!(D, 1 // 4, D, diag_alg)) return _mul_herm!(sqrtA, Vs) @@ -41,9 +42,11 @@ end function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig) check_input(squareroot!, A, sqrtA, alg) D, V = eig_full!(A, alg.eig_alg) + λ = diagview(D) + atol = _resolve_domain_atol(λ, alg) # a real result requires the spectrum to stay off the negative real axis - eltype(A) <: Real && _clamp_domain_eigenvalues!(D, alg.domain_atol) - diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol) + eltype(A) <: Real && _clamp_domain_eigenvalues!(λ, atol) + diag_alg = DiagonalAlgorithm(; domain_atol = atol) return _apply_eig!(sqrtA, V, squareroot!(D, D, diag_alg)) end @@ -53,10 +56,8 @@ function squareroot!(A::AbstractMatrix, sqrtA, alg::DiagonalAlgorithm) check_input(squareroot!, A, sqrtA, alg) λ = diagview(sqrtA) copy!(λ, diagview(A)) - if eltype(λ) <: Real - atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ)) - _clamp_domain_eigenvalues!(λ, atol) - end + # `sqrt(0) = 0`, so the domain includes its boundary + eltype(λ) <: Real && _clamp_domain_eigenvalues!(λ, _resolve_domain_atol(λ, alg)) λ .= sqrt.(λ) return sqrtA end diff --git a/src/interface/logarithm.jl b/src/interface/logarithm.jl index b84499c14..68904e149 100644 --- a/src/interface/logarithm.jl +++ b/src/interface/logarithm.jl @@ -13,8 +13,9 @@ whose eigenvalues have imaginary part in `(-π, π]`. The scalar type of the output matches that of the input. As a consequence, a real matrix with eigenvalues on the negative real axis, for which the principal logarithm is complex, leads to a `DomainError`; pass a complex matrix to obtain the principal value. A matrix with (numerically) zero eigenvalues has no logarithm and also leads to a `DomainError`. -Both checks use a tolerance `domain_atol`, which can be specified for the algorithms -that support it and defaults to [`default_domain_atol`](@ref). +Both checks use a tolerance `domain_atol`, which defaults to [`default_domain_atol`](@ref). As the +origin is excluded from the domain, there is no boundary to clamp onto and raising `domain_atol` +rejects more matrices rather than fewer; see [Domain considerations](@ref sec_matrixfunction_domain). !!! note The bang method `logarithm!` optionally accepts the output structure and possibly destroys the input matrix `A`. diff --git a/src/interface/matrixfunctions.jl b/src/interface/matrixfunctions.jl index 39527ba42..d2e41abdc 100644 --- a/src/interface/matrixfunctions.jl +++ b/src/interface/matrixfunctions.jl @@ -2,9 +2,15 @@ # MATRIX FUNCTION ALGORITHMS # ================================ """ - MatrixFunctionViaLA() + MatrixFunctionViaLA(; domain_atol=-1) Algorithm type to denote computing a function of a matrix `A` via the implementation of `LinearAlgebra`. +For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`logarithm`](@ref)), +`domain_atol` specifies the absolute tolerance on the imaginary part of the result below which a +complex result is attributed to rounding rather than to a domain violation, with a negative value +denoting the default tolerance. Note that this differs from the eigenvalue tolerance of +[`MatrixFunctionViaEig`](@ref) and [`MatrixFunctionViaEigh`](@ref), as `LinearAlgebra` does not +expose the spectrum; see [Domain considerations](@ref sec_matrixfunction_domain). """ @algdef MatrixFunctionViaLA @@ -29,14 +35,15 @@ As this algorithm requires no LAPACK support, it also applies at arbitrary preci @algdef MatrixFunctionViaTaylor """ - MatrixFunctionViaEigh(eigh_alg; domain_atol=nothing) + MatrixFunctionViaEigh(eigh_alg; domain_atol=-1) Algorithm type for computing a function of a matrix by computing its hermitian eigenvalue decomposition and applying the function to the eigenvalues. The `eigh_alg` specifies which hermitian eigendecomposition implementation to use. For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`logarithm`](@ref)), -`domain_atol` specifies the absolute tolerance below which out-of-domain eigenvalues are treated -as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default -tolerance [`default_domain_atol`](@ref). +`domain_atol` specifies the absolute tolerance within which eigenvalues that violate the domain are +treated as rounding artifacts, with a negative value denoting the default tolerance +[`default_domain_atol`](@ref). Depending on the function, such eigenvalues are clamped onto the +domain boundary or rejected; see [Domain considerations](@ref sec_matrixfunction_domain). """ struct MatrixFunctionViaEigh{A <: AbstractAlgorithm} <: AbstractAlgorithm eigh_alg::A @@ -52,14 +59,15 @@ function Base.show(io::IO, alg::MatrixFunctionViaEigh) end """ - MatrixFunctionViaEig(eig_alg; domain_atol=nothing) + MatrixFunctionViaEig(eig_alg; domain_atol=-1) Algorithm type for computing a function of a matrix by computing its eigenvalue decomposition and applying the function to the eigenvalues. The `eig_alg` specifies which eigendecomposition implementation to use. For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`logarithm`](@ref)), -`domain_atol` specifies the absolute tolerance below which out-of-domain eigenvalues are treated -as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default -tolerance [`default_domain_atol`](@ref). +`domain_atol` specifies the absolute tolerance within which eigenvalues that violate the domain are +treated as rounding artifacts, with a negative value denoting the default tolerance +[`default_domain_atol`](@ref). Depending on the function, such eigenvalues are clamped onto the +domain boundary or rejected; see [Domain considerations](@ref sec_matrixfunction_domain). """ struct MatrixFunctionViaEig{A <: AbstractAlgorithm} <: AbstractAlgorithm eig_alg::A diff --git a/src/interface/power.jl b/src/interface/power.jl index 6232b6e22..2844256be 100644 --- a/src/interface/power.jl +++ b/src/interface/power.jl @@ -18,10 +18,12 @@ The scalar type of the output matches that of the input. As a consequence, for fractional `p`, a real matrix with eigenvalues on the negative real axis, for which the principal power is complex, leads to a `DomainError`; pass a complex matrix to obtain the principal value. -Real eigenvalues that are negative within a tolerance `domain_atol` are treated as -rounding artifacts and clamped to zero, where `domain_atol` can be specified for the -algorithms that support it and defaults to [`default_domain_atol`](@ref). -For negative fractional `p`, (numerically) zero eigenvalues also lead to a `DomainError`. +For any `p < 0`, integer or fractional, (numerically) zero eigenvalues also lead to a `DomainError`. + +Both checks use a tolerance `domain_atol`, which defaults to [`default_domain_atol`](@ref). For +`p > 0` it clamps eigenvalues that are negative within the tolerance onto zero, so that raising it +accepts more matrices, whereas for `p < 0` it is a rejection radius around the origin and raising it +rejects more; see [Domain considerations](@ref sec_matrixfunction_domain). !!! note The bang method `power!` optionally accepts the output structure and diff --git a/src/interface/squareroot.jl b/src/interface/squareroot.jl index 295cb4cbe..e7be7d28b 100644 --- a/src/interface/squareroot.jl +++ b/src/interface/squareroot.jl @@ -14,9 +14,9 @@ The scalar type of the output matches that of the input. As a consequence, a real matrix with eigenvalues on the negative real axis, for which the principal square root is complex, leads to a `DomainError`; pass a complex matrix to obtain the principal value. -Real eigenvalues that are negative within a tolerance `domain_atol` are treated as -rounding artifacts and clamped to zero, where `domain_atol` can be specified for the -algorithms that support it and defaults to [`default_domain_atol`](@ref). +Real eigenvalues that are negative within a tolerance `domain_atol` are treated as rounding +artifacts and clamped to zero, so that raising `domain_atol` accepts more matrices. It defaults to +[`default_domain_atol`](@ref); see [Domain considerations](@ref sec_matrixfunction_domain). !!! note The bang method `squareroot!` optionally accepts the output structure and diff --git a/test/matrixfunctions/power.jl b/test/matrixfunctions/power.jl index c738ff593..ba09bc4d9 100644 --- a/test/matrixfunctions/power.jl +++ b/test/matrixfunctions/power.jl @@ -70,7 +70,7 @@ if !is_buildkite TestSuite.test_power_trivial(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),); test_spectrum) TestSuite.test_power_reference(AT, m; test_hermitian = !(T in GenericFloats)) - TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) + TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),)) end end @@ -95,7 +95,7 @@ if CUDA.functional() TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_trivial(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),)) - TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) + TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),)) end end @@ -117,6 +117,6 @@ if AMDGPU.functional() TestSuite.test_power_algs(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_trivial(AT, m, (DiagonalAlgorithm(),)) TestSuite.test_power_hermitian(AT, m, (DiagonalAlgorithm(),)) - TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),); test_singular = true) + TestSuite.test_power_domain(AT, md, (DiagonalAlgorithm(),)) end end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 5dc29ec9a..593d57f73 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -122,6 +122,18 @@ function instantiate_posdef_matrix(T, sz) return project_hermitian!(A * A') + I end +# Rebuild `alg` with an explicit `domain_atol`, so that the domain tests can exercise the tolerance +# without every caller having to spell out the inner decomposition algorithm a second time. +with_domain_atol(alg::MatrixFunctionViaEig, atol) = MatrixFunctionViaEig(alg.eig_alg; domain_atol = atol) +with_domain_atol(alg::MatrixFunctionViaEigh, atol) = MatrixFunctionViaEigh(alg.eigh_alg; domain_atol = atol) +with_domain_atol(::MatrixAlgebraKit.DiagonalAlgorithm, atol) = DiagonalAlgorithm(; domain_atol = atol) +with_domain_atol(::MatrixFunctionViaLA, atol) = MatrixFunctionViaLA(; domain_atol = atol) + +# A tolerance generous enough to admit an eigenvalue at `-√eps`. For `MatrixFunctionViaLA` it bounds +# the imaginary part of the result rather than the spectrum, which sits on a coarser scale. +domain_test_atol(::MatrixAlgebraKit.AbstractAlgorithm, R) = cbrt(eps(R)) +domain_test_atol(::MatrixFunctionViaLA, R) = one(R) / 2 + # Hermitian with the prescribed (real) spectrum `λ`, for probing the domain boundary. function instantiate_hermitian_spectrum(T, sz, λ) A = instantiate_matrix(T, sz) diff --git a/test/testsuite/matrixfunctions/logarithm.jl b/test/testsuite/matrixfunctions/logarithm.jl index e94bc4968..f7736ab27 100644 --- a/test/testsuite/matrixfunctions/logarithm.jl +++ b/test/testsuite/matrixfunctions/logarithm.jl @@ -94,10 +94,34 @@ function test_logarithm_domain( λzero[1] = zero(R) @test_throws DomainError logarithm(instantiate_hermitian_spectrum(T, sz, λzero), alg) - # roundoff-scale negative eigenvalue: clamped onto the boundary, which is still singular + # roundoff-scale negative eigenvalue: too close to the excluded boundary to be admissible λtiny = collect(R, 1:n) - λtiny[1] = -10 * eps(R) - @test_throws DomainError logarithm(instantiate_hermitian_spectrum(T, sz, λtiny), alg) + λtiny[1] = -eps(R) + Atiny = instantiate_hermitian_spectrum(T, sz, λtiny) + @test_throws DomainError logarithm(Atiny, alg) + + # here `domain_atol` is a rejection radius rather than a clamping radius, so raising it + # cannot rescue such a matrix: it only widens the range of eigenvalues that are rejected. + # `λwide` is in domain at the default tolerance yet rejected at the wider one. + wide_alg = with_domain_atol(alg, domain_test_atol(alg, R)) + @test_throws DomainError logarithm(Atiny, wide_alg) + λwide = collect(R, 1:n) + λwide[1] = -sqrt(eps(R)) + Awide = instantiate_hermitian_spectrum(T, sz, λwide) + @test_throws DomainError logarithm(Awide, wide_alg) + + if eltype(T) <: Real || hermitian_output + # at the default tolerance the same eigenvalue is out of domain for being negative + # rather than for being zero, and the reported value is the one the decomposition + # produced rather than a clamped stand-in + @test_throws DomainError logarithm(Awide, alg) + err = try + logarithm(Awide, alg) + catch e + e + end + @test err isa DomainError && !iszero(err.val) + end end end diff --git a/test/testsuite/matrixfunctions/power.jl b/test/testsuite/matrixfunctions/power.jl index 6276f23d2..eab9f86c1 100644 --- a/test/testsuite/matrixfunctions/power.jl +++ b/test/testsuite/matrixfunctions/power.jl @@ -1,5 +1,5 @@ using TestExtras -using LinearAlgebra: LinearAlgebra, I, SingularException +using LinearAlgebra: LinearAlgebra, I using MatrixAlgebraKit: ishermitian, one! # `power` takes the exponent as a second positional argument, and splits into an integer branch @@ -147,13 +147,13 @@ end # `logarithm`, clamping a roundoff-negative eigenvalue onto zero is harmless for a positive # fractional exponent, since `0^p` is well defined there. # -# `hermitian_output = true` and `supports_domain_atol = false`: see `squareroot.jl` and -# `logarithm.jl`. `test_singular = true` only where the zero eigenvalue is exactly representable -# (the `Diagonal` path); the `eig`-based kernels test `any(iszero, λ)` on *computed* eigenvalues, -# which a numerically singular matrix does not satisfy. +# `hermitian_output = true`: see `squareroot.jl`. `supports_domain_atol = false`: see `logarithm.jl`. +# A negative exponent excludes the origin from the domain whether or not it is an integer, so a +# numerically singular matrix is a `DomainError` for `p = -1` just as it is for fractional `p < 0`, +# on every algorithm that inspects the spectrum. function test_power_domain( T::Type, sz, algs; - hermitian_output = false, supports_domain_atol = true, test_singular = false, kwargs... + hermitian_output = false, supports_domain_atol = true, kwargs... ) R = real(eltype(T)) n = sz isa Tuple ? first(sz) : sz @@ -180,18 +180,31 @@ function test_power_domain( # roundoff-scale negative eigenvalue: clamped onto the boundary, which is fine for p > 0 λtiny = collect(R, 1:n) - λtiny[1] = -10 * eps(R) + λtiny[1] = -eps(R) Atiny = instantiate_hermitian_spectrum(T, sz, λtiny) powAtiny = @testinferred power(Atiny, half, alg) @test eltype(powAtiny) == eltype(Atiny) @test powAtiny * powAtiny ≈ Atiny atol = sqrt(eps(R)) - # a negative fractional exponent additionally requires a nonzero spectrum + # a negative exponent additionally requires a nonzero spectrum, integer or not λzero = collect(R, 1:n) λzero[1] = zero(R) Azero = instantiate_hermitian_spectrum(T, sz, λzero) @test_throws DomainError power(Azero, -half, alg) - test_singular && @test_throws SingularException power(Azero, -1, alg) + @test_throws DomainError power(Azero, -1, alg) + + # an explicit `domain_atol` widens the clamping radius for a positive fractional exponent, + # and the rejection radius for a negative one + λwide = collect(R, 1:n) + λwide[1] = -sqrt(eps(R)) + Awide = instantiate_hermitian_spectrum(T, sz, λwide) + wide_alg = with_domain_atol(alg, domain_test_atol(alg, R)) + if eltype(T) <: Real || hermitian_output + @test_throws DomainError power(Awide, half, alg) + end + powAwide = @testinferred power(Awide, half, wide_alg) + @test powAwide * powAwide ≈ Awide atol = sqrt(sqrt(eps(R))) + @test_throws DomainError power(Awide, -half, wide_alg) end end diff --git a/test/testsuite/matrixfunctions/squareroot.jl b/test/testsuite/matrixfunctions/squareroot.jl index 1e2b2b453..16e285d3d 100644 --- a/test/testsuite/matrixfunctions/squareroot.jl +++ b/test/testsuite/matrixfunctions/squareroot.jl @@ -94,11 +94,25 @@ function test_squareroot_domain(T::Type, sz, algs; hermitian_output = false, kwa # roundoff-scale negative eigenvalue: clamped onto the boundary rather than rejected λclamp = collect(R, 1:n) - λclamp[1] = -10 * eps(R) + λclamp[1] = -eps(R) Aclamp = instantiate_hermitian_spectrum(T, sz, λclamp) sqrtAclamp = @testinferred squareroot(Aclamp, alg) @test eltype(sqrtAclamp) == eltype(Aclamp) @test sqrtAclamp * sqrtAclamp ≈ Aclamp atol = sqrt(eps(R)) + + # an eigenvalue beyond every default tolerance is out of domain, while an explicit + # `domain_atol` admits it after all + λwide = collect(R, 1:n) + λwide[1] = -sqrt(eps(R)) + Awide = instantiate_hermitian_spectrum(T, sz, λwide) + if eltype(T) <: Real || hermitian_output + @test_throws DomainError squareroot(Awide, alg) + end + wide_alg = with_domain_atol(alg, domain_test_atol(alg, R)) + sqrtAwide = @testinferred squareroot(Awide, wide_alg) + @test eltype(sqrtAwide) == eltype(Awide) + # accepting is backward stable, but only to the size of the eigenvalue that was discarded + @test sqrtAwide * sqrtAwide ≈ Awide atol = sqrt(sqrt(eps(R))) end end