2  R を用いた行列の計算

この教材では今後、分散共分散行列など行列を用いた説明も出てくる。 そこで、ここでは線形代数の基本的な事柄として、 行列、固有値、固有ベクトル、行列の対角化について述べる。 また、Rにおける行列計算の方法について述べる。

キーワード

行列、固有値、固有ベクトル、対角化、正定値行列

2.1 記述統計量

\(n\)人の身長など１種類のデータがあるとする。 その値 を\(x_1\)、\(x_2\)、\(\cdots\)、\(x_n\)とする。 通常、観測されるデータは母集団全てのデータが得られるわけではなく、 その一部である。得られたデータを標本（サンプル）という。 この標本をもとに母集団の性質を知りたいと考える。 データ数 \(n\) のことをサンプルサイズという。 母集団の性質を知るために標本データから計算される量を 統計量という。基本的な統計量として 最小値、最大値、中央値、平均値、分散などがある。 Chapter 5 で述べるように、観測されたデータを分析する場合には その背景にある確率分布を想定することが多い。 大文字で表した時にはある確率分布の特徴を推定するための 推定する式としての意味がある。 また、単に観測値から計算だけを説明したい時には小文字を使う。 この章では小文字で使い、分布を踏まえて議論するときには 大文字を使う。

Table 2.1: 基本統計量とよく用いられる文字

統計量	変数
平均	\(\bar{x}\)、\(\mu_x\)
期待値	\(E(X)\)
標準偏差	\(\sigma_x\)
分散	\(V(X)\)、\(\sigma_x^2\)
偏差二乗和	\(S_{xx}\)
共分散	\({\rm Cov}(x,y)\)、\(\sigma_{xy}\)
相関係数	\({\rm Cor}(x,y)\)、\(\rho_{xy}\)

データの総和をサンプルサイズ\(n\) で割ったものを平均（mean または average）という。

\[\begin{eqnarray} \mu_x &=& \frac{1}{n} \sum_{i=1}^{n}x_{i} \end{eqnarray}\]

\((x_{i}-\mu_{x})\)を偏差という。偏差の総和は\(0\)となる。 つまり、この\(\bar{x}\)を基準に値を見直すと正のものと負のものに 別れて、全て足すと\(0\) になってバランスが取れる。つまり重心に なっていることがわかる。 また、偏差を二乗したものの総和を偏差二乗和 と言い、\(S_{xx}\)とする。

\[\begin{eqnarray} S_{xx} &=& {\sum_{i=1}^{n}(x_{i}-\bar{x})^2} \end{eqnarray}\]

である。右辺を展開すると

\[\begin{eqnarray} S_{xx} &=& \sum_{i=1}^{n}x_{i}^2 -n \bar{x}^2 \end{eqnarray}\]

となる。そのばらつき具合を表すものを分散（variance）という。 分散はこれを \(n\) または \(n-1\) で割った値が用いられる。 R で分散を計算する場合には \(n-1\)で割った値を用いることも多い。 これを不偏分散という。

\[\begin{eqnarray} \sigma_{x}^2 &=& \frac{S_{xx}}{n-1} \end{eqnarray}\]

である。また、分散の\(1/2\)乗した統計量 標準偏差(standard deviation) \(\sigma_x\) という。 また、同じ人数\(n\)人のデータ \(y_1\)、\(y_2\)、\(\cdots\)、\(y_n\) があるとする。このとき、両者の関係を調べる統計量として 共分散がある。

\[\begin{eqnarray} S_{xy}&=& \sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y}) \\ &=& \sum_{i=1}^{n}x_iy_i -n\bar{x}\bar{y} \\ {\rm Cov}_{xy} &=& \sigma_{xy}= \frac{S_{xy}}{n-1} \end{eqnarray}\]

を共分散という。共分散 共分散の値の大きさはそれぞれの項目のばらつきの大きさに依存する。 そこで、共分散を二つの標準偏差によって割ると、 \(-1\) から\(1\) までの範囲となり、 両者の関係の強さを測る指標となる。これを（積率）相関係数という。

