Skip to content
Prev 198747 / 398506 Next

matrix^(-1/2)

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: