On 11-03-2012, at 18:18, Berend Hasselman wrote: > > On 11-03-2012, at 17:52, Spencer Graves wrote: > >> 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)) >> > > I have just tested and it should be this > > library(expm) > > expm((-1/2)*logm(A)) > expm((-3/2)*logm(A)) >
Indeed. It's getting on and I made a muddle. Berend ______________________________________________ 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.