Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
DB == Douglas Bates ba...@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 cbe...@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 cbe...@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 cbe...@tajo.ucsd.edu wrote: On
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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 cbe...@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,]21 [2,]12 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,]21 [2,]12 On Sat, Sep 26, 2009 at 6:24 PM, Charles C. Berry cbe...@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 mimosa1...@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@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@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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. --
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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 cbe...@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 cbe...@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 cbe...@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 mimosa1...@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-h...@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-h...@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ r-h...@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
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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,]10 [2,]11 logm(expm(M)) [,1] [,2] [1,]00 [2,]10 And thanks to Vincent, Martin, et al for this package! :-) Chuck On Sun, Sep 27, 2009 at 11:44 AM, Charles C. Berry cbe...@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 cbe...@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 cbe...@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 mimosa1...@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-h...@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-h...@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ r-h...@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
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
Try: expm( - M) On Sat, Sep 26, 2009 at 5:06 PM, Mimosa Zeus mimosa1...@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@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@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.
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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 mimosa1...@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@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@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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.
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
OK. Try this: library(Matrix) M - matrix(c(2, 1, 1, 2), 2); M [,1] [,2] [1,]21 [2,]12 # log of expm(M) is original matrix M with(eigen(expm(M)), vectors %*% diag(log(values)) %*% t(vectors)) [,1] [,2] [1,]21 [2,]12 On Sat, Sep 26, 2009 at 6:24 PM, Charles C. Berry cbe...@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 mimosa1...@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@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@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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.
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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,]21 [2,]12 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,]21 [2,]12 On Sat, Sep 26, 2009 at 6:24 PM, Charles C. Berry cbe...@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 mimosa1...@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@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@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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.
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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 cbe...@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 cbe...@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 mimosa1...@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@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@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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.
Re: [R] implementation of matrix logarithm (inverse of matrix exponential)
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) 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. 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 cbe...@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,]21 [2,]12 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,]21 [2,]12 On Sat, Sep 26, 2009 at 6:24 PM, Charles C. Berry cbe...@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 mimosa1...@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@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@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@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@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.