logm {expm}

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

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

expm

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"))