[R] Matrix negative fraction power
spencer.graves at structuremonitoring.com
Sun Mar 11 17:52:40 CET 2012
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:
If compute time becomes an issue, I would want to do numerical
comparisons of the results from this with the results from the your
Hope this helps.
p.s. I found the above in part by using findFn('expm') after
On 3/11/2012 9:14 AM, Joshua Wiley wrote:
> On Sun, Mar 11, 2012 at 8:56 AM, Peter Langfelder
> <peter.langfelder at gmail.com> wrote:
>> On Sun, Mar 11, 2012 at 1:46 AM, Ebrahim Jahanshiri
>> <e.jahanshiri at gmail.com> wrote:
>>> Dear list,
>>> I understand that to raise matrix A to power (-1/2) we should use something
>>> like this:
>>> [from previous discussions:
>>> 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!
>> R-help at r-project.org mailing list
>> 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
More information about the R-help