行列演算
今回は行列(matrix)の扱いについて書きます。行列オブジェクトをどう定義するか、積の求め方、逆行列をどうやって求めるか、固有値・固有ベクトルをどうやって求めるかについて解説します。途中、行列の転置の仕方、det の求め方、アダマール積の求め方についても触れます。
行列を定義
matrix()の引数に値、行数、列数を入れて作成する
> A <- matrix(1:9, 3, 3)
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
このように1:9で１〜9の連続した整数を指定して3行３列の行列を埋めます。
ただし、デフォルトでは列優先になってしまうので、それが嫌なら転置するかbyrow=TRUE を指定してください。(下記@WolfMoon さんのコメント参照)
> A <- t(A) # 転置
> A
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
> A <- matrix(1:9, 3, 3, byrow=TRUE)
> A
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
これで行優先(すなわち第1行に1〜3、第２行に4〜6... )になります。
まったくの余談ですが、R の行列オブジェクトでは、行や列に名前をつけることができます。
> colnames(A) <- c("a", "b", "c") # 列名の変更
> rownames(A) <- c("r1", "r2", "r3") # 行名の変更
> A
a b c
r1 1 2 3
r2 4 5 6
r3 7 8 9
例えば、次のような定義の仕方もありです。
> C <- matrix(c(1:3, 10:12, 15:17), 3, 3)
> C
[,1] [,2] [,3]
[1,] 1 10 15
[2,] 2 11 16
[3,] 3 12 17
> D <- matrix(rep(1, 9), 3, 3)
> D
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 1
[3,] 1 1 1
> E <- matrix(c(rep(1, 5), rep(2, 4)), 3, 3)
> E
[,1] [,2] [,3]
[1,] 1 1 2
[2,] 1 1 2
[3,] 1 2 2
くれぐれも(デフォルトでは)列優先であることに注意してください。
その他、まず最初に空の行列を定義しておいて、後から数値を代入していくという方法もよく用いられます。
> A <- matrix(NA, 3, 3)
> A
[,1] [,2] [,3]
[1,] NA NA NA
[2,] NA NA NA
[3,] NA NA NA
これで「空っぽ」の$3\times 3$行列が定義されるので、値を代入してきます。
> A[1, ] <- 1:3
> A[2, ] <- 4:6
> A[3, ] <- 7:9
> A
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
積
上記のmatrix Aと次で定義されるmatrix Bを用います。
> B <- matrix(1:6, 3, 2)
> B
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
これらの行列の積$AB$は、以下で計算されます。
> A %*% B
[,1] [,2]
[1,] 14 32
[2,] 32 77
[3,] 50 122
積が定義できない場合、以下のようなエラーが発生します。
> B %*% A
B %*% A でエラー: 適切な引数ではありません
以上のように、いわゆる「行列の積」を行う演算子は%*%です。
*は「行列の積」ではなく、アダマール積(同じ位置の要素同士の積)を求める演算子ですので、注意してください。
> A * C # アダマール積
[,1] [,2] [,3]
[1,] 1 20 45
[2,] 8 55 96
[3,] 21 96 153
逆行列
まず、逆行列が存在するか、$\mathrm{det}(A)$ を調べましょう。
> det(A)
[1] 6.661338e-16
Aはダメっぽいです。代わりに以下の行列Fを用います。
> F
[,1] [,2] [,3]
[1,] 1 2 9
[2,] 4 5 5
[3,] 7 8 2
Fの逆行列を求めるには、solve()を用います。
> solve(F)
[,1] [,2] [,3]
[1,] 10 -22.66667 11.66667
[2,] -9 20.33333 -10.33333
[3,] 1 -2.00000 1.00000
ちょっと桁が長いですね。例えば以下のようにすれば、短く表示することができます。
> library(tidyverse)
> solve(F) %>% print(digits=3)
[,1] [,2] [,3]
[1,] 10 -22.7 11.7
[2,] -9 20.3 -10.3
[3,] 1 -2.0 1.0
%>% はライブラリtidyverse で定義されるパイプ・オペレータです。
固有値・固有ベクトル
ここでは以下のmatrix G を用います。
> G
[,1] [,2] [,3]
[1,] 0 -1 0
[2,] -1 0 0
[3,] 0 2 4
固有値・固有ベクトルを求めるにはeigen() を用います。
> eigen(G)
eigen() decomposition
$values
[1] 4 1 -1
$vectors
[,1] [,2] [,3]
[1,] 0 -0.6396021 0.6804138
[2,] 0 0.6396021 0.6804138
[3,] 1 -0.4264014 -0.2721655
結果の読み方ですが、まず固有値として$$\lambda = 4, 1, -1$$の3つが求まりました。各固有値と固有ベクトルの対応は以下の通りです。
$\lambda=4$に属する固有ベクトルは
\left(\begin{array}{c}
0\\
0\\
-1
\end{array}\right)
$\lambda=1$に属する固有ベクトルは
\left(\begin{array}{r}
-0.6396021 \\
0.6396021 \\
-0.4264014 \\
\end{array}\right)
$\lambda=-1$に属する固有ベクトルは
\left(\begin{array}{r}
0.6804138 \\
0.6804138 \\
-0.2721655 \\
\end{array}\right)
となります。つまり、eigen(G)$vectors は3つの列ベクトル(縦ベクトル)の集まりで、左から順に、固有値$\lambda=4, 1, -1$に属する固有ベクトルです。
参考文献
トップページはこちら