Matrix Logarithm

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! .$$

logm(x, method = c("Higham08", "Eigen"),
<!-- %      order = 8, trySym = TRUE, -->
     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.

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.

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.

See also

expm

Author

Method "Higham08" was implemented by Michael Stadelmann as part of his master thesis in mathematics, at ETH Zurich; the "Eigen" method by Christophe Dutang.

Examples

m <- diag(2)
logm(m)
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
expm(logm(m))
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1

## 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))
#> Warning: Inverse scaling did not work (t = 1).
#> The matrix logarithm may not exist for this matrix.Setting m = 3 arbitrarily.
#>      [,1]     [,2]
#> [1,]    0        0
#> [2,]    0 -9576994
stopifnot( all.equal(D0, expm(L.)) )

## A matrix for which clearly no logm(.) exists:
(m <- cbind(1:2, 1))
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    2    1
(l.m <- try(logm(m))) ## all NA {Warning in sqrt(S[ij, ij]) : NaNs produced}
#> Warning: NaNs produced
#> Warning: NA/NaN from  || Tr - I ||  after 1 step.
#> The matrix logarithm may not exist for this matrix.
#>      [,1] [,2]
#> [1,]  NaN  NaN
#> [2,]  NaN  NaN
## 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"))
#>          [,1]      [,2]
#> [1,] 1.414214 0.7071068
#> [2,] 1.414214 1.4142136