［R］前進代入法で連立方程式を解く

forwardsolve関数を使う。以下は、連立方程式を行列で表記し、上三角行列の係数行列をA、定数項行列をy、求める解をxとおき、x = Ayと考えてxを求めた例。関数だけではなく、solve関数を使った例、逆行列を求めて正規方程式を解いた例も示してあるが、得られた解は同じであることがわかる。

> set.seed(7)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxaa[upper.tri(mxaa)] <- 0
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    9    0    0    0
[2,]    4    8    0    0
[3,]    2    4    2    0
[4,]    1    9    3    1
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxy)
     [,1]
[1,]    6
[2,]    1
[3,]    9
[4,]    3
> forwardsolve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> solve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> print(solve(t(mxaa) %*% mxaa) %*% (t(mxaa) %*% mxy))
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667

forwardsolve関数は、与えられた行列の下三角の成分だけを使って解く。以下のとおり、すべての成分を使って解くsolve関数とは得られた解が違うことがわかる。

> set.seed(7)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    9    3    2    7
[2,]    4    8    5    1
[3,]    2    4    2    5
[4,]    1    9    3    1
> print(mxy)
     [,1]
[1,]    6
[2,]    1
[3,]    9
[4,]    3
> forwardsolve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> solve(mxaa, mxy)
            [,1]
[1,] -0.97516556
[2,]  0.06125828
[3,]  0.49337748
[4,]  1.94370861

［R］相関行列を求める

cor関数を使う。以下は、あらかじめ用意した行列を使って相関行列を求めた例。参考に手計算でも求めている。組み込み関数であるvar関数は、標本分散ではなく不偏分散（平均値からの偏差の平方和を「標本数-1」で割った値）が求まることに注意。

> # 「情報量規準による統計解析入門」表12.1(p.154)
> # 1行1標本
> print(xx)
      [,1] [,2] [,3] [,4] [,5]
 [1,]  165   53   86   56   92
 [2,]  160   47   84   52   92
 [3,]  166   55   86   64   89
 [4,]  164   56   90   60   95
 [5,]  168   55   87   56   87
 [6,]  164   54   87   57   92
 [7,]  168   54   94   58   97
 [8,]  169   55   88   57   92
 [9,]  169   53   86   58   93
[10,]  166   56   84   57   90
> # 組み込み関数corにより求めた相関行列
> print(cor(xx))
            [,1]       [,2]      [,3]         [,4]         [,5]
[1,]  1.00000000  0.6173717 0.3620673  0.368208142 -0.043119609
[2,]  0.61737166  1.0000000 0.3341712  0.665307665 -0.107443062
[3,]  0.36206726  0.3341712 1.0000000  0.280441930  0.685251832
[4,]  0.36820814  0.6653077 0.2804419  1.000000000 -0.006370564
[5,] -0.04311961 -0.1074431 0.6852518 -0.006370564  1.000000000
> 
> # 以下、手計算により相関行列を求める、nに標本数(行列の行数)を代入
> n <- nrow(xx)
> # 各列の平均
> dmean <- apply(xx, 2, mean)
> # 各列の標本分散(平均からの偏差の平方和をnで割る、不偏分散ではない)
> dvar <- apply(xx, 2, function(x) (length(x) - 1) / length(x) * var(x))
> # 各列の標本分散の平方根
> dstde <- sqrt(dvar)
> # 列ごとに平均を引いて標本分散で割って標準化する
> xx0 <- sweep(xx, 2, dmean, FUN = "-")
> xx0 <- sweep(xx0, 2, dstde, FUN = "/")
> # 相関行列を求める
> rr <- t(xx0) %*% xx0 / n
> print(rr)
            [,1]       [,2]      [,3]         [,4]         [,5]
[1,]  1.00000000  0.6173717 0.3620673  0.368208142 -0.043119609
[2,]  0.61737166  1.0000000 0.3341712  0.665307665 -0.107443062
[3,]  0.36206726  0.3341712 1.0000000  0.280441930  0.685251832
[4,]  0.36820814  0.6653077 0.2804419  1.000000000 -0.006370564
[5,] -0.04311961 -0.1074431 0.6852518 -0.006370564  1.000000000

［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

［R］下三角行列を作成する

upper.tri関数に行列を指定すると、対角成分より上の成分（各成分を(i,j)と表記すればi<jの成分）だけTRUEを返す。これを利用してその成分に0を代入すれば、既存の行列を簡単に下三角行列（対角成分より上の成分はすべて0の行列）に変換することができる。

> set.seed(5)
> mx <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    2    1    9    3
[2,]    7    7    1    6
[3,]    9    5    3    3
[4,]    3    8    5    2
> print(upper.tri(mx))
      [,1]  [,2]  [,3]  [,4]
[1,] FALSE  TRUE  TRUE  TRUE
[2,] FALSE FALSE  TRUE  TRUE
[3,] FALSE FALSE FALSE  TRUE
[4,] FALSE FALSE FALSE FALSE
> mx[upper.tri(mx)] <- 0
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    2    0    0    0
[2,]    7    7    0    0
[3,]    9    5    3    0
[4,]    3    8    5    2

