GitHub

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τA

Compute the exponential of the square matrix A or τ * A,

Note

The bang method exponential! optionally accepts the output structure and possibly destroys the input matrix A. Always use the return value of the function as it may not always be possible to use the provided expA as output.

source
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.
source
MatrixAlgebraKit.MatrixFunctionViaEig — Type
MatrixFunctionViaEig(eig_alg)

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.

source
MatrixAlgebraKit.MatrixFunctionViaEigh — Type
MatrixFunctionViaEigh(eigh_alg)

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.

source