最小二乗法による実験データ解析（東京大学出版会） - R入門

目次

コレスキー分解を行う（p.81）

以下は、式(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

---