［R］上三角行列を作成する

lower.tri関数に行列を指定すると、対角成分より下の成分（各成分を(i,j)と表記すればi>jの成分）だけTRUEを返す。これを利用してその成分に0を代入すれば、既存の行列を簡単に上三角行列（対角成分より下の成分はすべて0の行列）に変換することができる。

> set.seed(4)
> mx <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    6    8    9    1
[2,]    1    3    1    9
[3,]    3    7    7    4
[4,]    3    9    3    5
> print(lower.tri(mx))
      [,1]  [,2]  [,3]  [,4]
[1,] FALSE FALSE FALSE FALSE
[2,]  TRUE FALSE FALSE FALSE
[3,]  TRUE  TRUE FALSE FALSE
[4,]  TRUE  TRUE  TRUE FALSE
> mx[lower.tri(mx)] <- 0
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    6    8    9    1
[2,]    0    3    1    9
[3,]    0    0    7    4
[4,]    0    0    0    5

［R］QR分解を行う

qr関数を使う。分解した行列はqr関数の戻り値からそれぞれqr.Q関数、qr.R関数を使うことで得ることができる。

m行n列の行列A(m >= n)が与えられたとき、Q^T Q = Iを満たす行列Qと上三角行列（対角成分より下の成分はすべて0の行列）Rを用いてA = Q Rと分解することを、行列AのQR分解という。QR分解はcomplete型とreduced型の2種類があり、それぞれ分解した行列の行数と列数が異なる。

complete型 A(m×n) = Q(m×m) R(m×n)
reduced型 A(m×n) = Q(m×n) R(n×n)

以下は適当な行列を作成してQR分解した例。Q^T Q = Iを満たしていることと、Q R = Aと元に戻ることも計算している。qr.Q関数とqr.R関数にそれぞれcompleteオプションにTRUEを指定すると、complete型の計算を行う。何も指定しなければreduced型の計算結果が返される。

> m <- 4
> n <- 3
> aa <- matrix(0.0, m, n)
> for (i in 1:m) for (j in 1:n) aa[i, j] <- (-1) ^ i * choose(2 * i + j, i)
> print(aa)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330
> 
> # complete型 A(m×n) = Q(m×m) R(m×n)
> r <- qr(aa)
> qq <- qr.Q(r, complete = TRUE)  # m×m
> rr <- qr.R(r, complete = TRUE)  # m×n
> print(qq)
            [,1]       [,2]       [,3]      [,4]
[1,] -0.02286814  0.3344769  0.7843075 0.5219808
[2,]  0.07622713 -0.5474374 -0.2854080 0.7829713
[3,] -0.26679495  0.7243109 -0.5399946 0.3355591
[4,]  0.96046183  0.2526086 -0.1086731 0.0434984
> print(rr)
         [,1]       [,2]        [,3]
[1,] 131.1869 217.872381 341.0782902
[2,]   0.0000   2.936931   9.3501598
[3,]   0.0000   0.000000  -0.4176769
[4,]   0.0000   0.000000   0.0000000
> print(t(qq) %*% qq)
              [,1]          [,2]          [,3]          [,4]
[1,]  1.000000e+00 -8.326673e-17 -2.220446e-16 -7.632783e-17
[2,] -8.326673e-17  1.000000e+00  9.714451e-17  2.775558e-17
[3,] -2.220446e-16  9.714451e-17  1.000000e+00  1.257675e-16
[4,] -7.632783e-17  2.775558e-17  1.257675e-16  1.000000e+00
> print(qq %*% rr)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330
> 
> # reduced型 A(m×n) = Q(m×n) R(n×n)
> r <- qr(aa)
> qq <- qr.Q(r)  # m×n
> rr <- qr.R(r)  # n×n
> print(qq)
            [,1]       [,2]       [,3]
[1,] -0.02286814  0.3344769  0.7843075
[2,]  0.07622713 -0.5474374 -0.2854080
[3,] -0.26679495  0.7243109 -0.5399946
[4,]  0.96046183  0.2526086 -0.1086731
> print(rr)
         [,1]       [,2]        [,3]
[1,] 131.1869 217.872381 341.0782902
[2,]   0.0000   2.936931   9.3501598
[3,]   0.0000   0.000000  -0.4176769
> print(t(qq) %*% qq)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -8.326673e-17 -2.220446e-16
[2,] -8.326673e-17  1.000000e+00  9.714451e-17
[3,] -2.220446e-16  9.714451e-17  1.000000e+00
> print(qq %*% rr)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330

qr関数にLAPACKオプションでTRUEを指定すると（デフォルトはFALSE)、計算にLAPACKを使いピボット選択が行われる（Using LAPACK (including in the complex case) uses column pivoting and does not attempt to detect rank-deficient matrices.）。このピボットの情報は、qr関数の戻り値のpivotというベクトルに含まれている。

> print(aa)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330
> r <- qr(aa, LAPACK = TRUE)
> print(r$pivot)
[1] 3 1 2
> qq <- qr.Q(r)
> rr <- qr.R(r)
> print(qq)
            [,1]       [,2]       [,3]
