[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