\[\begin{eqnarray} {\rm Cor}(x,y) &= & \rho_{xy} = \frac{{\rm Cov}(x,y)} {\sigma_x \sigma_y} \end{eqnarray}\]

相関係数は正の値だけでなく負の値もとる。値が \(1\)や \(-1\) に近いほど「相関が高い」、値が\(0\) に近いほど「相関がない」という。

A	B
118cm	128cm
119cm	129cm
121cm	130cm
122cm	131cm
170cm	132cm

グループA、Bの身長 {#tbl-base1}

この身長の平均値、中央値、分散、標準偏差をそれぞれRを使って表してみよう。 R では 平均は「mean()」、中央値は「median()」、分散は「var()」、 標準偏差は「sd()」で求めることができる。 複数の要素を持つデータ（ベクトルデータ）を一つのまとまりとして表現するには 「c()」(combine、または concatenate)を使う。まず、5人分のデータをそれぞれ x1、x2 として表そう。 この値をもとに R で計算するとつぎのようにすぐに結果が得られる。 この「mean」や「var」のように決められた処理を行なって結果を返す一連の命令群のことを関数という。 関数はある値やデータを入力して結果を返す。そこで今後は、どういったデータを入力するとどういった結果が得られ、それをどのように解釈するか、どういうことに注意するのか、といったことに着目して説明する。

> x1 <- c(118,119,121,122,170)
> x2 <- c(128,129,130,131,132)
> mean(x1)

[1] 130

> var(x1)

[1] 502.5

> cor(x1,x2)

[1] 0.7547198

Table 2.2: 主な Rの統計量関数

統計量	説明
sum	合計、例　sum(a,b,c)
mean	平均値、中央値は median
max	最大値、最小値は min
var	分散、行列を与えると分散共分散行列を計算する
sd	標準偏差、それぞれの列の標準偏差を求める
cor	相関、行列を与えると相関行列を求める

2.2 行列

\(n\) 個の変数があるとき、これを次のように \(m\) 個の線形に変換することを考える。この講義では各成分は 実数であるものとする。

\[\begin{eqnarray*} \label{eqn:matrix} y_1 &=& a_{11}x_1+a_{12}x_2 +\cdots+ a_{1n}x_n \\ y_2 &=& a_{21}x_1+a_{22}x_2 +\cdots+ a_{2n}x_n \\ \nonumber \\ &&\vdots \\ y_m &=& a_{m1}x_1+a_{m2}x_2 +\cdots+ a_{mn}x_n \end{eqnarray*}\]

これを

\[\begin{eqnarray} \label{eqn:matrix1} \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_m \end{array} \right) &=& \left( \begin{array}{cccc} a_{11}&a_{12}&\cdots&a_{1n} \\ a_{21}&a_{22}&\cdots&a_{2n} \\ \vdots&\vdots&\ddots &\vdots \\ a_{m1}&a_{m2}&\cdots&a_{mn} \\ \end{array} \right ) \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right) \end{eqnarray}\]

のように表す。 また、

\[\begin{eqnarray} {\mathbf y} = \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_m \end{array} \right) \; ,\; & A = \left( \begin{array}{cccc} a_{11}&a_{12}&\cdots&a_{1n} \\ a_{21}&a_{22}&\cdots&a_{2n} \\ \vdots&\vdots&\ddots &\vdots \\ a_{m1}&a_{m2}&\cdots&a_{mn} \\ \end{array} \right ) \; &, \; {\mathbf x} = \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right) \end{eqnarray}\]

とすると、 \({\mathbf y}= A{\mathbf x}\) と表すことができる。行列を 列ベクトルで表したが、これを横に並べた行ベクトル で表現することもある。 列ベクトルを行ベクトルに 変換する操作のことを転置といい、 \({\mathbf x}^{T}\) と書く。行列も同様に

