matrix function tolerances are hard

Author: lkdvos

docs/src/changelog.md

@@ -23,8 +23,6 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 ### Added
 
 - New matrix functions `squareroot`, `logarithm` and `power` (with both integer and fractional exponents), supporting the `MatrixFunctionViaLA`, `MatrixFunctionViaEig`, `MatrixFunctionViaEigh` and `DiagonalAlgorithm` algorithms ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)).
-- The scalar type of the output of these matrix functions matches that of the input, and out-of-domain eigenvalues (e.g. on the negative real axis for a real input) throw a `DomainError`; eigenvalues that violate the domain within a tolerance `domain_atol` (defaulting to `default_domain_atol`) are treated as rounding artifacts and clamped to the domain boundary ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)).
-- `MatrixFunctionViaEig` and `MatrixFunctionViaEigh` accept a `domain_atol` keyword argument to control this tolerance ([#261](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/261)).
 
 ### Changed
 

docs/src/user_interface/algorithms.md

@@ -103,7 +103,7 @@ The following algorithms for matrix functions are available.
 | Algorithm | Applicable matrix functions | Key keyword arguments |
 |:----------|:--------------------------|:----------------------|
 | [`MatrixFunctionViaTaylor`](@ref) | exponential | `tol`, `balance` |
-| [`MatrixFunctionViaLA`](@ref) | exponential, squareroot, logarithm, power | |
+| [`MatrixFunctionViaLA`](@ref) | exponential, squareroot, logarithm, power | `domain_atol` |
 | [`MatrixFunctionViaEig`](@ref) | exponential, squareroot, logarithm, power | `eig_alg`, `domain_atol` |
 | [`MatrixFunctionViaEigh`](@ref) | exponential, squareroot, logarithm, power | `eigh_alg`, `domain_atol` |
 

docs/src/user_interface/matrix_functions.md

@@ -38,14 +38,56 @@ MatrixAlgebraKit.MatrixFunctionViaEig
 MatrixAlgebraKit.MatrixFunctionViaEigh
 ```
 
-## Domain considerations
+## [Domain considerations](@id sec_matrixfunction_domain)
 
 The functions below ([`squareroot`](@ref), [`logarithm`](@ref) and [`power`](@ref) with fractional powers) are only defined for matrices whose eigenvalues avoid (part of) the negative real axis, and their principal values are complex whenever eigenvalues on that axis are present.
 In MatrixAlgebraKit, we aim to keep type stability, and thus the scalar type of the output always matches that of the input.
 As such, a real matrix with eigenvalues on the negative real axis leads to a `DomainError`.
 You should pass a complex matrix instead to obtain the complex principal value.
-To avoid spurious errors for eigenvalues that lie on the negative real axis only because of rounding errors (e.g. a positive semidefinite matrix with a tiny negative eigenvalue), eigenvalues within an absolute tolerance `domain_atol` of the domain boundary are clamped onto it.
-This tolerance defaults to [`default_domain_atol`](@ref) and can be specified explicitly for the algorithms that support it, e.g. `MatrixFunctionViaEigh(eigh_alg; domain_atol=...)`.
+To avoid spurious errors for eigenvalues that lie on the negative real axis only because of rounding errors (e.g. a positive semidefinite matrix with a tiny negative eigenvalue), an absolute tolerance `domain_atol` decides which eigenvalues are treated as rounding artifacts.
+It can be specified for every algorithm, e.g. `MatrixFunctionViaEigh(eigh_alg; domain_atol=...)` or `squareroot(A; domain_atol=...)`, and defaults to [`default_domain_atol`](@ref).
+
+### The tolerance has two opposite roles
+
+Whether the boundary point `λ = 0` belongs to the domain differs per function, and this flips what raising `domain_atol` does.
+
+| Function | Domain for a real result | `domain_atol` is | Raising it |
+|:---------|:-------------------------|:-----------------|:-----------|
+| [`squareroot`](@ref) | `λ ∉ ℝ₋`, boundary included | a clamping radius: eigenvalues within `domain_atol` of the negative real axis are moved onto it | accepts more matrices |
+| [`power`](@ref), fractional `p > 0` | as `squareroot`, since `0^p = 0` | as `squareroot` | accepts more matrices |
+| [`logarithm`](@ref) | `λ ∉ ℝ₋ ∪ {0}`, boundary excluded | a rejection radius: eigenvalues within `domain_atol` of the origin are `DomainError`s, since `log(0)` does not exist | rejects more matrices |
+| [`power`](@ref), fractional `p < 0` | as `logarithm`, since `0^p` diverges | as `logarithm` | rejects more matrices |
+| [`power`](@ref), integer `p < 0` | invertible matrices | a rejection radius around the origin | rejects more matrices |
+| [`power`](@ref), integer `p ≥ 0` | unrestricted | unused | — |
+| [`exponential`](@ref) | unrestricted | unused | — |
+
+So for `logarithm` and negative `power`, raising `domain_atol` never rescues a matrix that was rejected: there is no boundary point to clamp onto, and the tolerance only widens the neighbourhood of the origin that is rejected.
+Raising it past a small negative eigenvalue merely turns a "negative eigenvalue" `DomainError` into a "numerically zero eigenvalue" one; *lowering* it is what admits an eigenvalue that is small but genuinely nonzero.
+
+Clamping is backward stable, but only to the size of the eigenvalue that was discarded, and the forward error it incurs is *not* of that size.
+For `squareroot`, clamping an eigenvalue at `-δ` perturbs the result by `O(√δ)`, so an accepted result computed at the default tolerance can differ from the exact principal value by considerably more than the tolerance itself.
+
+### What the tolerance is measured on
+
+For the eigenvalue-decomposition-based algorithms, `domain_atol` is the distance of an eigenvalue to the domain boundary, and its default reflects how accurately that algorithm obtains the spectrum.
+
+| Algorithm | `domain_atol` measures | Default |
+|:----------|:-----------------------|:--------|
+| [`DiagonalAlgorithm`](@ref) | the eigenvalue itself, which is exact input here | `n * eps * maximum(abs, λ)` |
+| [`MatrixFunctionViaEigh`](@ref) | the eigenvalue, accurate to the backward error of a stable hermitian eigensolver | `n * eps * maximum(abs, λ)` |
+| [`MatrixFunctionViaEig`](@ref) | the eigenvalue, whose accuracy is additionally limited by the conditioning of the eigenvectors | `defaulttol(λ) * maximum(abs, λ)` |
+| [`MatrixFunctionViaLA`](@ref) | the imaginary part of the *result* | `defaulttol * norm(f(A), Inf)` |
+
+The first two use the same rule as `LinearAlgebra.sqrt(::Hermitian; rtol = eps(T) * size(A, 1))`, so for hermitian input MatrixAlgebraKit and `LinearAlgebra` accept and reject the same matrices.
+`MatrixFunctionViaEig` is deliberately looser, since a defective eigenvalue of multiplicity `k` is only resolved to `eps^(1/k)`.
+
+`MatrixFunctionViaLA` is the exception, because `LinearAlgebra` never exposes the spectrum: it decides internally whether a real result exists and returns a complex matrix when it does not.
+The tolerance therefore bounds the imaginary part of the result instead, which is a different quantity on a different scale — for `squareroot`, an eigenvalue at `-δ` shows up as an imaginary part of order `√δ`.
+Such a tolerance is not optional there: `LinearAlgebra` casts a fractional power back to a real matrix only when its imaginary part vanishes identically, so `A^p` of a real matrix with complex-conjugate eigenvalues is complex even when the spectrum stays well clear of the negative real axis.
+
+!!! warning "`MatrixFunctionViaLA` cannot detect a singular matrix"
+    Because it inspects only the result, `MatrixFunctionViaLA` does not enforce the nonzero-eigenvalue condition of `logarithm` and of `power` with a negative exponent: `LinearAlgebra` treats a numerically zero eigenvalue as in-domain and returns a finite result whose entries are merely large.
+    Use [`MatrixFunctionViaEig`](@ref) or [`MatrixFunctionViaEigh`](@ref) if you need such input rejected.
 
 ```@docs; canonical=false
 MatrixAlgebraKit.default_domain_atol

ext/MatrixAlgebraKitGenericSchurExt.jl

@@ -26,7 +26,7 @@ for default_f_algorithm in (
         :default_power_algorithm,
     )
     @eval function MatrixAlgebraKit.$default_f_algorithm(
-            type::Type{T}; domain_atol::Union{Nothing, Real} = nothing, kwargs...
+            type::Type{T}; domain_atol::Real = -1.0, kwargs...
         ) where {T <: StridedMatrix{<:GSFloat}}
         eig_alg = MatrixAlgebraKit.default_eig_algorithm(type; kwargs...)
         return MatrixFunctionViaEig(eig_alg; domain_atol)

src/common/defaults.jl

@@ -44,13 +44,22 @@ Default tolerance for deciding to warn if the provided `A` is not hermitian.
 default_hermitian_tol(A) = eps(norm(A, Inf))^(3 / 4)
 
 """
-    default_domain_atol(λ)
+    default_domain_atol(λ, alg)
 
 Default absolute tolerance for deciding when the eigenvalues `λ` should be considered
 to lie outside of the domain of a matrix function, e.g. on the negative real axis for
 [`squareroot`](@ref) and [`logarithm`](@ref) of a real matrix.
+
+The tolerance has to absorb the error with which `alg` obtained `λ`, and that error differs by
+orders of magnitude between the algorithms, so the default is algorithm-dependent.
+See [Domain considerations](@ref sec_matrixfunction_domain) for the resulting values.
 """
-default_domain_atol(λ) = defaulttol(λ) * maximum(abs, λ; init = abs(zero(eltype(λ))))
+function default_domain_atol end
+
+# roundoff scales only with number of elements and spectrum - e.g. hermitian precision
+# conditioning is less strict
+_roundoff_domain_atol(λ) = length(λ) * eps(real(float(one(eltype(λ))))) * maximum(abs, λ; init = abs(zero(eltype(λ))))
+_conditioning_domain_atol(λ) = defaulttol(λ) * maximum(abs, λ; init = abs(zero(eltype(λ))))
 
 
 const DEFAULT_FIXGAUGE = Ref(true)

src/implementations/logarithm.jl

@@ -19,30 +19,34 @@ initialize_output(::typeof(logarithm!), A::AbstractMatrix, ::AbstractAlgorithm)
 # --------------
 function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaLA)
     check_input(logarithm!, A, logA, alg)
-    isempty(alg.kwargs) || throw(ArgumentError("`MatrixFunctionViaLA` does not accept keyword arguments for `logarithm`"))
-    # `LinearAlgebra.log` of a real matrix is real whenever the principal logarithm is,
-    # so a complex result with a real output signals a genuine domain violation
+    domain_atol = _la_domain_atol(alg, logarithm!)
+    # `LinearAlgebra.log` of a real matrix is real whenever the principal logarithm is. Note that a
+    # (numerically) zero eigenvalue goes undetected here, as documented in the manual.
     logAc = LinearAlgebra.log(A)
     if eltype(logAc) <: Complex && !(eltype(logA) <: Complex)
-        throw(_realness_domainerror(logarithm!))
+        _la_project_real!(logA, logAc, domain_atol, logarithm!)
+    else
+        copy!(logA, logAc)
     end
-    copy!(logA, logAc)
     return logA
 end
 
 function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh)
     check_input(logarithm!, A, logA, alg)
     D, V = eigh_full!(A, alg.eigh_alg)
