diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 58f704b6e..e7c3ea011 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -22,6 +22,8 @@ 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)). + ### Changed ### Deprecated @@ -30,6 +32,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 diff --git a/docs/src/user_interface/algorithms.md b/docs/src/user_interface/algorithms.md index d90b76074..a9677d48f 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 | `domain_atol` | +| [`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..439ec54d5 100644 --- a/docs/src/user_interface/matrix_functions.md +++ b/docs/src/user_interface/matrix_functions.md @@ -37,3 +37,85 @@ MatrixAlgebraKit.MatrixFunctionViaLA MatrixAlgebraKit.MatrixFunctionViaEig MatrixAlgebraKit.MatrixFunctionViaEigh ``` + +## [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), 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 +``` + +## 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 +``` diff --git a/ext/MatrixAlgebraKitGenericSchurExt.jl b/ext/MatrixAlgebraKitGenericSchurExt.jl index c34f1abc4..d04ca1a2b 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::Real = -1.0, 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..541b3f0f9 100644 --- a/src/common/defaults.jl +++ b/src/common/defaults.jl @@ -43,6 +43,24 @@ 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(λ, 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. +""" +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/exponential.jl b/src/implementations/exponential.jl index 3e9c92412..bd19824cf 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 @@ -135,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]) @@ -174,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 new file mode 100644 index 000000000..be035c924 --- /dev/null +++ b/src/implementations/logarithm.jl @@ -0,0 +1,65 @@ +# Inputs +# ------ +copy_input(::typeof(logarithm), A::AbstractMatrix) = _matrixfunction_copy_input(A) + +function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, alg::AbstractAlgorithm) + return _matrixfunction_check_input(A, logA, alg) +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) + 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) + _la_project_real!(logA, logAc, domain_atol, logarithm!) + else + copy!(logA, logAc) + end + 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 = _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 && _check_domain_eigenvalues(λ, atol, false) + diag_alg = DiagonalAlgorithm(; domain_atol = atol) + return _apply_eig!(logA, V, logarithm!(D, D, diag_alg)) +end + +# Diagonal logic +# -------------- +function logarithm!(A::AbstractMatrix, logA, alg::DiagonalAlgorithm) + check_input(logarithm!, A, logA, alg) + λ = diagview(logA) + copy!(λ, diagview(A)) + 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) + eltype(λ) <: Real && _check_domain_eigenvalues(λ, atol, false) + λ .= log.(λ) + return logA +end diff --git a/src/implementations/matrixfunctions.jl b/src/implementations/matrixfunctions.jl new file mode 100644 index 000000000..631ee3873 --- /dev/null +++ b/src/implementations/matrixfunctions.jl @@ -0,0 +1,191 @@ +# 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 +# ------------------------------------------------------------- +# `_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, 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. " * advice + ) + ) +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. " * + "Decrease `domain_atol` if the eigenvalue is genuine and should not be treated as zero." + ) + ) +end + +# 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(λ))) + λ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) + _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 exponent. +function _check_nonzero_eigenvalues(λ, atol::Real) + amin = minimum(abs, λ; init = typemax(real(eltype(λ)))) + amin <= atol && throw_zero_eigenvalue(amin, atol) + return λ +end + +# 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/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 new file mode 100644 index 000000000..a6583a789 --- /dev/null +++ b/src/implementations/power.jl @@ -0,0 +1,86 @@ +# Inputs +# ------ +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) + return _matrixfunction_check_input(A, powA, alg) +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) + 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) + _la_project_real!(powA, powAc, domain_atol, power!) + else + copy!(powA, powAc) + end + return powA +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 = _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) + 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) + λ = 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 + 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 + +# Diagonal logic +# -------------- +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) + copy!(λ, diagview(A)) + # 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) + if eltype(λ) <: Real && !isinteger(p) + p < 0 ? _check_domain_eigenvalues(λ, atol, false) : _clamp_domain_eigenvalues!(λ, atol) + end + end + λ .= λ .^ p + return powA +end 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 new file mode 100644 index 000000000..490970fc2 --- /dev/null +++ b/src/implementations/squareroot.jl @@ -0,0 +1,63 @@ +# Inputs +# ------ +copy_input(::typeof(squareroot), A::AbstractMatrix) = _matrixfunction_copy_input(A) + +function check_input(::typeof(squareroot!), A::AbstractMatrix, sqrtA, alg::AbstractAlgorithm) + return _matrixfunction_check_input(A, sqrtA, alg) +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) + 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) + 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 = _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) +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!(λ, atol) + diag_alg = DiagonalAlgorithm(; domain_atol = atol) + return _apply_eig!(sqrtA, V, squareroot!(D, D, diag_alg)) +end + +# Diagonal logic +# -------------- +function squareroot!(A::AbstractMatrix, sqrtA, alg::DiagonalAlgorithm) + check_input(squareroot!, A, sqrtA, alg) + λ = diagview(sqrtA) + copy!(λ, diagview(A)) + # `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 new file mode 100644 index 000000000..68904e149 --- /dev/null +++ b/src/interface/logarithm.jl @@ -0,0 +1,38 @@ +# 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 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`. + 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..d2e41abdc 100644 --- a/src/interface/matrixfunctions.jl +++ b/src/interface/matrixfunctions.jl @@ -1,10 +1,16 @@ # ================================ -# EXPONENTIAL ALGORITHMS +# MATRIX FUNCTION ALGORITHMS # ================================ """ - MatrixFunctionViaLA() + MatrixFunctionViaLA(; domain_atol=-1) -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`. +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,31 +35,49 @@ As this algorithm requires no LAPACK support, it also applies at arbitrary preci @algdef MatrixFunctionViaTaylor """ - MatrixFunctionViaEigh(eigh_alg) + 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 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 + 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) + alg.domain_atol < 0 || print(io, "; domain_atol=", alg.domain_atol) return print(io, ")") end """ - MatrixFunctionViaEig(eig_alg) + 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 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 + 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) + alg.domain_atol < 0 || 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..2844256be --- /dev/null +++ b/src/interface/power.jl @@ -0,0 +1,47 @@ +# 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 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 +complex matrix to obtain the principal value. +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 + 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..e7be7d28b --- /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, 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 + 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 diff --git a/test/README.md b/test/README.md index 18eb02a1d..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. @@ -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 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/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/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..ba09bc4d9 --- /dev/null +++ b/test/matrixfunctions/power.jl @@ -0,0 +1,122 @@ +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_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)) + 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_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) + 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_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(),)) + 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_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(),)) + 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_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(),)) + 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/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 36ae68304..593d57f73 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -102,6 +102,52 @@ 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 + +# 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) + 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 +163,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..93f0ce7d6 --- /dev/null +++ b/test/testsuite/matrixfunctions/exponential.jl @@ -0,0 +1,178 @@ +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 + + # `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 + +# `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..f7736ab27 --- /dev/null +++ b/test/testsuite/matrixfunctions/logarithm.jl @@ -0,0 +1,145 @@ +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: too close to the excluded boundary to be admissible + λtiny = collect(R, 1:n) + λ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 + +# 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..eab9f86c1 --- /dev/null +++ b/test/testsuite/matrixfunctions/power.jl @@ -0,0 +1,233 @@ +using TestExtras +using LinearAlgebra: LinearAlgebra, I +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 +# 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 + +# `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; + 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`: 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, 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] = -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 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_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 + +# 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..16e285d3d --- /dev/null +++ b/test/testsuite/matrixfunctions/squareroot.jl @@ -0,0 +1,133 @@ +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] = -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 + +# 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