以下は、式(4.115)に基に作成したコレスキー分解を行う関数。行列 B を上三角行列 L の転置行列と上三角行列の積に分解している。細かなチェックは行っていないので注意。
B = R^T R
choldecomp <- function(mxbb) {
m <- ncol(mxbb)
mxrr <- matrix(0., m, m)
mxrr[1, 1] <- sqrt(mxbb[1, 1])
for (j in 2:m) {
mxrr[1, j] <- mxbb[1, j] / mxrr[1, 1]
}
for (i in 2:m) {
k <- 1:(i - 1)
mxrr[i, i] <- sqrt(mxbb[i, i] - sum(mxrr[1:(i - 1), i] ^ 2))
for (j in i:m) {
mxrr[i, j] <- (mxbb[i, j]
- sum(mxrr[k, i] * mxrr[k, j])) / mxrr[i, i]
}
}
return(mxrr)
}
以下、動作確認。
> d <- c(1, 1, 1, 1, 1, 2, 3, 4, 1, 3, 6, 10, 1, 4, 10, 20)
> mxbb <- matrix(d, 4, 4, byrow = TRUE)
> mxrr <- choldecomp(mxbb)
> cat("B\n")
B
> print(mxbb)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 3 4
[3,] 1 3 6 10
[4,] 1 4 10 20
> cat("R^T\n")
R^T
> print(t(mxrr))
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 1 1 0 0
[3,] 1 2 1 0
[4,] 1 3 3 1
> cat("R\n")
R
> print(mxrr)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 0 1 2 3
[3,] 0 0 1 3
[4,] 0 0 0 1
> cat("R^T T\n")
R^T T
> print(t(mxrr) %*% mxrr)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 3 4
[3,] 1 3 6 10
[4,] 1 4 10 20
Rへの組み込み関数cholでも計算してみる。上の結果（mxrr, R）と一致していることがわかる。
> base::chol(mxbb)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 0 1 2 3
[3,] 0 0 1 3
[4,] 0 0 0 1
