Cholesky Decomposition

> cholesky <- function(A) {
+     m <- dim(A)[1]
+     if (m > 1) {
+         i <- 1
+         j <- i + 1
+         A.new <- A[j:m, j:m] - A[j:m, i] %*% t(A[j:m, i])/A[i,
+             i]
+         R <- diag(m)
+         R[1, 1] <- sqrt(A[1, 1])
+         R[2:m, 1] <- A[2:m, 1]/sqrt(A[1, 1])
+         res <- R[, 1]
+     }
+     else {
+         res <- sqrt(A[1, 1])
+         A.new <- 0
+     }
+     return(list(R = res, A = A.new))
+ }

> set.seed(12345)
> R1 <- matrix(rnorm(9), 3, 3)
> A <- t(R1) %*% R1
> M <- A
> R <- diag(m)
> for (i in 1:m) {
+     cholt <- choltest(M)
+     X <- cholt$R
+     M <- cholt$A
+     R[i:m, i] <- X
+ }
> round(R %*% t(R) - A, 3)
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0