Rで「あてはめ」ができるまでの道筋　行列演算

# 行列演算
今回は行列(matrix)の扱いについて書きます。行列オブジェクトをどう定義するか、積の求め方、逆行列をどうやって求めるか、固有値・固有ベクトルをどうやって求めるかについて解説します。途中、行列の転置の仕方、det の求め方、アダマール積の求め方についても触れます。

## 行列を定義
matrix()の引数に値、行数、列数を入れて作成する
```R
> 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 さんのコメント参照)
```R
> 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 の行列オブジェクトでは、行や列に名前をつけることができます。
```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
```

例えば、次のような定義の仕方もありです。
```R
> 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
```
くれぐれも(デフォルトでは)列優先であることに注意してください。

その他、まず最初に空の行列を定義しておいて、後から数値を代入していくという方法もよく用いられます。
```R
> A <- matrix(NA, 3, 3)
> A
     [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]   NA   NA   NA
[3,]   NA   NA   NA
```
これで「空っぽ」の$3\times 3$行列が定義されるので、値を代入してきます。
```R
> 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を用います。
```R
> B <- matrix(1:6, 3, 2)
> B
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
```

これらの行列の積$AB$は、以下で計算されます。
```R
> A %*% B
     [,1] [,2]
[1,]   14   32
[2,]   32   77
[3,]   50  122
```
積が定義できない場合、以下のようなエラーが発生します。
```R
> B %*% A
 B %*% A でエラー:  適切な引数ではありません 
```

以上のように、いわゆる「行列の積」を行う演算子は`%*%`です。
`*`は「行列の積」ではなく、アダマール積(同じ位置の要素同士の積)を求める演算子ですので、注意してください。

```R
> A * C # アダマール積
     [,1] [,2] [,3]
[1,]    1   20   45
[2,]    8   55   96
[3,]   21   96  153
```

## 逆行列
まず、逆行列が存在するか、$\mathrm{det}(A)$ を調べましょう。
```R
> det(A)
[1] 6.661338e-16
```
Aはダメっぽいです。代わりに以下の行列Fを用います。
```R
> F
     [,1] [,2] [,3]
[1,]    1    2    9
[2,]    4    5    5
[3,]    7    8    2
```
Fの逆行列を求めるには、solve()を用います。
```R
> solve(F)
     [,1]      [,2]      [,3]
[1,]   10 -22.66667  11.66667
[2,]   -9  20.33333 -10.33333
[3,]    1  -2.00000   1.00000
```

ちょっと桁が長いですね。例えば以下のようにすれば、短く表示することができます。
```R
> 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 を用います。

```R
> G
     [,1] [,2] [,3]
[1,]    0   -1    0
[2,]   -1    0    0
[3,]    0    2    4
```
固有値・固有ベクトルを求めるにはeigen() を用います。
```R
> 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$に属する固有ベクトルは
```math
\left(\begin{array}{c}
0\\
0\\
-1
\end{array}\right)
```
$\lambda=1$に属する固有ベクトルは
```math
\left(\begin{array}{r}
-0.6396021 \\
0.6396021 \\
-0.4264014 \\
\end{array}\right)
```

$\lambda=-1$に属する固有ベクトルは
```math
\left(\begin{array}{r}
0.6804138 \\
0.6804138 \\
-0.2721655 \\
\end{array}\right)
```
となります。つまり、`eigen(G)$vectors` は3つの列ベクトル(縦ベクトル)の集まりで、左から順に、固有値$\lambda=4, 1, -1$に属する固有ベクトルです。
## 参考文献
[biostatistics](https://stats.biopapyrus.jp/r/basic/matrix.html)

トップページは[こちら](https://qiita.com/sflabtech/items/562dd78418926f00a0d5)