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

Douglas Bates bates at stat.wisc.edu
Mon Sep 28 00:25:39 CEST 2009


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.

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




More information about the R-help mailing list