\[\begin{eqnarray} \label{eqn:trans} A^{T} &=& \left( \begin{array}{cccc} a_{11} & a_{21} & \cdots &a_{m1} \\ a_{12} & a_{22} & \cdots &a_{m2} \\ \vdots & \vdots & \ddots &\vdots \\ \vdots & \vdots & \ddots &\vdots \\ a_{1n} & a_{2n} & \cdots &a_{mn} \\ \end{array} \right ) \end{eqnarray}\]

のように縦と横の成分を入れ替えることを考えることができ、 この行列を転置行列という。転置行列は、\(n \times m\) 行列となる。 \(n\) と \(m\) の値が等しいとき、行列は正方形の形になる。この行列を正方行列という。正方行列において、対角成分\(a_{ij}=a_{ji}\) のとき、転置行列が元の行列と一致する。この行列のことを対称行列 という。

2.3 行列の演算

行列の和や差については、同じサイズの行列に対して、そのそれぞれの成分の和や差を計算する。行列の積については、\(A\) を \(l \times m\) 行列、\(B\) を \(m \times n\) とすると、行列の積 \(AB\) を計算することができ、 \(l \times n\) 行列 となる。 このとき、\(AB\) の \(ij\) 成分を \(c_{ij}\) とすると、\(A\) の \(i\) 行目の \(m\) 個の成分 \(a_{ik}\) と \(B\) の\(j\) 列目の\(m\) 個の成分 \(b_{kj}\) を掛け合わせたものを全て足すことで 計算される。

