If my memory is correct, the archives of this list contains
several discussions of round off error problems associated with
different methods for computing things like this. The "Matrix" package
(part of the base distribution) contains a function "expm", whose help
file says, "The expm package contains newer (partly faster and more
accurate) algorithms for expm() and includes logm and sqrtm." This
suggests to me that the most numerically stable way to get an arbitrary
power p of a matrix M in R might be as follows:
expm(p*logm(M))
If compute time becomes an issue, I would want to do numerical
comparisons of the results from this with the results from the your
other method.
Hope this helps.
Spencer
p.s. I found the above in part by using findFn('expm') after
"library(sos)".
On 3/11/2012 9:14 AM, Joshua Wiley wrote:
On Sun, Mar 11, 2012 at 8:56 AM, Peter Langfelder
<peter.langfel...@gmail.com> wrote:
On Sun, Mar 11, 2012 at 1:46 AM, Ebrahim Jahanshiri
<e.jahansh...@gmail.com> wrote:
Dear list,
I understand that to raise matrix A to power (-1/2) we should use something
like this:
eigen(A)$vectors%*%diag(1/sqrt(eigen(A)$values))%*%t(eigen(A)$vectors)
[from previous discussions:
http://r.789695.n4.nabble.com/matrix-power-td900335.html]
But this will only do it for negative sqrt of the matrix not for other
fraction powers like (-3/2).
Not sure why you think this won't work for -3/2 - simply use
eigen(A)$values^(-3/2) instead of the 1/sqrt(eigen(A)$values) and
you're good to go. Generalizations to other powers are left as
exercise for the reader :)
also please be kind to your poor computer and store and reuse
eigen(A), that is not trivial to compute!
Peter
______________________________________________
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 Technology Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph: 408-655-4567
web: www.structuremonitoring.com
______________________________________________
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.
