［R］LU分解を行う

matrixcalcパッケージのlu.decomposition関数を使う。n次の正方行列Aが与えられたとき、下三角行列Lと上三角行列Uを用いてA = L Uと分解することを、行列AのLU分解という。戻り値はLが下三角行列、Uが上三角行列。lu.decomposition関数はピボット選択は行わないため、戻り値に置換行列はない。以下は4次の正方行列のパスカル行列をLU分解した例。最後に、元の行列に戻るかどうか計算している。

> n <- 4
> aa <- matrix(1, n, n)
> for (i in 2:n) for (j in 2:n) aa[i, j] <- aa[i, j - 1] + aa[i - 1, j]
> print(aa)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20
> library(matrixcalc)
> r <- lu.decomposition(aa)
> print(r$L)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    1    0
[4,]    1    3    3    1
> print(r$U)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    0    1    2    3
[3,]    0    0    1    3
[4,]    0    0    0    1
> print(r$L %*% r$U)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20