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行列から成分を取り出す
行ごと列ごとに角括弧（ブラケット）をつかってインデックス（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 番目の引数を参照)線形代数
転置行列を求める
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 }}
行列のトレース（固有和）を求める
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行列の掛け算
%*% 演算子を使う。最後の例のとおり、掛けられる行列の列数と掛ける行列の行数が同じでなければ計算することはできない。
> 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 でエラー: 適切な引数ではありません逆行列を求める
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