[1,] -0.01465387  0.2996579  0.7984525
[2,]  0.06154627 -0.5360454 -0.3095536
[3,] -0.24618510  0.7547242 -0.5071335
[4,]  0.96715574  0.2307639 -0.0972919
> print(rr)
         [,1]       [,2]        [,3]
[1,] 341.2067 131.137526 217.8708796
[2,]   0.0000  -3.598527  -3.0434558
[3,]   0.0000   0.000000   0.1310637
> print(t(qq) %*% qq)
             [,1]          [,2]          [,3]
[1,] 1.000000e+00  5.551115e-17  1.387779e-16
[2,] 5.551115e-17  1.000000e+00 -9.020562e-17
[3,] 1.387779e-16 -9.020562e-17  1.000000e+00
> print(aa[, r$pivot])
     [,1] [,2] [,3]
[1,]   -5   -3   -4
[2,]   21   10   15
[3,]  -84  -35  -56
[4,]  330  126  210
> print(qq %*% rr)
     [,1] [,2] [,3]
[1,]   -5   -3   -4
[2,]   21   10   15
[3,]  -84  -35  -56
[4,]  330  126  210

［R］LU分解を行う

Matrixパッケージのlu関数を使う。LU分解したそれぞれの行列を取り出すには、同パッケージに含まれているexpand関数を使う。

n次の正方行列Aが与えられたとき、下三角行列Lと上三角行列Uを用いてA = L Uと分解することを、行列AのLU分解という。lu関数は置換行列Pを用いてP A = L Uと分解する。戻り値はLが下三角行列、Uが上三角行列、Pが置換行列。以下は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(Matrix)
> r <- expand(lu(aa))
> print(r$L)
4 x 4 Matrix of class "dtrMatrix" (unitriangular)
     [,1]      [,2]      [,3]      [,4]     
[1,] 1.0000000         .         .         .
[2,] 1.0000000 1.0000000         .         .
[3,] 1.0000000 0.6666667 1.0000000         .
[4,] 1.0000000 0.3333333 1.0000000 1.0000000
> print(r$U)
4 x 4 Matrix of class "dtrMatrix"
     [,1]       [,2]       [,3]       [,4]      
[1,]  1.0000000  1.0000000  1.0000000  1.0000000
[2,]          .  3.0000000  9.0000000 19.0000000
[3,]          .          . -1.0000000 -3.6666667
[4,]          .          .          .  0.3333333
> print(r$P)
4 x 4 sparse Matrix of class "pMatrix"
[1,] | . . .
[2,] . . . |
[3,] . . | .
[4,] . | . .
> print(solve(r$P) %*% r$L %*% r$U)
4 x 4 Matrix of class "dgeMatrix"
     [,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］パスカル行列を求める

パスカル行列は正定値対称行列であり、下三角コレスキー因子（パスカル行列P = L L'を満たす下三角行列L）を求めて計算することで作成している。以下は、lapply関数を使用して5次の正方行列であるパスカル行列を作成した例。

> n <- 5
> lis <- lapply(0:(n - 1), function(x) choose(x, 0:x))
> ll <- diag(0, n)
> for (i in 1:n) for (j in 1:i) ll[i, j] <- lis[[i]][j]
> print(ll)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    1    1    0    0    0
[3,]    1    2    1    0    0
[4,]    1    3    3    1    0
[5,]    1    4    6    4    1
> pp <- ll %*% t(ll)
> print(pp)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    2    3    4    5
[3,]    1    3    6   10   15
[4,]    1    4   10   20   35
[5,]    1    5   15   35   70

［R］行列の要素を抜き出してベクトルにする

as.vector関数を使う。行列の列ベクトルごとに順に要素を抜き出すため、行ベクトルごとに抜き出したい場合は、引数に与える行列を転置行列にしておくこと。

> mx <- matrix(1:12, nrow = 3, ncol = 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> v <- as.vector(mx)
> print(v)
 [1]  1  2  3  4  5  6  7  8  9 10 11 12
> v <- as.vector(t(mx))
> print(v)
 [1]  1  4  7 10  2  5  8 11  3  6  9 12

［R］行列のアダマール積を求める

*演算子を使う。それぞれ行列を指定する。アダマール積の定義（行の数と列の数が同じ行列同士の各成分を掛ける）から、二つの行列の行数と列数が一致しない場合はエラーが発生する。

> mxxx <- matrix(c(1, 2, 2, 2, 3, 5, 3, 4, 6), nrow = 3, ncol = 3)
> mxxx
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    4
[3,]    2    5    6
> mxxx * mxxx
     [,1] [,2] [,3]
[1,]    1    4    9
[2,]    4    9   16
[3,]    4   25   36
> mxaa <- matrix(c(1, 2), nrow = 1, ncol = 2)
> mxbb <- matrix(c(1, 2), nrow = 2, ncol = 1)
> mxaa
     [,1] [,2]
[1,]    1    2
> mxbb
     [,1]
[1,]    1
[2,]    2
> mxaa * mxbb
mxaa * mxbb でエラー: 適切な配列ではありません