CurrentModule = MatrixAlgebraKit
CollapsedDocStrings = true
Another class of matrix algebra methods consists of calculating some function of a single input A.
In order to streamline these functions, they all follow a similar common code pattern.
For a given function f, this consists of the following methods:
f(A; kwargs...) -> F...
f!(A, [F]; kwargs...) -> F...Here, the input matrix is always the first argument, and optionally the output can be provided as well.
The keywords are algorithm-specific, and can be used to influence the behavior of the algorithms.
For a full description of how to select and configure algorithms, see [Algorithm Selection](@ref sec_algorithmselection).
Importantly, for generic code patterns it is recommended to always use the output F explicitly, since some implementations may not be able to reuse the provided memory.
Additionally, the f! method typically assumes that it is allowed to destroy the input A, and making use of the contents of A afterwards should be deemed as undefined behavior.
The exponential of a square matrix A is used in many scientific applications, as it arises in the solution of an autonomous linear differential equation.
The default algorithm MatrixFunctionViaTaylor is a pure-Julia scaling-and-squaring evaluation of the Taylor series.
As it requires no LAPACK support, it also applies to generic data types at arbitrary precision.
Alternatively, an implementation based on a Padé approximation is available in LinearAlgebra, and can be accessed by the algorithm MatrixFunctionViaLA.
The exponential can also be calculated by first calculating the (hermitian) eigenvalue decomposition, and then computing the scalar exponential of the diagonal elements.
This strategy is implemented via the algorithms MatrixFunctionViaEig and MatrixFunctionViaEigh, and call eig_full and eigh_full, respectively.
Additionally, in order to calculate exp(τ * A), the function exponential can be called with (τ, A), using the same algorithms as before.
exponential
MatrixAlgebraKit.MatrixFunctionViaTaylor
MatrixAlgebraKit.MatrixFunctionViaLA
MatrixAlgebraKit.MatrixFunctionViaEig
MatrixAlgebraKit.MatrixFunctionViaEigh
The functions below (squareroot, logarithm and power 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.
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 |
λ ∉ ℝ₋, boundary included |
a clamping radius: eigenvalues within domain_atol of the negative real axis are moved onto it |
accepts more matrices |
power, fractional p > 0 |
as squareroot, since 0^p = 0 |
as squareroot |
accepts more matrices |
logarithm |
λ ∉ ℝ₋ ∪ {0}, boundary excluded |
a rejection radius: eigenvalues within domain_atol of the origin are DomainErrors, since log(0) does not exist |
rejects more matrices |
power, fractional p < 0 |
as logarithm, since 0^p diverges |
as logarithm |
rejects more matrices |
power, integer p < 0 |
invertible matrices | a rejection radius around the origin | rejects more matrices |
power, integer p ≥ 0 |
unrestricted | unused | — |
exponential |
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.
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 |
the eigenvalue itself, which is exact input here | n * eps * maximum(abs, λ) |
MatrixFunctionViaEigh |
the eigenvalue, accurate to the backward error of a stable hermitian eigensolver | n * eps * maximum(abs, λ) |
MatrixFunctionViaEig |
the eigenvalue, whose accuracy is additionally limited by the conditioning of the eigenvectors | defaulttol(λ) * maximum(abs, λ) |
MatrixFunctionViaLA |
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 or MatrixFunctionViaEigh if you need such input rejected.
MatrixAlgebraKit.default_domain_atol
The principal square root of a square matrix A is the unique square root whose eigenvalues have nonnegative real part.
It is computed by the function squareroot, where the default algorithm MatrixFunctionViaLA wraps the Schur-based implementation of LinearAlgebra, and the eigenvalue-decomposition-based algorithms MatrixFunctionViaEig and MatrixFunctionViaEigh are available as well.
squareroot
The principal logarithm 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, with the same algorithm choices as squareroot.
logarithm
Matrix powers A^p for real p are computed by the function power, 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.
power