[R] implementation of matrix logarithm (inverse of matrix exponential)
Charles C. Berry
cberry at tajo.ucsd.edu
Sun Sep 27 18:44:32 CEST 2009
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
More information about the R-help
mailing list