[R] matrix^(-1/2)

[spencerg]

  A question, a comment, and an alternative answer to matrix^(-1/2):


QUESTION:


What's the status of the "expm" package, mentioned in the email you 
cited from Martin Maechler, dated Apr 5 19:52:09 CEST 2008? I tried both 
install.packages('expm') and 
install.packages("expm",repos="http://R-Forge.R-project.org"), and got 
"package 'expm' is not available" in both cases.


COMMENT:


The solution proposed by Venables rests on Sylvester's matrix theorem, 
which essentially says that if a matrix A is diagonalizable with 
eigenvalue decomposition eigA <- eigen(A) and f: D → C with D ⊂ C be a 
function for which f(A) is well defined 
(http://en.wikipedia.org/wiki/Sylvester%27s_matrix_theorem), then f(A) = 
with(eigA, vectors %*% diag(f(values)) %*% solve(vectors)). Maechler and 
others have noted that this can be one of the least accurate and most 
computationally expensive ways to compute f(A).


ALTERNATIVE ANSWER:


For A^(-1/2), if A is symmetric and nonnegative definite, then 
solve(chol(A)) would be a very good way to compute it.


Hope this helps,
Spencer


David Winsemius wrote:
>
> On Oct 31, 2009, at 9:33 PM, David Winsemius wrote:
>
>>
>> On Oct 31, 2009, at 4:39 PM, Kajan Saied wrote:
>>
>>> Dear R-Help Team,
>>>
>>> as a R novice I have a (maybe for you very simple question), how do 
>>> I get
>>> the following solved in R:
>>>
>>> Let R be a n x n matrix:
>>>
>>> \mid R\mid^{-\frac{1}{2}}
>>>
>>> solve(A) gives me the inverse of the matrix R, however not the 
>>> ^(-1/2) of
>>> the matrix...
>>
>> GIYF: (and Bill Venables if friendly, too.)
>>
>> http://www.lmgtfy.com/?q=powers+of+matrix+r-project
>
> I had assumed that the first hit I got:
>
> https://stat.ethz.ch/pipermail/r-help/2008-April/160662.html
>
> ... would be the first hit anybody got, but that's not necessarily 
> true now and especially for the future. And further searching within 
> the results produced this more recent Maechler posting:
>
> https://stat.ethz.ch/pipermail/r-devel/2008-April/048969.html
>
> For the Mac users, there appears to be no binary, but the source 
> compiles without error on a 64-bit version of R 2.10.0:
>
> install.packages("expm",repos="http://R-Forge.R-project.org", 
> type="source")
>
> #The suggested code throws an error, so my very minor revision would be:
>
> library(expm)
> ?"%^%"
>


-- 
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