expm (version 0.999-4)

logm: 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 "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".

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.

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.

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, http://en.wikipedia.org/wiki/Matrix_logarithm.

See Also

expm

Examples

# NOT RUN {
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 or even error {on Solaris}
# }
# NOT RUN {
## The "Eigen" method  ``works''  but wrongly :
expm(logm(m, "Eigen"))
# }