[R] implementation of matrix logarithm (inverse of matrix exponential)

Charles C. Berry cberry at tajo.ucsd.edu
Mon Sep 28 04:59:52 CEST 2009


On Sun, 27 Sep 2009, Douglas Bates wrote:

> There is a logm function in the expm package in the expm project on
> R-forge.  See http://expm.R-forge.R-project.org/
>
> Martin was the person who added that function so I will defer to his
> explanations of what it does.  I know he has been traveling and it may
> be a day or two before he can get to this thread.
>

Douglas,

Thanks for this note.

The logm() function in the expm package does the trick.

It seems to properly handle block structures for non-diagonalizable 
matrices:

> require(expm)
> M <- matrix(c(0,1,0,0),2)
> expm(M)
      [,1] [,2]
[1,]    1    0
[2,]    1    1
> logm(expm(M))
      [,1] [,2]
[1,]    0    0
[2,]    1    0

And thanks to Vincent, Martin, et al for this package!

:-)

Chuck


> On Sun, Sep 27, 2009 at 11:44 AM, Charles C. Berry <cberry at tajo.ucsd.edu> wrote:
>> On Sat, 26 Sep 2009, spencerg wrote:
>>
>>>    Sylvester's formula
>>> (http://en.wikipedia.org/wiki/Sylvester%27s_formula) applies to a square
>>> matrix A = S L solve(S), where L = a diagonal matrix and S = matrix of
>>> eigenvectors.  Let "f" be an analytic function [for which f(A) is well
>>> defined].  Then f(A) = S f(L) solve(S).
>>>    We can code this as follows:
>>> sylvester <- function(x, f){
>>>  n <- nrow(x)
>>>  eig <- eigen(x)
>>>  vi <- solve(eig$vectors)
>>>  with(eig, (vectors * rep(f(values), each=n)) %*% vi)
>>> }
>>>
>>>
>>> logm <- function(x)sylvester(x, log)
>>>
>>>
>>> Example:
>>> A <- matrix(1:4, 2)
>>> eA <- expm(A)
>>> logm(eA)
>>>
>>>
>>>          With Chuck Berry's example, we get the following:
>>> M <- matrix( c(0,1,0,0), 2 )
>>> sylvester(M, log)
>>
>>
>> The case I gave would be
>>
>>        sylvester( as.matrix( expm( M ) ), log )
>>
>> for which the perfectly sensible answer is M,  not what appears here:
>>
>>
>>> Error in solve.default(eig$vectors) :
>>>  system is computationally singular: reciprocal condition number =
>>>  1.00208e-292
>>>
>>>
>>>          This is a perfectly sensible answer in this case.  We get the
>>> same result from "sylvester(M, exp)", though "expm(M)" works fine.
>>>          A better algorithm for this could be obtains by studying the code
>>> for "expm" in the Matrix package and the references in the associated help
>>> page.
>>
>> I doubt that code already in R will handle cases requiring Jordan blocks for
>> evaluation of the matrix logarithm (which cases arise in the context of
>> discrete state, continuous time Markov chains) without requiring one to
>> built that code more or less from scratch.
>>
>> I'd be happy to hear that this is not so.
>>
>> HTH,
>>
>> Chuck
>>
>>>
>>>     Hope this helps. Spencer
>>>
>>>
>>> Gabor Grothendieck wrote:
>>>>
>>>>  Often one uses matrix logarithms on symmetric positive definite
>>>>  matrices so the assumption of being symmetric is sufficient in many
>>>>  cases.
>>>>
>>>>  On Sat, Sep 26, 2009 at 7:28 PM, Charles C. Berry <cberry at tajo.ucsd.edu>
>>>>  wrote:
>>>>
>>>>>  On Sat, 26 Sep 2009, Gabor Grothendieck wrote:
>>>>>>>>  OK. Try this:
>>>>>>>>>>>  library(Matrix)
>>>>>>>  M <- matrix(c(2, 1, 1, 2), 2); M
>>>>>>>>>     [,1] [,2]
>>>>>>  [1,]    2    1
>>>>>>  [2,]    1    2
>>>>>>>>>  Right. expm( M ) is diagonalizable.
>>>>>>  But for
>>>>>>  M <- matrix( c(0,1,0,0), 2 )
>>>>>>  you get the wrong result.
>>>>>>  Maybe I should have added that I do not see the machinery in R for >
>>>>>>  dealing
>>>>>  with Jordan blocks.
>>>>>>  HTH,
>>>>>>  Chuck
>>>>>>>>>>>  # log of expm(M) is original matrix M
>>>>>>>  with(eigen(expm(M)), vectors %*% diag(log(values)) %*% t(vectors))
>>>>>>>>>     [,1] [,2]
>>>>>>  [1,]    2    1
>>>>>>  [2,]    1    2
>>>>>>>>>>  On Sat, Sep 26, 2009 at 6:24 PM, Charles C. Berry > >
>>>>>>>>>>  <cberry at tajo.ucsd.edu>
>>>>>>  wrote:
>>>>>>>>>  On Sat, 26 Sep 2009, Gabor Grothendieck wrote:
>>>>>>>>>>>>>>  Try:
>>>>>>>>>>>>  expm( - M)
>>>>>>>>>>>  Mimosa probably meant say 'the inverse function'.
>>>>>>>>>>  I do not see one in R.
>>>>>>>>>>  Chuck
>>>>>>>>>>>>>>  On Sat, Sep 26, 2009 at 5:06 PM, Mimosa Zeus
>>>>>>>>>>>>>> <mimosa1879 at yahoo.fr>
>>>>>>>>  wrote:
>>>>>>>>>>>>>  Dear R users,
>>>>>>>>>>>>>>  Does anyone has implemented the inverse of the
>>>>>>>>>>>>>> matrix > > > > >  exponential (expm
>>>>>>>>>  in the package Matrix)?
>>>>>>>>>>>>>>  In Matlab, there're logm and expm, there's only expm
>>>>>>>>>>>>>> in R.
>>>>>>>>>  Cheers
>>>>>>>>>  Mimosa
>>>>>>>>>>>>>>>>>>>>>>>>         [[alternative HTML
>>>>>>>>>>>>>>>>>>>>>>>> version deleted]]
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>  ______________________________________________
>>>>>>>>>  R-help at r-project.org mailing list
>>>>>>>>>  https://stat.ethz.ch/mailman/listinfo/r-help
>>>>>>>>>  PLEASE do read the posting guide
>>>>>>>>>  http://www.R-project.org/posting-guide.html
>>>>>>>>>  and provide commented, minimal, self-contained, reproducible >
>>>>>>>>>>>>>  code.
>>>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>>>  ______________________________________________
>>>>>>>>  R-help at r-project.org mailing list
>>>>>>>>  https://stat.ethz.ch/mailman/listinfo/r-help
>>>>>>>>  PLEASE do read the posting guide
>>>>>>>>  http://www.R-project.org/posting-guide.html
>>>>>>>>  and provide commented, minimal, self-contained, reproducible
>>>>>>>> code.
>>>>>>>>>>>>>>>  Charles C. Berry                            (858)
>>>>>>>>>>>>>>> 534-2098
>>>>>>>                                             Dept of
>>>>>>> Family/Preventive
>>>>>>>  Medicine
>>>>>>>  E mailto:cberry at tajo.ucsd.edu               UC San Diego
>>>>>>>  http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego
>>>>>>>  92093-0901
>>>>>>>>>>>>>>>  ______________________________________________
>>>>>>  R-help at r-project.org mailing list
>>>>>>  https://stat.ethz.ch/mailman/listinfo/r-help
>>>>>>  PLEASE do read the posting guide
>>>>>>  http://www.R-project.org/posting-guide.html
>>>>>>  and provide commented, minimal, self-contained, reproducible code.
>>>>>>>>>  Charles C. Berry                            (858) 534-2098
>>>>>                                             Dept of Family/Preventive
>>>>>  Medicine
>>>>>  E mailto:cberry at tajo.ucsd.edu               UC San Diego
>>>>>  http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego >
>>>>>  92093-0901
>>>>>>>
>>>>  ______________________________________________
>>>>  R-help at r-project.org mailing list
>>>>  https://stat.ethz.ch/mailman/listinfo/r-help
>>>>  PLEASE do read the posting guide
>>>>  http://www.R-project.org/posting-guide.html
>>>>  and provide commented, minimal, self-contained, reproducible code.
>>>>
>>>>
>>>
>>>
>>> --
>>> Spencer Graves, PE, PhD
>>> President and Chief Operating Officer
>>> Structure Inspection and Monitoring, Inc.
>>> 751 Emerson Ct.
>>> San José, CA 95126
>>> ph:  408-655-4567
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>> Charles C. Berry                            (858) 534-2098
>>                                            Dept of Family/Preventive
>> Medicine
>> E mailto:cberry at tajo.ucsd.edu               UC San Diego
>> http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901
>>
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901



More information about the R-help mailing list