-    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
+    diag_alg = DiagonalAlgorithm(; domain_atol = _resolve_domain_atol(diagview(D), alg))
     return _apply_eigh!(logA, V, logarithm!(D, D, diag_alg))
 end
 
 function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig)
     check_input(logarithm!, A, logA, alg)
     D, V = eig_full!(A, alg.eig_alg)
+    λ = diagview(D)
+    atol = _resolve_domain_atol(λ, alg)
+    _check_nonzero_eigenvalues(λ, atol)
     # a real result requires the spectrum to stay off the negative real axis
-    eltype(A) <: Real && _clamp_domain_eigenvalues!(D, alg.domain_atol)
-    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
+    eltype(A) <: Real && _check_domain_eigenvalues(λ, atol, false)
+    diag_alg = DiagonalAlgorithm(; domain_atol = atol)
     return _apply_eig!(logA, V, logarithm!(D, D, diag_alg))
 end
 
@@ -52,11 +56,10 @@ function logarithm!(A::AbstractMatrix, logA, alg::DiagonalAlgorithm)
     check_input(logarithm!, A, logA, alg)
     λ = diagview(logA)
     copy!(λ, diagview(A))
-    atol = something(get(alg.kwargs, :domain_atol, nothing), default_domain_atol(λ))
+    atol = _resolve_domain_atol(λ, alg)
+    # `log(0)` does not exist, so the origin is excluded from the domain and nothing is clamped
     _check_nonzero_eigenvalues(λ, atol)
-    if eltype(λ) <: Real
-        _clamp_domain_eigenvalues!(λ, atol)
-    end
+    eltype(λ) <: Real && _check_domain_eigenvalues(λ, atol, false)
     λ .= log.(λ)
     return logA
 end