There are a couple of ways to do matrix square roots in R. The first is through eigendecomposition:
# a is the matrix for which we want to square root:
a.eig <- eigen(a)
a.sqrt <- a.eig$vectors %*% diag(sqrt(a.eig$values)) %*% solve(a.eig$vectors)
Of course, this only works for positive semidefinite matrices that are diagonalizable (no normal Jordan forms, etc.)
Another way is to use the Denman-Beavers algorithm:
denman.beavers <- function(mat,maxit=50) {
stopifnot(nrow(mat) == ncol(mat))
niter <- 0
y <- mat
z <- diag(rep(1,nrow(mat)))
for (niter in 1:maxit) {
y.temp <- 0.5*(y+solve(z))
z <- 0.5*(z+solve(y))
y <- y.temp
}
return(list(sqrt=y,sqrt.inv=z))
}
This can probably be improved with better convergence criteria, but so far it works well and gives both the square root and inverse square root, both of which are required for many applications.
