Skip to content

Commit bbd0c06

Browse files
lkdvoslkdvos
committed
feat: add squareroot, logarithm and power matrix functions
Following the `exponential` template, each function supports `MatrixFunctionViaLA` (wrapping the LinearAlgebra implementation), `MatrixFunctionViaEig`/`MatrixFunctionViaEigh` (eigendecomposition-based, also serving generic eltypes via GenericSchur/GenericLinearAlgebra), and `DiagonalAlgorithm`. `power(A, p)` uses the two-argument `@functiondef` form with an exponent-type-based split between integer and fractional powers. The output scalar type always matches the input: out-of-domain eigenvalues (negative real axis for real input, zero eigenvalues for `logarithm` and negative fractional powers) throw a `DomainError`, while eigenvalues that violate the domain within a tolerance `domain_atol` (new keyword of `MatrixFunctionViaEig(h)`, defaulting to the new `default_domain_atol`) are clamped as rounding artifacts.
1 parent eb1e927 commit bbd0c06

11 files changed

Lines changed: 545 additions & 6 deletions

File tree

‎ext/MatrixAlgebraKitGenericSchurExt.jl‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,18 @@ function MatrixAlgebraKit.default_exponential_algorithm(
2121
return MatrixFunctionViaEig(eig_alg)
2222
end
2323

24+
for default_f_algorithm in (
25+
:default_squareroot_algorithm, :default_logarithm_algorithm,
26+
:default_power_algorithm,
27+
)
28+
@eval function MatrixAlgebraKit.$default_f_algorithm(
29+
type::Type{T}; domain_atol::Union{Nothing, Real} = nothing, kwargs...
30+
) where {T <: StridedMatrix{<:GSFloat}}
31+
eig_alg = MatrixAlgebraKit.default_eig_algorithm(type; kwargs...)
32+
return MatrixFunctionViaEig(eig_alg; domain_atol)
33+
end
34+
end
35+
2436
function geev!(::GS, A::AbstractMatrix, Dd::AbstractVector, V::AbstractMatrix; kwargs...)
2537
D, Vmat = GenericSchur.eigen!(A)
2638
copyto!(Dd, D)

‎src/MatrixAlgebraKit.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,9 @@ export left_polar!, right_polar!
3131
export left_orth, right_orth, left_null, right_null
3232
export left_orth!, right_orth!, left_null!, right_null!
3333
export exponential, exponential!
34+
export squareroot, squareroot!
35+
export logarithm, logarithm!
36+
export power, power!
3437

3538
export Householder, Native_HouseholderQR, Native_HouseholderLQ
3639
export DivideAndConquer, SafeDivideAndConquer, QRIteration, Bisection, Jacobi, SVDViaPolar
@@ -115,6 +118,9 @@ include("interface/schur.jl")
115118
include("interface/polar.jl")
116119
include("interface/orthnull.jl")
117120
include("interface/exponential.jl")
121+
include("interface/squareroot.jl")
122+
include("interface/logarithm.jl")
123+
include("interface/power.jl")
118124

119125
include("implementations/projections.jl")
120126
include("implementations/truncation.jl")
@@ -128,6 +134,10 @@ include("implementations/schur.jl")
128134
include("implementations/polar.jl")
129135
include("implementations/orthnull.jl")
130136
include("implementations/exponential.jl")
137+
include("implementations/matrixfunctions.jl")
138+
include("implementations/squareroot.jl")
139+
include("implementations/logarithm.jl")
140+
include("implementations/power.jl")
131141

132142
include("common/gauge.jl") # needs to be defined after the functions are
133143

‎src/common/defaults.jl‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,15 @@ Default tolerance for deciding to warn if the provided `A` is not hermitian.
4343
"""
4444
default_hermitian_tol(A) = eps(norm(A, Inf))^(3 / 4)
4545

46+
"""
47+
default_domain_atol(λ)
48+
49+
Default absolute tolerance for deciding when the eigenvalues `λ` should be considered
50+
to lie outside of the domain of a matrix function, e.g. on the negative real axis for
51+
[`squareroot`](@ref) and [`logarithm`](@ref) of a real matrix.
52+
"""
53+
default_domain_atol(λ) = defaulttol(λ) * maximum(abs, λ; init = abs(zero(eltype(λ))))
54+
4655

4756
const DEFAULT_FIXGAUGE = Ref(true)
4857

‎src/implementations/logarithm.jl‎

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
# Inputs
2+
# ------
3+
function copy_input(::typeof(logarithm), A::AbstractMatrix)
4+
return copy!(similar(A, float(eltype(A))), A)
5+
end
6+
copy_input(::typeof(logarithm), A::Diagonal) = Diagonal(float.(diagview(A)))
7+
8+
function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, alg::AbstractAlgorithm)
9+
m = LinearAlgebra.checksquare(A)
10+
@check_size(logA, (m, m))
11+
@check_scalar(logA, A)
12+
return nothing
13+
end
14+
15+
function check_input(::typeof(logarithm!), A::AbstractMatrix, logA, ::DiagonalAlgorithm)
16+
m = LinearAlgebra.checksquare(A)
17+
@assert isdiag(A)
18+
@assert logA isa Diagonal
19+
@check_size(logA, (m, m))
20+
@check_scalar(logA, A)
21+
return nothing
22+
end
23+
24+
# Algorithm selection
25+
# -------------------
26+
logarithm!(A::AbstractMatrix, alg::DefaultAlgorithm) = logarithm!(A, select_algorithm(logarithm!, A, nothing; alg.kwargs...))
27+
logarithm!(A::AbstractMatrix, out, alg::DefaultAlgorithm) = logarithm!(A, out, select_algorithm(logarithm!, A, nothing; alg.kwargs...))
28+
29+
# Outputs
30+
# -------
31+
initialize_output(::typeof(logarithm!), A::AbstractMatrix, ::AbstractAlgorithm) = A
32+
33+
# Implementation
34+
# --------------
35+
function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaLA)
36+
check_input(logarithm!, A, logA, alg)
37+
isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `logarithm`"))
38+
result = LinearAlgebra.log(A)
39+
_copy_result!(logarithm!, logA, result)
40+
return logA
41+
end
42+
43+
function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh)
44+
check_input(logarithm!, A, logA, alg)
45+
D, V = eigh_full!(A, alg.eigh_alg)
46+
λ = diagview(D)
47+
atol = something(alg.domain_atol, default_domain_atol(λ))
48+
_check_nonzero_eigenvalues(λ, atol)
49+
_clamp_domain_eigenvalues!(λ, atol)
50+
λ .= log.(λ)
51+
VD = V * D
52+
mul!(logA, VD, V')
53+
return project_hermitian!(logA)
54+
end
55+
56+
function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig)
57+
check_input(logarithm!, A, logA, alg)
58+
D, V = eig_full!(A, alg.eig_alg)
59+
λ = diagview(D)
60+
atol = something(alg.domain_atol, default_domain_atol(λ))
61+
_check_nonzero_eigenvalues(λ, atol)
62+
if eltype(A) <: Real
63+
_clamp_domain_eigenvalues!(λ, atol)
64+
λ .= log.(λ)
65+
VD = V * D
66+
logAc = rdiv!(VD, LinearAlgebra.lu!(V))
67+
return logA .= real.(logAc)
68+
else
69+
λ .= log.(λ)
70+
logA .= V .* transpose(λ)
71+
return rdiv!(logA, LinearAlgebra.lu!(V))
72+
end
73+
end
74+
75+
# Diagonal logic
76+
# --------------
77+
function logarithm!(A::AbstractMatrix, logA, alg::DiagonalAlgorithm)
78+
check_input(logarithm!, A, logA, alg)
79+
λ = diagview(logA)
80+
copyto!(λ, diagview(A))
81+
atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ))
82+
_check_nonzero_eigenvalues(λ, atol)
83+
if eltype(λ) <: Real
84+
_clamp_domain_eigenvalues!(λ, atol)
85+
end
86+
λ .= log.(λ)
87+
return logA
88+
end

‎src/implementations/matrixfunctions.jl‎

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
# Shared helpers for matrix functions with a restricted domain
2+
# -------------------------------------------------------------
3+
4+
# Clamp real eigenvalues that are negative within `atol` (rounding artifacts) to zero,
5+
# and throw a `DomainError` for eigenvalues that are genuinely negative, since then the
6+
# result cannot be expressed with the same (real) scalar type.
7+
function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real)
8+
for i in eachindex(λ)
9+
x = λ[i]
10+
if x < -atol
11+
throw(
12+
DomainError(
13+
x,
14+
"The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " *
15+
"Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact."
16+
)
17+
)
18+
elseif x < 0
19+
λ[i] = zero(x)
20+
end
21+
end
22+
return λ
23+
end
24+
25+
# Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative
26+
# real axis obstruct a real result; complex-conjugate pairs do not.
27+
function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real)
28+
for i in eachindex(λ)
29+
x = λ[i]
30+
if abs(imag(x)) <= atol && real(x) < 0
31+
if real(x) < -atol
32+
throw(
33+
DomainError(
34+
x,
35+
"The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " *
36+
"Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact."
37+
)
38+
)
39+
else
40+
λ[i] = zero(x)
41+
end
42+
end
43+
end
44+
return λ
45+
end
46+
47+
# Reject (numerically) zero eigenvalues for functions that are undefined there,
48+
# e.g. `logarithm` and `power` with a negative fractional power.
49+
function _check_nonzero_eigenvalues(λ, atol::Real)
50+
for x in λ
51+
if abs(x) <= atol
52+
throw(
53+
DomainError(
54+
x,
55+
"The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined."
56+
)
57+
)
58+
end
59+
end
60+
return λ
61+
end
62+
63+
# For `MatrixFunctionViaLA`, domain violations surface as a complex result from
64+
# `LinearAlgebra`; detect this after the fact to preserve type-stable output.
65+
# `LinearAlgebra` sometimes computes a real result in complex arithmetic (e.g. fractional
66+
# powers via a complex Schur decomposition), so rounding-level imaginary components do
67+
# not signal a domain violation and are simply discarded.
68+
function _copy_result!(f, out::AbstractMatrix, result::AbstractMatrix)
69+
if eltype(result) <: Complex && !(eltype(out) <: Complex)
70+
# base the tolerance on the working precision, which is the lower of the two:
71+
# `LinearAlgebra` may internally promote, e.g. `Float32` fractional powers are
72+
# computed with `Float64`-eltype results but `Float32`-level accuracy
73+
atol = max(defaulttol(out), defaulttol(result)) * norm(result, Inf)
74+
for x in result
75+
if abs(imag(x)) > atol
76+
throw(
77+
DomainError(
78+
f,
79+
"The result of this matrix function applied to the given real matrix is complex (eigenvalues on the negative real axis). " *
80+
"Pass a complex matrix to obtain the principal value, or use `MatrixFunctionViaEigh`/`MatrixFunctionViaEig` with a suitable " *
81+
"`domain_atol` if the offending eigenvalues are rounding artifacts."
82+
)
83+
)
84+
end
85+
end
86+
out .= real.(result)
87+
else
88+
copy!(out, result)
89+
end
90+
return out
91+
end

‎src/implementations/power.jl‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
# Inputs
2+
# ------
3+
function copy_input(::typeof(power), A::AbstractMatrix, p::Real)
4+
return copy!(similar(A, float(eltype(A))), A), p
5+
end
6+
copy_input(::typeof(power), A::Diagonal, p::Real) = Diagonal(float.(diagview(A))), p
7+
8+
function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, alg::AbstractAlgorithm)
9+
m = LinearAlgebra.checksquare(A)
10+
@check_size(powA, (m, m))
11+
@check_scalar(powA, A)
12+
return nothing
13+
end
14+
15+
function check_input(::typeof(power!), A::AbstractMatrix, p::Real, powA, ::DiagonalAlgorithm)
16+
m = LinearAlgebra.checksquare(A)
17+
@assert isdiag(A)
18+
@assert powA isa Diagonal
19+
@check_size(powA, (m, m))
20+
@check_scalar(powA, A)
21+
return nothing
22+
end
23+
24+
# Algorithm selection
25+
# -------------------
26+
power!(A::AbstractMatrix, p::Real, alg::DefaultAlgorithm) = power!(A, p, select_algorithm(power!, (A, p), nothing; alg.kwargs...))
27+
power!(A::AbstractMatrix, p::Real, out, alg::DefaultAlgorithm) = power!(A, p, out, select_algorithm(power!, (A, p), nothing; alg.kwargs...))
28+
29+
# Outputs
30+
# -------
31+
initialize_output(::typeof(power!), A::AbstractMatrix, p::Real, ::AbstractAlgorithm) = A
32+
33+
# Implementation
34+
# --------------
35+
function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaLA)
36+
check_input(power!, A, p, powA, alg)
37+
isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `power`"))
38+
result = A^p
39+
_copy_result!(power!, powA, result)
40+
return powA
41+
end
42+
43+
function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh)
44+
check_input(power!, A, p, powA, alg)
45+
D, V = eigh_full!(A, alg.eigh_alg)
46+
λ = diagview(D)
47+
if isinteger(p)
48+
p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0))
49+
λ .= λ .^ p
50+
VD = V * D
51+
mul!(powA, VD, V')
52+
else
53+
atol = something(alg.domain_atol, default_domain_atol(λ))
54+
p < 0 && _check_nonzero_eigenvalues(λ, atol)
55+
_clamp_domain_eigenvalues!(λ, atol)
56+
# `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction
57+
λ .= λ .^ (p / 2)
58+
Vs = rmul!(V, D)
59+
mul!(powA, Vs, Vs')
60+
end
61+
return project_hermitian!(powA)
62+
end
63+
64+
function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig)
65+
check_input(power!, A, p, powA, alg)
66+
D, V = eig_full!(A, alg.eig_alg)
67+
λ = diagview(D)
68+
if isinteger(p)
69+
p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0))
70+
else
71+
atol = something(alg.domain_atol, default_domain_atol(λ))
72+
p < 0 && _check_nonzero_eigenvalues(λ, atol)
73+
eltype(A) <: Real && _clamp_domain_eigenvalues!(λ, atol)
74+
end
75+
if eltype(A) <: Real
76+
λ .= λ .^ p
77+
VD = V * D
78+
powAc = rdiv!(VD, LinearAlgebra.lu!(V))
79+
return powA .= real.(powAc)
80+
else
81+
λ .= λ .^ p
82+
powA .= V .* transpose(λ)
83+
return rdiv!(powA, LinearAlgebra.lu!(V))
84+
end
85+
end
86+
87+
# Diagonal logic
88+
# --------------
89+
function power!(A::AbstractMatrix, p::Real, powA, alg::DiagonalAlgorithm)
90+
check_input(power!, A, p, powA, alg)
91+
λ = diagview(powA)
92+
copyto!(λ, diagview(A))
93+
if isinteger(p)
94+
p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0))
95+
else
96+
atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ))
97+
p < 0 && _check_nonzero_eigenvalues(λ, atol)
98+
eltype(λ) <: Real && _clamp_domain_eigenvalues!(λ, atol)
99+
end
100+
λ .= λ .^ p
101+
return powA
102+
end

0 commit comments

Comments
 (0)