Matrix functions
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. 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.
Exponential
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.
MatrixAlgebraKit.exponential — Function
exponential(A; kwargs...) -> expA
exponential(A, alg::AbstractAlgorithm) -> expA
exponential!(A, [expA]; kwargs...) -> expA
exponential!(A, [expA], alg::AbstractAlgorithm) -> expA
exponential((τ, A); kwargs...) -> expτA
exponential((τ, A), alg::AbstractAlgorithm) -> expτA
exponential!((τ, A), [expA]; kwargs...) -> expτA
exponential!((τ, A), [expA], alg::AbstractAlgorithm) -> expτACompute the exponential of the square matrix A or τ * A,
MatrixAlgebraKit.MatrixFunctionViaTaylor — Type
MatrixFunctionViaTaylor(; tol=eps, balance=true, estimate_order=4)Algorithm type to denote finding the exponential of A through a pure-Julia scaling-and-squaring evaluation of its Taylor series, following Fasi & Higham (2018). The truncation order and the number of squarings are chosen to reach a relative accuracy tol, and the Taylor polynomial is evaluated with the Paterson–Stockmeyer scheme. When balance is true, A is first balanced by a diagonal similarity. estimate_order sets how many powers of A are formed up front to sharpen the norm estimate via the Al-Mohy–Higham quantities ‖Aᵖ‖^(1/p) (Al-Mohy & Higham, 2009); these powers are reused by the Paterson–Stockmeyer evaluation. As this algorithm requires no LAPACK support, it also applies at arbitrary precision.
References
- A. H. Al-Mohy and N. J. Higham, "A New Scaling and Squaring Algorithm for the Matrix Exponential", SIAM J. Matrix Anal. Appl., 31(3), 970–989, 2009.
MatrixAlgebraKit.MatrixFunctionViaLA — Type
MatrixFunctionViaLA(; domain_atol=-1)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 and logarithm), 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 and MatrixFunctionViaEigh, as LinearAlgebra does not expose the spectrum; see Domain considerations.
MatrixAlgebraKit.MatrixFunctionViaEig — Type
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 and logarithm), 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. Depending on the function, such eigenvalues are clamped onto the domain boundary or rejected; see Domain considerations.
MatrixAlgebraKit.MatrixFunctionViaEigh — Type
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 and logarithm), 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. Depending on the function, such eigenvalues are clamped onto the domain boundary or rejected; see Domain considerations.
Domain considerations
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.
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 | λ ∉ ℝ₋, 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.
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 | 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.
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 — Function
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 and logarithm 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 for the resulting values.
Square root
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.
MatrixAlgebraKit.squareroot — Function
squareroot(A; kwargs...) -> sqrtA
squareroot(A, alg::AbstractAlgorithm) -> sqrtA
squareroot!(A, [sqrtA]; kwargs...) -> sqrtA
squareroot!(A, [sqrtA], alg::AbstractAlgorithm) -> sqrtACompute 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; see Domain considerations.
Logarithm
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.
MatrixAlgebraKit.logarithm — Function
logarithm(A; kwargs...) -> logA
logarithm(A, alg::AbstractAlgorithm) -> logA
logarithm!(A, [logA]; kwargs...) -> logA
logarithm!(A, [logA], alg::AbstractAlgorithm) -> logACompute 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. 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.
Power
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.
MatrixAlgebraKit.power — Function
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) -> powACompute 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. 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.