1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
|
root.matrix <- function(X)
{
if((ncol(X) == 1L)&&(nrow(X) == 1L)) return(sqrt(X))
else
{
X.eigen <- eigen(X, symmetric=TRUE)
if(any(X.eigen$values < 0)) stop("matrix is not positive semidefinite")
sqomega <- sqrt(diag(X.eigen$values))
V <- X.eigen$vectors
V <- V %*% sqomega %*% t(V)
dimnames(V) <- dimnames(X)
return(V)
}
}
solveCrossprod <- function(X, method = c("qr", "chol", "solve")) {
switch(match.arg(method),
"qr" = chol2inv(qr.R(qr(X))),
"chol" = chol2inv(chol(crossprod(X))),
"solve" = solve(crossprod(X)))
}
|
