logm {expm} | R Documentation |
Matrix Logarithm
Description
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
e^L := \sum_{r=0}^\infty 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 |
tol |
a given tolerance used to check if |
Details
The exponential of a matrix is defined as the infinite Taylor series
e^M = \sum_{k = 1}^\infty \frac{M^k}{k!}.
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.
The Matrix Logarithm is very nicely defined by Wikipedia, https://en.wikipedia.org/wiki/Matrix_logarithm.
See Also
Examples
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 <- try(logm(m))) ## all NA {Warning in sqrt(S[ij, ij]) : NaNs produced}
## on r-patched-solaris-x86, additionally gives
## Error in solve.default(X[ii, ii] + X[ij, ij], S[ii, ij] - sumU) :
## system is computationally singular: reciprocal condition number = 0
## Calls: logm ... logm.Higham08 -> rootS -> solve -> solve -> solve.default
if(!inherits(l.m, "try-error")) stopifnot(is.na(l.m))
## The "Eigen" method ``works'' but wrongly :
expm(logm(m, "Eigen"))