[R] Matrix negative fraction power

Spencer Graves 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:


             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.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:
>>>
>>> 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 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 Technology Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph:  408-655-4567
web:  www.structuremonitoring.com



More information about the R-help mailing list