[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