R入門
行列
最終更新：
r-intro
目次
基本操作
行列を作成する
matrix関数を使う。行列の成分を格納したベクトルを第一引数に与えて、nrowオプションで行数を、ncolオプションに列数を与える。デフォルトでは列番号が小さい列から順に要素が格納される。これを行番号が小さい行から順に格納するにはbyrowオプションにTRUEを指定する。
> d <- 1:12
> mx <- matrix(d, 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
> mx <- matrix(d, ncol = 4)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> mx <- matrix(d, ncol = 5)
警告メッセージ:
matrix(d, ncol = 5) で:
データ長 [12] が列数 [5] を整数で割った、もしくは掛けた値ではありません
> mx <- matrix(d, ncol = 4, byrow = TRUE)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
> mx <- matrix(1:4, nrow = 2, ncol = 4, byrow = TRUE)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 1 2 3 4
> mx <- matrix(0., 2, 5)
> print(mx)
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0例のとおりに、何も指定しなければ第二引数が行列、第三引数が列数になる。ベクトルの要素の数で行列どちらかが決まればもう片方も決まるため、行または列の省略は可能。ただし、ベクトルの要素の数に応じた正しい値でなければならない。指定した行数×列数の要素数に対してベクトルの要素数が足らないときは、そのベクトルの要素数が行数×列数の約数であれば、自動的に繰り返し使われる。
行列を作る
行列はmatrix関数を使用して、ベクトルから作ることができる。以下は、2行3列の行列を新規に作成している。
> d <- 1:6
> mx <- matrix(d, 2, 3)
> mx
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6matrix関数の第一引数には、成分（並べられている各々の数）に代入する値が含まれているベクトルを指定し、特に指定をしなければ第二引数には行数を、第三引数には列数を指定する。行数はnrowオプション、列数はncolオプションで明示的に指定することができ（順不同）、例えば、以下のようにしても同じ。
> mx <- matrix(d, ncol = 3, nrow = 2)
> mx
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6特にオプションを指定しない場合、第一引数に与えたベクトルの要素について、列方向（縦方向）を優先に値を埋めるように代入する。これは見た目の直感となかなか合わない。列方向ではなく、行方向（横方向）を優先に値を埋めるように代入させるには、byrowオプションにTRUEを与える。
> mx <- matrix(d, 2, 3, byrow = TRUE)
> mx
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6行列を作る
matrix関数を使う。以下のように引数を指定する
matrix(成分を含むベクトル, 行数, 列数)行数と列数はnrowオプションとncolオプションで明示的に指定することができる。また、指定したベクトルの要素数に応じて、行または列数だけの指定でもう一方が決まる場合は、指定を省略できる。
> matrix(1:6, 2, 3)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> matrix(1:6, nrow = 2, ncol = 3)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> matrix(1:6, ncol = 3, nrow = 2)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> matrix(1:6, nrow = 3, ncol = 2)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
> matrix(1:6, ncol = 2)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6指定したベクトルの要素数と行数や列数が揃わない場合は、エラーが発生する。
> matrix(1:8, nrow = 2, ncol = 2)
[,1] [,2]
[1,] 1 3
[2,] 2 4
警告メッセージ:
matrix(1:8, nrow = 2, ncol = 2) で:
data length differs from size of matrix: [8 != 2 x 2]
> matrix(1:8, nrow = 3)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 1
警告メッセージ:
matrix(1:8, nrow = 3) で:
data length [8] is not a sub-multiple or multiple of the number of rows [3]matrix関数は指定したベクトルの要素を、行列の一列目の成分から順番に埋めようとする。これを一行目の成分から埋めるようにするには、byrowオプションにTRUEを指定する。
> matrix(1:6, nrow = 2, ncol = 3)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> matrix(1:6, nrow = 2, ncol = 3, byrow = TRUE)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6行列から成分を取り出す
行ごと列ごとに角括弧（ブラケット）をつかってインデックス（1～）を指定する。インデックスは1から始まることに注意。
> mx <- matrix(1:6, 3, 2)
> print(mx)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
> mx[1, 2]
[1] 4
> mx[3, 2]
[1] 6
> mx[4, 2]
mx[4, 2] でエラー: 添え字が許される範囲外です最後の例のように、行列の各長さより大きい値を指定すると、エラーが発生する。
行列から特定の行か特定の列のみ取り出す
インデックスを指定する際に、数値を与えずに空欄のままにする。以下の例ではそれぞれ2行目のみ、3列目のみを取り出している。
> mx <- matrix(1:9, 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> mx[2, ]
[1] 2 5 8
> mx[, 3]
[1] 7 8 9
> is.vector(mx[, 3])
[1] TRUE最後のとおりに、1行（列）だけ取り出したものはベクトルになる。これを行列のままにする場合はdropオプションにFALSEを与える。
> mx[2, , drop = FALSE]
[,1] [,2] [,3]
[1,] 2 5 8複数行（列）取り出したい場合は、インデックスをベクトルで与える。以下は2行目と3行目のみを取り出した例。
> mx[c(2, 3), ]
[,1] [,2] [,3]
[1,] 2 5 8
[2,] 3 6 9行列をベクトルに変換する
as.vector関数を使用する。
> mx1 <- matrix(6, nrow = 1, ncol = 1)
> mx2 <- matrix(1:9, nrow = 3, ncol = 3)
> mx1
[,1]
[1,] 6
> as.vector(mx1)
[1] 6
> mx2
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> as.vector(mx2)
[1] 1 2 3 4 5 6 7 8 9
> as.vector(t(mx2))
[1] 1 4 7 2 5 8 3 6 91列目の要素すべて、2列目の要素すべて、・・・、という順にベクトルの要素になる。これを1行目の要素すべて、2行目の要素すべて、・・・、としたい場合は、最後の例のとおりにt関数を使って元の行列の転置行列を指定する。
行列から非対角成分だけを取り出す
row関数とcol関数を使うと、ベクトル形式で簡単に取り出すことができる。
> mx <- matrix(1:16, 4, 4)
> mx
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 2 6 10 14
[3,] 3 7 11 15
[4,] 4 8 12 16
> mx[row(mx) != col(mx)]
[1] 2 3 4 5 7 8 9 10 12 13 14 15応用として、上の例で2、7、12（対角成分のそれぞれ1行下の成分）を取り除くには、以下のようにする。
> mx[row(mx) != col(mx) + 1]
[1] 1 3 4 5 6 8 9 10 11 13 14 15 16行列に行を追加する
rbind関数を使う。最後の例のとおり、列数が異なる場合はエラーが発生する。
> mx1 <- matrix(1:9, ncol = 3)
> mx2 <- matrix(10:15, ncol = 3)
> print(mx1)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> print(mx2)
[,1] [,2] [,3]
[1,] 10 12 14
[2,] 11 13 15
> mx3 <- rbind(mx1, mx2)
> print(mx3)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
[4,] 10 12 14
[5,] 11 13 15
> mx4 <- matrix(16:19, ncol = 4)
> print(mx4)
[,1] [,2] [,3] [,4]
[1,] 16 17 18 19
> rbind(mx1, mx4)
rbind(mx1, mx4) でエラー:
行列の列数は一致していなければなりません (2 番目の引数を参照)行列に列を追加する
cbind関数を使う。最後の例のとおり、行数が異なる場合はエラーが発生する。
> mx1 <- matrix(1:9, nrow = 3)
> mx2 <- matrix(10:15, nrow = 3)
> print(mx1)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> print(mx2)
[,1] [,2]
[1,] 10 13
[2,] 11 14
[3,] 12 15
> mx3 <- cbind(mx1, mx2)
> print(mx3)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 4 7 10 13
[2,] 2 5 8 11 14
[3,] 3 6 9 12 15
> mx4 <- matrix(16:19, nrow = 4)
> print(mx4)
[,1]
[1,] 16
[2,] 17
[3,] 18
[4,] 19
> cbind(mx1, mx4)
cbind(mx1, mx4) でエラー:
行列の行数は一致していなければなりません (2 番目の引数を参照)演算
行列の和と差を求める
+演算子または-演算子を使う。行列の和と差の定義（行の数と列の数が同じ行列同士の各成分を足す、引く）から、二つの行列の行数と列数が一致しない場合はエラーが発生する。
> mxaa1 <- matrix(c(1, 6, 2, 4, 3, 2), nrow = 2, ncol = 3)
> mxbb1 <- matrix(c(1, 2, 1, 3, 0, 8), nrow = 2, ncol = 3)
> mxaa1
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 6 4 2
> mxbb1
[,1] [,2] [,3]
[1,] 1 1 0
[2,] 2 3 8
> mxaa1 + mxbb1
[,1] [,2] [,3]
[1,] 2 3 3
[2,] 8 7 10
> mxaa2 <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)
> mxbb2 <- matrix(c(7, 5, 9, 1), nrow = 2, ncol = 2)
> mxaa2
[,1] [,2]
[1,] 1 3
[2,] 2 4
> mxbb2
[,1] [,2]
[1,] 7 9
[2,] 5 1
> mxaa2 - mxbb2
[,1] [,2]
[1,] -6 -6
[2,] -3 3
> mxaa1 + mxbb2
mxaa1 + mxbb2 でエラー: 適切な配列ではありません
> mxaa1 - mxbb2
mxaa1 - mxbb2 でエラー: 適切な配列ではありません行列の積を求める
%*%演算子を使う。行列の乗法の定義より、被乗数（左側）の行列の列数と、乗数（右側）の行列の行数が同じでなければ求めることはできずにエラーが発生する。
> aa <- matrix(1:6, 3, 2)
> bb <- matrix(1:8, 2, 4)
> print(aa)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
> print(bb)
[,1] [,2] [,3] [,4]
[1,] 1 3 5 7
[2,] 2 4 6 8
> print(aa %*% bb)
[,1] [,2] [,3] [,4]
[1,] 9 19 29 39
[2,] 12 26 40 54
[3,] 15 33 51 69
> print(bb %*% aa)
bb %*% aa でエラー: 適切な引数ではありません行列の積を求める
%*% 演算子を使う。最後の例のとおり、掛けられる行列の列数と掛ける行列の行数が同じでなければ計算することはできない。
> mx1 <- matrix(c(-1, 8, 0, 1, 3, -5), 2, 3)
> mx2 <- matrix(c(3, -2, 4, 1, 0, 2), 2, 3)
> mx3 <- matrix(c(4, 2, 1, -3, 0, 5, 7, -1, 0), 3, 3)
> print(mx1)
[,1] [,2] [,3]
[1,] -1 0 3
[2,] 8 1 -5
> print(mx2)
[,1] [,2] [,3]
[1,] 3 4 0
[2,] -2 1 2
> print(mx3)
[,1] [,2] [,3]
[1,] 4 -3 7
[2,] 2 0 -1
[3,] 1 5 0
> mx1 %*% mx3
[,1] [,2] [,3]
[1,] -1 18 -7
[2,] 29 -49 55
> mx2 %*% mx3
[,1] [,2] [,3]
[1,] 20 -9 17
[2,] -4 16 -15
> mx1 %*% mx2
mx1 %*% mx2 でエラー: 適切な引数ではありません行列の積を求める
行列Aと行列Bの積ABを求めるには、%*%演算子を使う。行列の乗法は、左側の行列Aの列の数と右側の行列Bの行の数が一致している場合のみ求めることができる（そのような定義になっている）。最後の例のとおり、左右の列数と行数が一致しない場合はエラーが発生する。
> mxaa1 <- matrix(c(1, 2, 2, 2, 3, 5, 3, 4, 6), nrow = 3, ncol = 3)
> mxbb1 <- matrix(c(1, 0, 0, 0, 1, 1, 1, 0, 1), nrow = 3, ncol = 3)
> mxaa1
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 3 4
[3,] 2 5 6
> mxbb1
[,1] [,2] [,3]
[1,] 1 0 1
[2,] 0 1 0
[3,] 0 1 1
> mxaa1 %*% mxbb1
[,1] [,2] [,3]
[1,] 1 5 4
[2,] 2 7 6
[3,] 2 11 8
> mxaa2 <- matrix(c(1, 0, 0, 1, 1, 0), nrow = 2, ncol = 3)
> mxbb2 <- matrix(c(1, 0, 0, 2, 1, -2), nrow = 3, ncol = 2)
> mxaa2
[,1] [,2] [,3]
[1,] 1 0 1
[2,] 0 1 0
> mxbb2
[,1] [,2]
[1,] 1 2
[2,] 0 1
[3,] 0 -2
> mxaa2 %*% mxbb2
[,1] [,2]
[1,] 1 0
[2,] 0 1
> mxaa1 %*% mxaa2
mxaa1 %*% mxaa2 でエラー: 適切な引数ではありません行列のスカラー倍を求める
*演算子を使う。それぞれスカラー（数や関数のこと）と行列を指定する。
> a <- 3
> mxxx <- matrix(c(1, 2, 2, 2, 3, 5, 3, 4, 6), nrow = 3, ncol = 3)
> a
[1] 3
> mxxx
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 3 4
[3,] 2 5 6
> a * mxxx
[,1] [,2] [,3]
[1,] 3 6 9
[2,] 6 9 12
[3,] 6 15 18行列のアダマール積を求める
*演算子を使う。それぞれ行列を指定する。アダマール積の定義（行の数と列の数が同じ行列同士の各成分を掛ける）から、二つの行列の行数と列数が一致しない場合はエラーが発生する。
> 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 でエラー: 適切な配列ではありません線形代数
転置行列を求める
t関数を使う。
> mx <- matrix(1:9, 3, 3, byrow = TRUE)
> print(mx)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
> t(mx)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> mx <- matrix(1:6, 2, 3, byrow = TRUE)
> print(mx)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> t(mx)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6単位行列を作る
diag関数を使う。引数を一つだけ指定した場合は、n次の単位行列を作成する。
> diag(3)
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> diag(4)
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1対角行列を作成する
diag関数を使う。第1引数に対角要素の値、第2引数に求める正方行列の次数を指定する。
diag(4, 2)
[,1] [,2][1,] 4 0 [2,] 0 4
diag(5, 3)
[,1] [,2] [,3][1,] 5 0 0 [2,] 0 5 0 [3,] 0 0 5
diag(c(2, 4), 4)
[,1] [,2] [,3] [,4][1,] 2 0 0 0 [2,] 0 4 0 0 [3,] 0 0 2 0 [4,] 0 0 0 4 }}
上三角行列を作成する
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下三角行列を作成する
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には行列Aのトレース（固有和）Tr(A)を求めるための関数が標準では組み込まれていない。diag関数とsum関数を組み合わせることで、計算を実現することができる。以下の例では、Tr(A) = 1+5+9=15である。
> mx <- matrix(1:9, nrow = 3, byrow = TRUE)
> mx
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
> tr <- sum(diag(mx))
> tr
[1] 15行列式を求める
det関数を使う。
> d <- c(1, 2, 3, 4)
> mx <- matrix(d, 2, 2, byrow = TRUE)
> mx
[,1] [,2]
[1,] 1 2
[2,] 3 4
> det(mx)
[1] -2行列式を求める
det関数を使う。
> mx <- matrix(3, 1, 1)
> print(mx)
[,1]
[1,] 3
> det(mx)
[1] 3
> mx <- matrix(1:4, 2, 2, byrow = TRUE)
> print(mx)
[,1] [,2]
[1,] 1 2
[2,] 3 4
> det(mx)
[1] -2
> mx <- matrix(c(0, -4, 0, -1, 3, 6, -2, 0, 2, 1, 3, 0, -1, 5, -1, 2), 4, 4)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 0 3 2 -1
[2,] -4 6 1 5
[3,] 0 -2 3 -1
[4,] -1 0 0 2
> det(mx)
[1] 28逆行列を求める
solve関数を使う。
> mx <- matrix(3, 1, 1)
> print(mx)
[,1]
[1,] 3
> solve(mx)
[,1]
[1,] 0.3333333
> mx <- matrix(1:4, 2, 2, byrow = TRUE)
> print(mx)
[,1] [,2]
[1,] 1 2
[2,] 3 4
> solve(mx)
[,1] [,2]
[1,] -2.0 1.0
[2,] 1.5 -0.5
> mx <- matrix(c(3, 3, -1, -3, 2, -5, 1, 0, 1), 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 3 -3 1
[2,] 3 2 0
[3,] -1 -5 1
> solve(mx)
[,1] [,2] [,3]
[1,] 1.0 -1 -1.0
[2,] -1.5 2 1.5
[3,] -6.5 9 7.5正則ではない行列（行列式の値が0）の場合は、solve関数はエラーを返す。
> mx <- matrix(1:9, 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> det(mx)
[1] 0
> solve(mx)
solve.default(mx) でエラー:
Lapack routine dgesv: システムは正確に特異です: U[3,3] = 0一般逆行列を求める
MASSパッケージのginv関数を使う。
> library(MASS)
> mx <- matrix(1:16, 4, 4, byrow = TRUE)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
[4,] 13 14 15 16
> det(mx)
[1] 4.733165e-30
> solve(mx)
solve.default(mx) でエラー:
システムは数値的に特異です: 条件数の逆数 = 2.95756e-18
> ginv(mx)
[,1] [,2] [,3] [,4]
[1,] -0.2850 -0.1450 -0.0050 0.1350
[2,] -0.1075 -0.0525 0.0025 0.0575
[3,] 0.0700 0.0400 0.0100 -0.0200
[4,] 0.2475 0.1325 0.0175 -0.0975一般逆行列の定義を満たしているか否か確認。
> mxi <- ginv(mx)
> mx %*% mxi %*% mx
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
[4,] 13 14 15 16当然、正方行列でなくとも求めることができる。
> mx <- matrix(1:8, 2, 4, byrow = TRUE)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
> det(mx)
determinant.matrix(x, logarithm = TRUE, ...) でエラー:
'x' は正方行列でなければなりません
> solve(mx)
solve.default(mx) でエラー: 'a' (2 x 4) は正方行列である必要があります
> ginv(mx)
[,1] [,2]
[1,] -0.550 2.500000e-01
[2,] -0.225 1.250000e-01
[3,] 0.100 -2.081668e-17
[4,] 0.425 -1.250000e-01
> mxi <- ginv(mx)
> mx %*% mxi %*% mx
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8ginv関数はムーア・ペンローズ一般逆行列を数値的に求める関数のため、求まった一般逆行列は近似値であることに注意。ムーア・ペンローズ一般逆行列の定義の4つの式について、比較をした例は以下のとおり。すべての要素がTRUEにならないことがわかる。
> mx %*% mxi %*% mx == mx
[,1] [,2] [,3] [,4]
[1,] FALSE FALSE TRUE TRUE
[2,] FALSE FALSE FALSE TRUE
> mxi %*% mx %*% mxi == mxi
[,1] [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,] FALSE FALSE
[4,] FALSE FALSE
> t(mx %*% mxi) == mx %*% mxi
[,1] [,2]
[1,] TRUE FALSE
[2,] FALSE TRUE
> t(mxi %*% mx) == mxi %*% mx
[,1] [,2] [,3] [,4]
[1,] TRUE FALSE FALSE FALSE
[2,] FALSE TRUE FALSE FALSE
[3,] FALSE FALSE TRUE FALSE
[4,] FALSE FALSE FALSE TRUE連立一次方程式を解く
solve関数を使う。以下の連立一次方程式を解く（解はx=3,y=2）。
x + 2y = 7
3x - 4y = 1> mxaa <- matrix(c(1, 2, 3, -4), 2, 2, byrow = TRUE)
> mxy <- matrix(c(7, 1), 2, 1)
> print(mxaa)
[,1] [,2]
[1,] 1 2
[2,] 3 -4
> print(mxy)
[,1]
[1,] 7
[2,] 1
> solve(mxaa, mxy)
[,1]
[1,] 3
[2,] 2同じく以下の連立一次方程式を解く（解はx=2,y=-3,z=-14）。
3x - 3y + z = 1
3x + 2y = 0
-1x - 5y + z = -1> mxaa <- matrix(c(3, 3, -1, -3, 2, -5, 1, 0, 1), 3, 3)
> mxy <- matrix(c(1, 0, -1), 3, 1)
> print(mxaa)
[,1] [,2] [,3]
[1,] 3 -3 1
[2,] 3 2 0
[3,] -1 -5 1
> print(mxy)
[,1]
[1,] 1
[2,] 0
[3,] -1
> solve(mxaa, mxy)
[,1]
[1,] 2
[2,] -3
[3,] -14計算
QR分解をする
行列XをQR分解で行列Qと行列Rに分解することを考える。
X = Q R行列QはQ^T Q=I（Iは単位行列）となる直交行列。行列Rは上三角行列（対角成分はすべて0以外で、対角成分より下の要素はすべて0）。
qr.Q関数、qr.R関数を使うと、それぞれQR分解した際の行列Qと行列Rが得られる。引数にはqrオブジェクトを与える必要があり、qr関数を使ってあらかじめ分解したい行列Xからqrオブジェクトを作成しておく必要がある。
> d <- c(1, 2, 3, 4, 5, 6, 9, 9, 8)
> x <- matrix(d, 3, 3, byrow = TRUE)
> x
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 9 9 8
> objqr <- qr(x)
> q <- qr.Q(objqr)
> q
[,1] [,2] [,3]
[1,] -0.1010153 -0.7184092 -0.6882472
[2,] -0.4040610 -0.6025367 0.6882472
[3,] -0.9091373 0.3476173 -0.2294157
> r <- qr.R(objqr)
> r
[,1] [,2] [,3]
[1,] -9.899495 -10.404571 -10.0005102
[2,] 0.000000 -1.320946 -2.9895090
[3,] 0.000000 0.000000 0.2294157Q^T Q = I となるか試す。
> t(q) %*% q
[,1] [,2] [,3]
[1,] 1.000000e+00 -1.110223e-16 -1.665335e-16
[2,] -1.110223e-16 1.000000e+00 -9.714451e-17
[3,] -1.665335e-16 -9.714451e-17 1.000000e+00qr.X関数を使うことで、qrオブジェクトの元の行列を取り出すことができる。
> qr.X(objqr)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 9 9 8コレスキー分解する
chol関数を使う。
> mx0 <- matrix(c(1, 0, 1, 0, 2, 0, 1, 0, 3), 3, byrow = TRUE)
> mx0
[,1] [,2] [,3]
[1,] 1 0 1
[2,] 0 2 0
[3,] 1 0 3
> mx <- chol(mx0)
> mx
[,1] [,2] [,3]
[1,] 1 0.000000 1.000000
[2,] 0 1.414214 0.000000
[3,] 0 0.000000 1.414214得られた結果が正しいか、確認する。
> t(mx) %*% mx
[,1] [,2] [,3]
[1,] 1 0 1
[2,] 0 2 0
[3,] 1 0 3上記と同じ計算（A T A）はcrossprod関数を使うと同じ結果を得ることができる。
> crossprod(mx)
[,1] [,2] [,3]
[1,] 1 0 1
[2,] 0 2 0
[3,] 1 0 3固有値と固有ベクトルを求める
eigen関数を使う。
> # 例1
> d <- c(1, 3, -2, -4)
> mx <- matrix(d, 2, 2)
> print(mx)
[,1] [,2]
[1,] 1 -2
[2,] 3 -4
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] -2 -1
$vectors
[,1] [,2]
[1,] 0.5547002 0.7071068
[2,] 0.8320503 0.7071068
> # 固有値（合計2つ）
> print(eig$value)
[1] -2 -1
> # 固有値 -2 のときの固有ベクトル
> print(eig$vector[, 1])
[1] 0.5547002 0.8320503
> # 固有値 -1 のときの固有ベクトル
> print(eig$vector[, 2])
[1] 0.7071068 0.7071068
>
> # 例2
> d <- c(1, 1, 1, 0, 1, 0, 2, 1, 0)
> mx <- matrix(d, 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 1 0 2
[2,] 1 1 1
[3,] 1 0 0
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] 2 1 -1
$vectors
[,1] [,2] [,3]
[1,] 0.5345225 0 -0.7071068
[2,] 0.8017837 1 0.0000000
[3,] 0.2672612 0 0.7071068
> # 固有値（合計3つ）
> print(eig$value)
[1] 2 1 -1
> # 固有値 2 のときの固有ベクトル
> print(eig$vector[, 1])
[1] 0.5345225 0.8017837 0.2672612
> # 固有値 1 のときの固有ベクトル
> print(eig$vector[, 2])
[1] 0 1 0
> # 固有値 -1 のときの固有ベクトル
> print(eig$vector[, 3])
[1] -0.7071068 0.0000000 0.7071068固有値は数値的に求めているため、基本的に戻り値は近似値である。そのため、解析的には整数で得られる場合でも戻り値が複素数になる場合がある。この場合、固有ベクトルも複素数になる。以下の例では、解析的には固有値は3と-1（二重根）と求められるが、戻り値は複素数になっている。虚部が無視できるくらい小さいため、虚部を切り捨てて実数ととして扱ってかまわないが、ケースバイケースであることに注意。
> d <- c(1, 2, -1, 2, 1, -1, -4, 4, -1)
> mx <- matrix(d, 3, 3)
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] 3+0i -1+0i -1-0i
$vectors
[,1] [,2] [,3]
[1,] -0.9045340+0i -7.071068e-01+0.000000e+00i -7.071068e-01+0.000000e+00i
[2,] -0.3015113+0i 7.071068e-01-0.000000e+00i 7.071068e-01+0.000000e+00i
[3,] 0.3015113+0i 0.000000e+00+2.637113e-09i 0.000000e+00-2.637113e-09i固有値と固有ベクトルを求める
eigen関数を使う。戻り値はリストで、順番にvaluesに固有値が、vectorsにその固有値のときの固有ベクトルが含まれる。
> mxaa <- matrix(c(3, 1, 1, 1, 2, 0, 1, 0, 2), nrow = 3)
> print(mxaa)
[,1] [,2] [,3]
[1,] 3 1 1
[2,] 1 2 0
[3,] 1 0 2
> lam <- eigen(mxaa)
> print(lam)
eigen() decomposition
$values
[1] 4 2 1
$vectors
[,1] [,2] [,3]
[1,] 0.8164966 0.0000000 0.5773503
[2,] 0.4082483 -0.7071068 -0.5773503
[3,] 0.4082483 0.7071068 -0.5773503この例では固有値と固有ベクトルは3つあり、例えば3番目の固有値と固有ベクトルを取り出すには、以下のようにする。
> lam$values[3]
[1] 1
> lam$vectors[, 3]
[1] 0.5773503 -0.5773503 -0.5773503特異値分解を行う
svd関数を使う。特異値分解の等式は以下のとおり。
A = U D V^TUとVは直交行列。DはAの特異値を対角成分にもつ対角行列。
> a <- matrix(1:9, 3, 3)
> a
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> r <- svd(a)
> r
$d
[1] 1.684810e+01 1.068370e+00 5.543107e-16
$u
[,1] [,2] [,3]
[1,] -0.4796712 0.77669099 0.4082483
[2,] -0.5723678 0.07568647 -0.8164966
[3,] -0.6650644 -0.62531805 0.4082483
$v
[,1] [,2] [,3]
[1,] -0.2148372 -0.8872307 0.4082483
[2,] -0.5205874 -0.2496440 -0.8164966
[3,] -0.8263375 0.3879428 0.4082483svd関数の戻り値はリスト。要素dに特異値が代入されている。
> r$d
[1] 1.684810e+01 1.068370e+00 5.543107e-16上記の恒等式における、行列U（戻り値では要素u）と行列V（戻り値では要素v）が直交行列であることを確認する。
> t(r$u) %*% r$u
[,1] [,2] [,3]
[1,] 1.000000e+00 1.110223e-16 5.551115e-17
[2,] 1.110223e-16 1.000000e+00 -1.110223e-16
[3,] 5.551115e-17 -1.110223e-16 1.000000e+00
> t(r$v) %*% r$v
[,1] [,2] [,3]
[1,] 1.000000e+00 0.000000e+00 -1.110223e-16
[2,] 0.000000e+00 1.000000e+00 1.110223e-16
[3,] -1.110223e-16 1.110223e-16 1.000000e+00上記の恒等式が成り立っているか確認する。
> r$u %*% diag(r$d) %*% t(r$v)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9行列Aの特異値は、A^T Aの正の固有値の正の平方根でも求められる。
> sqrt(eigen(t(a) %*% a)$values)
[1] 1.684810e+01 1.068370e+00 7.942757e-083つ目は実質0ということだが、高精度に計算されていないことがわかる。
余因子を求める
polyMatrixパッケージのcofactor関数を使う。以下、使用例。計算に用いた行列はここを参考にしており、先に余因子行列をadjoint関数を使用して求めている。
> d <- c(8, 3, 4, 1, 5, 9, 6, 7, 2)
> mx <- matrix(d, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 8 1 6
[2,] 3 5 7
[3,] 4 9 2
> adjoint(mx)
[,1] [,2] [,3]
[1,] -53 52 -23
[2,] 22 -8 -38
[3,] 7 -68 37
> cofactor(mx, 1, 2)
[1] 22
> cofactor(mx, 2, 3)
[1] -68余因子行列を求める
polyMatrixパッケージのadjoint関数を使う。以下、使用例。計算に用いた行列はここhttps://jp.mathworks.com/help/symbolic/adjoint.htmlを参考にしている。
> d <- c(8, 3, 4, 1, 5, 9, 6, 7, 2)
> mx <- matrix(d, 3)
> mx
[,1] [,2] [,3]
[1,] 8 1 6
[2,] 3 5 7
[3,] 4 9 2
> adjoint(mx)
[,1] [,2] [,3]
[1,] -53 52 -23
[2,] 22 -8 -38
[3,] 7 -68 37余因子行列の性質 A A~ = A~ A ~ = |A|・En を確認する。
> print(mx %*% adjoint(mx))
[,1] [,2] [,3]
[1,] -3.600000e+02 2.842171e-13 -5.684342e-14
[2,] -7.815970e-14 -3.600000e+02 -5.684342e-14
[3,] -8.704149e-14 1.421085e-13 -3.600000e+02
> print(adjoint(mx) %*% mx)
[,1] [,2] [,3]
[1,] -3.600000e+02 2.842171e-14 -1.136868e-13
[2,] 5.684342e-14 -3.600000e+02 5.684342e-14
[3,] 8.526513e-14 1.705303e-13 -3.600000e+02
> print(det(mx) * diag(nrow(mx)))
[,1] [,2] [,3]
[1,] -360 0 0
[2,] 0 -360 0
[3,] 0 0 -360上2つの計算結果の非対角成分は非常に小さいので0と見なしてよく、上記の性質を満たしていることが確認できた。
ヘッセ行列を求める
numDerivパッケージのhessian関数を使う。第一引数に目的関数を、第二引数に変数の値をベクトルで与える。
> library(numDeriv)
> f1 <- function(x) x[1] ^ 2 + 4 * x[2] ^ 2
> hessian(f1, c(0, 0))
[,1] [,2]
[1,] 2.000000e+00 -3.308722e-16
[2,] -3.308722e-16 8.000000e+00
> f2 <- function(x) x[1] ^ 2 + 2 * x[2] ^ 2
> hessian(f2, c(0, 0))
[,1] [,2]
[1,] 2.000000e+00 3.308722e-16
[2,] 3.308722e-16 4.000000e+00
> f3 <- function(x) 2 * x[1] - 4 * x[2] + x[1] ^ 2 + 2 * x[2] ^ 2 + 2 * x[1] * x[2]
> hessian(f3, c(0, 0))
[,1] [,2]
[1,] 2 2
[2,] 2 4勾配ベクトルを求める
numDerivパッケージのjacobian関数を使う。第一引数に目的関数を、第二引数に変数の値をベクトルで与える。
> library(numDeriv)
> f1 <- function(x) {6 - 0.1 * (x[1] - 5) ^ 2} * {6 - 0.04 * (x[2] - 5) ^ 2}
> jacobian(f1, c(8.25, 5.25))
[,1] [,2]
[1,] -3.898375 -0.098875
> f2 <- function(x) 2 * x[1] - 4 * x[2] + x[1] ^ 2 + 2 * x[2] ^ 2 + 2 * x[1] * x[2]
> jacobian(f2, c(1, 1))
[,1] [,2]
[1,] 6 2その他
行列の要素を抜き出してベクトルにする
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パスカル行列を求める
パスカル行列は正定値対称行列であり、下三角コレスキー因子（パスカル行列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 70LU分解を行う
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 20LU分解を行う
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 20QR分解を行う
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 330qr関数に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後退代入法で連立方程式を解く
backsolve関数を使う。以下は連立方程式を行列で表記し、上三角行列の係数行列をA、定数項行列をy、求める解をxとおき、x = Ayと考えてxを求めた例。backsolve関数だけではなく、solve関数を使った例、逆行列を求めて正規方程式を解いた結果も示したが、得られた解は同じであることがわかる。
> set.seed(6)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxaa[lower.tri(mxaa)] <- 0
> print(mxaa)
[,1] [,2] [,3] [,4]
[1,] 6 8 5 1
[2,] 0 9 1 3
[3,] 0 0 6 7
[4,] 0 0 0 3
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxy)
[,1]
[1,] 5
[2,] 7
[3,] 2
[4,] 7
> backsolve(mxaa, mxy)
[,1]
[1,] 2.0812757
[2,] 0.2654321
[3,] -2.3888889
[4,] 2.3333333
> solve(mxaa, mxy)
[,1]
[1,] 2.0812757
[2,] 0.2654321
[3,] -2.3888889
[4,] 2.3333333
> print(solve(t(mxaa) %*% mxaa) %*% (t(mxaa) %*% mxy))
[,1]
[1,] 2.0812757
[2,] 0.2654321
[3,] -2.3888889
[4,] 2.3333333backsolve関数は、与えられた行列の上三角の成分だけを使って解く。以下のとおり、すべての成分を使って解くsolve関数とは得られた解が違うことがわかる。
> set.seed(6)
> 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,] 6 8 5 1
[2,] 9 9 1 3
[3,] 3 9 6 7
[4,] 4 7 9 3
> print(mxy)
[,1]
[1,] 5
[2,] 7
[3,] 2
[4,] 7
> backsolve(mxaa, mxy)
[,1]
[1,] 2.0812757
[2,] 0.2654321
[3,] -2.3888889
[4,] 2.3333333
> solve(mxaa, mxy)
[,1]
[1,] 1.8801598
[2,] -1.3715047
[3,] 0.8322237
[4,] 0.5299601相関行列を求める
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前進代入法で連立方程式を解く
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.5416667forwardsolve関数は、与えられた行列の下三角の成分だけを使って解く。以下のとおり、すべての成分を使って解く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零行列を作成する
すべての成分が0である行列を零行列という。Rには零行列を作成する関数はないが、行列を作成するmatrix関数を使えば、簡単に作成することができる。第一引数に要素が一つだけのベクトルを指定すれば、すべての成分が値である行列が作成されるので、この機能を使えばよい。
> matrix(0.0, nrow = 2, ncol = 3)
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
> matrix(0.0, 2, 4)
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 0 0 0 0