Rで線形代数:コレスキー分解

Rでデータサイエンス

コレスキー分解

\(A\)が実正定値対称行列\((A^t=A,\,x^tAx>0\quad \forall x(\neq0)\in\mathbb{R}^n)\)ならば対角要素がすべて正の下三角行列\(L\)をえらんで\[A=LL^t\]と書ける。この分解は\(L\)が単位下三角行列でなくとも一意である。

  • 参考引用資料: 山本哲朗(2010),『行列解析の基礎』,サイエンス社,pp.63-64
library(dplyr)
func_matrix_decomposition_cholesky <- function(A) {
    n <- nrow(A)
    l <- matrix(nrow = n, ncol = n)
    for (ccc in 1:n) {
        # i=j part
        ll <- A[ccc, ccc]
        buf <- 0
        for (kkk in 1:ccc) {
            if (kkk < ccc) {
                buf <- l[ccc, kkk] * l[ccc, kkk] + buf
            }
        }
        l[ccc, ccc] <- sqrt(ll - buf)
        # i>j part 但し (ccc+1):n とせずに ccc:n としているため rrr=ccc の分、一回重複する。
        for (rrr in ccc:n) {
            ll <- A[rrr, ccc]
            buf <- 0
            for (kkk in 1:ccc) {
                if (kkk < ccc) {
                  buf <- l[rrr, kkk] * l[ccc, kkk] + buf
                }
            }
            l[rrr, ccc] <- (ll - buf)/l[ccc, ccc]
        }
    }
    l[is.na(l)] <- 0
    llt <- l %*% t(l)
    return(list(l = l, llt = llt, A = A))
}
# 例1
A <- c(8, 1, 5, 1, 3, 2, 5, 2, 5) %>%
    matrix(nrow = 3, byrow = T)
# eigen(A)$values > 0
func_matrix_decomposition_cholesky(A = A)
$l
          [,1]      [,2]     [,3]
[1,] 2.8284271 0.0000000 0.000000
[2,] 0.3535534 1.6955825 0.000000
[3,] 1.7677670 0.8109308 1.103355

$llt
     [,1] [,2] [,3]
[1,]    8    1    5
[2,]    1    3    2
[3,]    5    2    5

$A
     [,1] [,2] [,3]
[1,]    8    1    5
[2,]    1    3    2
[3,]    5    2    5
# 例2
A <- c(2, -1, 0, -1, 2, -1, 0, -1, 2) %>%
    matrix(nrow = 3, byrow = T)
# eigen(A)$values > 0
func_matrix_decomposition_cholesky(A = A)
$l
           [,1]       [,2]     [,3]
[1,]  1.4142136  0.0000000 0.000000
[2,] -0.7071068  1.2247449 0.000000
[3,]  0.0000000 -0.8164966 1.154701

$llt
     [,1] [,2] [,3]
[1,]    2   -1    0
[2,]   -1    2   -1
[3,]    0   -1    2

$A
     [,1] [,2] [,3]
[1,]    2   -1    0
[2,]   -1    2   -1
[3,]    0   -1    2
# 例3
A <- c(5, 7, 6, 5, 7, 10, 8, 7, 6, 8, 10, 9, 5, 7, 9, 10) %>%
    matrix(nrow = 4, byrow = T)
# eigen(A)$values > 0
func_matrix_decomposition_cholesky(A = A)
$l
         [,1]       [,2]     [,3]      [,4]
[1,] 2.236068  0.0000000 0.000000 0.0000000
[2,] 3.130495  0.4472136 0.000000 0.0000000
[3,] 2.683282 -0.8944272 1.414214 0.0000000
[4,] 2.236068  0.0000000 2.121320 0.7071068

$llt
     [,1] [,2] [,3] [,4]
[1,]    5    7    6    5
[2,]    7   10    8    7
[3,]    6    8   10    9
[4,]    5    7    9   10

$A
     [,1] [,2] [,3] [,4]
[1,]    5    7    6    5
[2,]    7   10    8    7
[3,]    6    8   10    9
[4,]    5    7    9   10
# 対称行列の作り方
# 正定値とは限らない
n <- 4
A <- matrix(data = 0, nrow = n, ncol = n)
num <- sample(seq(10), choose(n, 2))
A[lower.tri(A)] <- num
A <- t(A) + A
diag(A) <- sample(seq(10), n)
A
     [,1] [,2] [,3] [,4]
[1,]    4   10    3    7
[2,]   10    2    4    2
[3,]    3    4    9    9
[4,]    7    2    9    7

参考引用資料

  1. 山本哲朗(2010),『行列解析の基礎』,サイエンス社,pp.63-65
  2. http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/lecture/mic/text/02_linearsystem1.pdf
  3. https://www.math.utah.edu/~zwick/Classes/Fall2012_2270/Lectures/Lecture33_with_Examples.pdf
  4. https://nhigham.com/2020/07/21/what-is-a-symmetric-positive-definite-matrix/
  5. https://davetang.org/muse/2014/01/22/making-symmetric-matrices-in-r/

最終更新

Sys.time()
[1] "2024-03-27 14:54:17 JST"

R、Quarto、Package

R.Version()$version.string
[1] "R version 4.3.3 (2024-02-29 ucrt)"
quarto::quarto_version()
[1] '1.4.542'
packageVersion(pkg = "tidyverse")
[1] '2.0.0'

著者