[R] implementation of matrix logarithm (inverse of matrix exponential)
Martin Maechler
maechler at stat.math.ethz.ch
Mon Sep 28 13:01:33 CEST 2009
>>>>> "DB" == Douglas Bates <bates at stat.wisc.edu>
>>>>> on Sun, 27 Sep 2009 17:25:39 -0500 writes:
DB> There is a logm function in the expm package in the expm project on
DB> R-forge. See http://expm.R-forge.R-project.org/
DB> Martin was the person who added that function so I will defer to his
DB> explanations of what it does. I know he has been traveling and it may
DB> be a day or two before he can get to this thread.
Indeed, thank you, Doug.
Yes, I'd strongly recommend using the expm package's logm().
Computing Matrix Exponentials and Logs (and more) reliably is actually a
quite a science.
The recent addition of expm::logm() actually happened based on
a master thesis done here at ETH, and that itself was based
on *THE* reference on these topics,
Higham, N.~J. (2008). _Functions of Matrices: Theory and
Computation_; Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA.
and so is are part of latest methods for the matrix exponential,
expm().
One reason that I haven't released 'expm' to CRAN yet,
has been that it hasn't been entirely clear to me what the
default 'method' should be.
The (currently only one) method, uses in Matrix::expm() is
no longer state of the art (and yes, I had plans to also
"update" that).
An interesting side-question is actually how to deal with such
issues, of active research rendering "the state of the art" to
a moving target.
In Matlab, I see they just provide "the current" method, so are
not bug-back-compatible with previous versions of itself.
In R, I'd at least want to provide an optional 'method' argument
and keep former methods available.
The question still remain if it's okay to change the default
method, as the state of the art advances.
For the matrix functions expm(), logm(), etc., I'd say "yes",
the default should be allowed to change.
Martin Maechler, ETH Zurich
DB> 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.
>>
>>
DB> ______________________________________________
DB> R-help at r-project.org mailing list
DB> https://stat.ethz.ch/mailman/listinfo/r-help
DB> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
DB> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list