\[\begin{eqnarray} c_{ij} = a_{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{im}b_{mj} = \sum_{k=1}^{m} a_{ik} b_{kj}　 \end{eqnarray}\]

このように積が計算できるためには \(A\) の列の数と \(B\) の行の数が一致している必要があるため、\(AB\) が計算できても、\(BA\) は計算できないという場合もある。\(A\) も \(B\) も正方行列であれば、\(AB\) も \(BA\) は計算はできるがこの値が同じになるとは限らない。行列の積の転置を計算すると、 \((AB)^T\)=\(B^TA^T\) が成り立つ。ベクトルの積について、\({\mathbf x}\) と \({\mathbf y}\) がどちらも \(n\) 次元ベクトルであるとき、

\[\begin{eqnarray} {\mathbf x}^T \cdot {\mathbf y} = {\mathbf y}^T \cdot {\mathbf x} = x_1 y_1 + x_2 y_2 + \cdots +x_n y_n =\sum_{i=1}^{n} x_iy_i \end{eqnarray}\]

で定義される量を ベクトルの 内積という。自分自身との内積

\[\begin{eqnarray} {\mathbf x}^T \cdot {\mathbf x} &=& \sum_{i=1}^{n} x_i^2 =| {\mathbf x} |^2 \end{eqnarray}\]

はベクトルの大きさを表す量であり、これを\(L^2\)ノルムという。 次に、\(n \times n\)型の正方行列について考える。

\[\begin{eqnarray} I =\left( \begin{array}{cccc} 1 &0 & \cdots & 0\\ 0 &1 & \cdots & 0\\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 & \cdots & 1 \\ \end{array} \right) \end{eqnarray}\]

と対角成分だけが \(1\)のとき、同じ\(n \times n\) 型の正方行列 \(A\) に対して、\(IA=AI=A\) が成り立つ。これを単位行列という。ある行列 \(A\) に対して、　\(AB =BA =I\) が成り立つような \(B\)が存在するとき、その\(B\) のことを 逆行列といい、\(A^{-1}\) と表す。

2.4 固有値と固有ベクトル

行列の性質を表す値として 固有値 がある。正方行列\(A\)に対して、ある値(実数か複素数) \(\lambda\) が存在して、 \[ A{\mathbf x} =\lambda {\mathbf x} \]　が成り立つとき \(\lambda\)を 固有値 といい、\(x\) を固有ベクトルという。一般に \(n\times n\) の正方行列 \(A\) に対して固有値を求めるには 固有多項式と呼ばれる \(n\) 次の多項式を解くことになり、固有値は重解を別々に数えると 最大で \(n\) 個ある。今、値の固有値が \(n\) 個あって、 固有ベクトルが \(n\) 個あるとし、行列 \(P\) が

\[ \left\{ \begin{array}{ccc} A{\mathbf x_1} &=& \lambda_1 {\mathbf x_1} \\ A{\mathbf x_2} &=& \lambda_2 {\mathbf x_2} \\ &\vdots& \\ A{\mathbf x_n} &=& \lambda_n {\mathbf x_n} \end{array} \right. \] と書けるものとする。 固有ベクトルは向きだけしか決まらないので、今L2ノルムの大きさが \(1\) のものを選ぶものとする。このベクトルを並べて

\[\begin{eqnarray} P=({\mathbf x_1},{\mathbf x_2}, \cdots,{\mathbf x_n}) \end{eqnarray} \]

という行列を考えると、 \[ AP=P\left( \begin{array}{cccc} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 & \cdots & \lambda_n \end{array} \right) \] となる。この\(P\)が逆行列を持つとすると、

\[\begin{eqnarray} P^{-1}AP &=& \left( \begin{array}{cccc} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 & \cdots & \lambda_n \end{array} \right) \end{eqnarray}\]

と書くことができる。この操作のことを 対角化 という。 対角化された行列に全ての成分が \(1\) のベクトルを掛けると

\[\begin{eqnarray} \nonumber \left( \begin{array}{cccc} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \ddots & \ddots & 0 \\ 0 & 0 & \cdots & \lambda_n \end{array} \right) \left( \begin{array}{c} 1\\ 1\\ \vdots\\ 1 \end{array} \right) &=& \left( \begin{array}{c} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_n \end{array} \right) \end{eqnarray}\]

であり、固有値はそれぞれの成分を \(\lambda_i\)倍 する計算であると 考えられる。したがって、もしあるベクトル \({\mathbf y}\) が

\[\begin{eqnarray} {\mathbf y} &=& a_1 {\mathbf x_1}+ a_2 {\mathbf x_2}+ \cdots +a_n {\mathbf x_n} \end{eqnarray}\]

と表せたとすると、行列 \(A\) を掛ける変換は

\[\begin{eqnarray} A{\mathbf y} = \lambda_1 a_1 {\mathbf x_1}+ \lambda_2 a_2{\mathbf x_2}+ \cdots +\lambda_n a_n{\mathbf x_n} \end{eqnarray}\]

となる。これは各成分を\(\lambda_i\)倍する変換と考えることができる。 このように、対角化することで行列の特徴を把握することができる。 また、もし固有値に一つでも値が0のものがあるとそのベクトルの向きの成分は潰れてしまうということであり、変換によって移った点をもとの点に戻すことはできない。このとき、逆行列を持たないことがわかる。

2.5 正定値行列

\(A\) が対称行列 であるとする。対称行列\(A\) が 固有値 \(\lambda\) を持つとすると、固有値はすべて実数になり、固有ベクトルが 直行することが知られている。また、任意のベクトル\({\mathbf x}\) に対して　

\[\begin{eqnarray} {\mathbf x}^t A {\mathbf x} >0 \end{eqnarray}\]

となるとき、行列 \(A\) を正定値行列という。 である。

\[\begin{eqnarray} {\mathbf x}^t A {\mathbf x} \geq 0 \end{eqnarray}\]

のとき半正定値行列であるという。

\(n\)次の対称行列 \(A\) が、半正定値行列であることは

1. 固有値はすべて非負である。

2. ある直交行列 \(P\) と対角要素がすべて非負の対角行列 D を用いて，\(A = PDP^{T}\) と書くことができる。

3. ある\(n\)次の実行列 \(X\) を用いて、\(A = X^TX\) と書くことができる。

と同値である。分散共分散行列は半正定値行列である。

2.6 Rを用いた行列の計算

行列を作成するには matrix() を指定して、行の数と列の数を指定して行列を作ることができる。

> A <- matrix( c(1,2,2,6),ncol=2,nrow=2)

また、ベクトルや行列などを結合する関数として cbind と rbind という関数がある。cbindはベクトルを列(column)に 並べて結合（右に追加）、rbindは行に並べて結合（下に追加）する。

> x1 <- c(1,2)
> x2 <- c(4,6)
> B <- cbind(x1,x2)
> C <- rbind(x1,x2)
> B

     x1 x2
[1,]  1  4
[2,]  2  6

> C

   [,1] [,2]
x1    1    2
x2    4    6

Rでは、演算子が用意されている。主なものを表 (Table 2.3) に示す。 和や差に対しては、通常の文字の場合と同じように扱うことで成分同士の計算をしてくれる。そこで、行列の積の場合には \(\%*\%\) とする。\(A * B\) としてもエラーではなく、成分同士の掛け算になる。

Table 2.3: Rの行列操作

演算子	説明	例と意味
+	足し算	\(A+B\) 成分ごとの足し算
-	引き算	\(A-B\) 成分ごとの引き算
\(*\)	掛け算	A*B 成分ごとの掛け算
%*%	掛け算	A%*%B
eigen()	固有値	eigen(A) \(A\)は正方行列
solve()	逆行列	solve(A) \(A\)は正方行列
ncol()	列数	列の数を返す
nrow()	行数	行の数を返す
t()	転置	転置行列を返す

固有値を求めてみよう。eigen(A)と打つと次のように結果が表示される。　

> eigen(A)

eigen() decomposition
$values
[1] 6.7015621 0.2984379

$vectors
          [,1]       [,2]
[1,] 0.3310069 -0.9436283
[2,] 0.9436283  0.3310069

この場合２個の異なる固有値、固有ベクトルが求められた。固有値が $values、固有ベクトルが $vectorsである。 固有ベクトルは大きさが \(1\)になるように計算され、縦ベクトルとして表現されている。 また、eigen(A)$valuesとすると、固有値だけを取り出すということを意味する。

固有ベクトルをy1、y2としてみよう。 エクセルのシートのように、1列目だけを表したい場合には、 文字の後に[,1]、 逆に 1行目だけを表すには、[1,]を文字の後に付け加えればよい。

> y1 <- eigen(A)$vectors[,1]
> y2 <- eigen(A)$vectors[,2]
> P <- cbind(y1,y2)
> P

            y1         y2
[1,] 0.3310069 -0.9436283
[2,] 0.9436283  0.3310069

これで、行列の\(P\)までが表すことができた。 後は、\(P^{-1}AP\) を計算して対角成分が固有ベクトルになっていることを確認すればよい。 逆行列を計算するには solve() を用いる。これは連立方程式 \(AX=b\) を解く演算で、solve(A,b)とすると \(X\) を求めてくれる。 今、逆行列を P1 で表そう。

> P1 <- solve(P)
> Q <- P1 %*% A %*% P

結果を見ると

> P1

         [,1]      [,2]
y1  0.3310069 0.9436283
y2 -0.9436283 0.3310069

> Q

              y1           y2
y1  6.701562e+00 3.384065e-16
y2 -3.655004e-17 2.984379e-01

となる。厳密には一致していないが、対角化の値と ほぼ一致していることがわかる。 Rでは通常、数値は浮動小数として扱われ、符号、指数部分を含め、 有限の桁数でしか表現することができない。小数部分も２進数で表現されるので、 10進数で有限桁で表される小数も２進数では循環小数となることもある。 その場合、有限の桁数で打ち切ると誤差を持つこともある。たとえば

> 0.3-0.2-0.1

[1] -2.775558e-17

最後の値は本来 \(0\)であるのに異なった値が出力されている。ここで e-17は10の-17乗を意味し、出てきた答えは \(-0.000\cdots0002775558\)(小数点の前を含めて\(0\)の個数が17個)と 非常に小さい値ではあるが、正しい値にはなっていない。 有限桁で丸める関数として round という関数がある。

> round(0.3-0.2-0.1,16)

[1] 0

このようにコンピュータの計算は誤差を含むことがある。 こうした誤差は通常小さい値であるため実用上問題なく利用しているが、 時としてその違いが本質的となることもあり、 その場合には計算の手順などを工夫して対応することになる。

2.7 関数の引数と変数の範囲

関数を自分で定義して使用することもできる。function(){ } という形で作成する。()の中には変数を指定する。まず例で見てみよう。
第１章で述べたようにRのコマンドは一行で書かなくてもよい。数行に渡る場合には、+という文字が表示されるので、そのまま入力を続ける。 入力が終わるとプロンプト > が表示される。これによって、\(5\times 6\) を計算する関数が定義されたことになる。

> f1 <- function(){ 
+  a <- 5
+  b <- 6
+  a*b 
+ }
> f1

function(){ 
 a <- 5
 b <- 6
 a*b 
}

> f1()

[1] 30

上記の例では、f1 と関数名だけを入力すると、関数の中身が表示される。 実際に計算をさせるためには、f1() とする。すると 30 のように実際に計算がされる。 ここでは、もう少し詳しく関数について見てみよう。 関数で例えば a*b を計算させたいとする。先ほどの例は a、 bも 決まった値だった。今度は、bは固定し、a を毎回変化させることを考えてみよう。

> f2 <- function(a){
+  b <- 6
+  a*b 
+ }
> f2(5) 

[1] 30

ここでは、function(a)とすることで関数がaという 引数を持つことを示している。 また利用するときにはf2(5) のようにする。 f2 の場合、f2()とすると a の値が決まっていないためエラーとなる。

> f2()

Error in f2():  引数 "a" がありませんし、省略時既定値もありません 

もう一つの例を考えてみよう。

>  f3 <- function(a=2){ 
+  b <- 6
+  a*b 
+  }
> f3()

[1] 12

> f3(5)

[1] 30

f3の場合には、何も指定しない場合にはa=2 として計算し12 と出力し、また、f3(5) のように値を指定した場合には a=5 として計算する。 このように省略されたときの値（これを 省略時既定値という）を設定しておくこともできる。

2.8 まとめと展望

行列の積のように自分で手で行うには大変だと思う計算であっても R を用いると瞬時に値を計算することができる。 そこでは、行列を一つのまとまりとして見ることになる。 複雑な計算を行う場合には計算が正しく行えたかを確認する方法を作っておくのも大切だろう。行列を変換という観点からわかりやすく説明した本として [1] がある。また、 例題が多いものとして [2] もある。

参考文献

[1]. 平岡和幸、堀玄、“プログラミングのための線形代数”,オーム社,2004年

[2].Gilbert Strang（著）,松崎公紀、新妻 弘（翻訳） “ストラング線形代数イントロダクション”,近代科学社,2015年

演習

1. 次の R のスクリプトは単位円上の点の集まりを A によって移動させたものである。これを見ると Aという変換が どういう変換をしているのかをイメージしやすい （ベクトルを列ではなくて行にしているため、計算が 転値している）。seq(a,b,c) は aからbまでc刻みで数を発生させる。

N <- 100
theta <-seq(-pi,pi,2*pi/(N+1) )
p_before <- cbind(cos(theta),sin(theta))
plot(p_before)
A <- matrix(c(1,2,3,4),ncol=2,nrow=2)
p_after <- p_before %*% t(A) 
plot(p_after)

また、上の変換で\(A\) の代わりに \(P\)とすると \(P{\mathbf x}\)という計算が回転に対応していることがわかる。

2. 関数について

e1 <-function(a=5){
  b <- 3
  b*a
}

とすると e1() と e1(3) とするとどうなるか？

解答

それぞれ、15 、9　となる。