This function computes the (principal) matrix logarithm of a square matrix.
A logarithm of a matrix A is L such that A= e^L
(meaning A == expm(L)), see the documentation for the matrix
exponential, expm, which can be defined
as
sum[r = 0,..,Inf; L^r/r!] .
Usage
logm(x, method = c("Higham08", "Eigen"),
tol = .Machine$double.eps)
Arguments
x
a square matrix.
method
a string specifying the algorithmic method to be used.
The default uses the algorithm by Higham(2008).
The simple "Eigen" method tries to diagonalise the matrix
x; if that is not possible, it raises an error.
tol
a given tolerance used to check if x is
computationally singular when method = "Eigen".
Details
The exponential of a matrix is defined as the infinite Taylor series
exp(M) = I + M + M^2/2! + M^3/3! + ….
The matrix logarithm of A is a matrix M such that
exp(M) = A. Note that there typically are an infinite number
number of such matrices, and we compute the prinicipal matrix
logarithm, see the references.
Method "Higham08" works via “inverse scaling and
squaring”, and from the Schur decomposition, applying a matrix
square root computation. It is somewhat slow but also works for
non-diagonalizable matrices.
Value
A matrix ‘as x’ with the matrix logarithm of x,
i.e., all.equal( expm(logm(x)), x, tol) is typically true for
quite small tolerance tol.
Author(s)
Method "Higham08" was implemented by Michael Stadelmann as part of his
master thesis in mathematics, at ETH Zurich;
the "Eigen" method by Christophe Dutang.
References
Higham, N.~J. (2008).
Functions of Matrices: Theory and Computation;
Society for Industrial and Applied Mathematics, Philadelphia, PA, USA.
m <- diag(2)
logm(m)
expm(logm(m))
## Here, logm() is barely defined, and Higham08 has needed an amendment
## in order for not to loop forever:
D0 <- diag(x=c(1, 0.))
(L. <- logm(D0))
stopifnot( all.equal(D0, expm(L.)) )
## A matrix for which clearly no logm(.) exists:
(m <- cbind(1:2, 1))
(l.m <- logm(m)) ## all NA
stopifnot(is.na(l.m))
## The "Eigen" method ``works'' but wrongly :
expm(logm(m, "Eigen"))