feat: squareroot, logarithm and power

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

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).
 

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
+```

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)

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
 

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)
 

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

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