こんにちは！@tsum_tsumです。
今回は、主成分分析（PCA, Principal Component Analysis）について、数式をいっぱい交えて紹介していきたいと思います！
筆者は統計検定準1級までは合格しており、PCAはその範囲の1つになっています。公式本で勉強はするのですが、説明も短いし行間も広いので、よく分からないまま式だけ頑張って覚えたという状態になっています。数学においてこの状態は良くないと思い、改めて勉強しようと思った次第です。
そんな状態なので、誤字脱字、認識違い等あるかもしれません（もちろんできる限り無いよう、善処しております）。もし発見されましたら、お手数をかけてしまいますが、コメント欄でご指摘いただければと思います。
また、こちらの文章は統計検定準1級を受験しようとしている方にも役立つように書いているつもりでございます。勉強の際のご参考にしていただければ幸いです。
では前書きもこの辺りで、本編に進んでまいりましょう。
はじめに
まず、主成分分析とは何でしょうか？Wikipediaによると、下記のように記載されております。
相関のある多数の変数から相関のない少数で全体のばらつきを最もよく表す主成分と呼ばれる変数を合成する多変量解析の一手法[1]。
少し言葉が堅苦しいので、もう少しかみ砕いて説明していきましょう。
あらゆるデータは座標空間上の点として表現することができます。身長 $170cm$、体重 $65kg$ は 2次元の座標空間（座標平面）上の $(170, 65)$ という点で表現できる、といった感じです。
しかし、1データの持つ特徴量が10次元、100次元となれば、人間がデータの全体像を認識することはできません。そこで使われるのが、次元削減です。1データの持つ特徴量が10次元から2次元になれば、人間の目で確認することもできますし、計算量の削減にもつながります。
一方で、適当に次元削減してもいいわけではありません。もしあなたが身長のデータしか必要ないのであれば、体重のデータを無視することで1データの持つ特徴量を「身長」「体重」から「身長」に減らすことができますね。しかし、すべてのデータの分布を確認したいこともあるはずです。
そこで主成分分析の出番です。主成分分析はデータの分布全体を考慮しながら、情報をできるだけ失わずに次元を削減する手法です。 具体的には、下記のようにして計算されます。
- データのばらつきが最大になる方向に新しい軸を設定する
- その軸にデータを投影する
これでも分かりにくい方はこちらの画像をご覧ください。
Data points はランダムに生成された点で2次元の座標平面上に散らばって分布しています。このデータが最もばらつきが多い軸を探すのが主成分分析です。ちなみにこの画像では赤矢印（principal Component）が最もばらつきが大きい軸です。そして、最もばらつきが大きい軸にデータを写した画像が下記になります。
今回の記事では、
- なぜこのような軸（画像の赤矢印）を見つけることができるのか？
- どのような計算が行われているのか？
に注目して 数式をたくさん交えて 説明していきたいと思います。
主成分分析の計算
それでは、実際に計算していきましょう。とはいえ、いきなり任意の次元にすると、数式についていくのが大変になりそうですね。ということで簡単に2次元から1次元の主成分分析を考えていきましょう。ここでの計算の目的は上画像の赤矢印の導出です。
$N$ 個の2次元データを下記のように表します。
{\rm x}_i = (x_{i1}, x_{i2}), \ (i=1, \cdots, N).
これを1次元に圧縮することを考えます。
y_i = {\rm x}_i {\rm w} =
(x_{i1}, x_{i2}) \begin{pmatrix} w_1 \\ w_2 \end{pmatrix}.
ただし、${\rm w} = \begin{pmatrix} w_1 \ w_2 \end{pmatrix}^T$ としておきます。さらに、下記の記号も導入しておきます。
\begin{eqnarray}
{\rm y} = (y_1, \cdots, y_N)^T &=& {\rm X}{\rm w} \\
&=& \begin{pmatrix}
x_{11} & x_{12} \\
x_{21} & x_{22} \\
\vdots & \vdots \\
x_{N1} & x_{N2}
\end{pmatrix}
\begin{pmatrix}
w_1 \\
w_2 \\
\end{pmatrix}.
\end{eqnarray}
今、${\rm x}$ が2次元の元データに対し、$y$ が1次元に圧縮されたデータですね。では、$y$ の不偏分散 $s_y^2$ を求めていきましょう。すると、このように計算ができるかと思います。
\begin{eqnarray}
s_y^2 &=& \frac{1}{N-1}\sum_{i=1}^N (y_i - {\bar y})^2 \\
&=& \frac{1}{N-1}\sum_{i=1}^N \{
(x_{i1} - {\bar x}_{*1})w_1 + (x_{i2} - {\bar x}_{*2})w_2
\}^2 \\
&=& \frac{1}{N-1} \left\{
w_1^2 \sum_{i=1}^N (x_{i1} - {\bar x}_{*1})^2
+ 2w_1w_2 \sum_{i=1}^N (x_{i1} - {\bar x}_{*1})(x_{i2} - {\bar x}_{*2})
+ w_2 \sum_{i=1}^N (x_{i2} - {\bar x}_{*2})^2
\right\} \\
&=& (w_1 \ w_2)
\begin{pmatrix} s_{11} & s_{12} \\ s_{21} & s_{22} \end{pmatrix}
\begin{pmatrix} w_1 \\ w_2 \end{pmatrix} \\
&=& {\rm w}^T {\rm S} {\rm w}
\end{eqnarray}
ただし、
{\bar y}=\frac{1}{N}\sum_{i=1}^N y_i
は $y$ の標本平均で、
\begin{eqnarray}
{\bar x}_{*1} &=& \frac{1}{N}\sum_{i=1}^N x_{i1} , \\
{\bar x}_{*2} &=& \frac{1}{N}\sum_{i=1}^N x_{i2}
\end{eqnarray}
はそれぞれ、${\rm x}$ の第一成分、第二成分の標本平均を表します。また、$s_{jk} \ (j,k \in \{1,2\})$ は（不偏）分散または共分散
s_{jk} = \frac{1}{N-1}\sum_{i=1}^N (x_{ij} - {\bar x}_{*j})(x_{ik} - {\bar x}_{*k})
を表しています。よって、${\rm S}$ は分散共分散行列を表しています。
さて、ここまで $s_y^2$ の計算をしていましたが、その目的は何だったでしょうか？それは、最もばらつきが大きくなる軸（上画像の赤矢印）を求めることでしたね。ということは、$s_y^2$ の最大値を求めればいいということが分かります。そして、$s_y^2$ が最大値になるときの ${\rm w}$ が上画像の赤矢印になります！
最大値の求め方
今求めたいのは $s_y^2$ が最大値を取るような、${\rm w}$ の値でしたね。数式で表現すると以下のようになります。
s_y^2 = \underset{{\rm w}}{\max} ({\rm w}^T {\rm S} {\rm w}).
ただし、このままではいくらでも大きな $w$ を取れば $s_y^2$ を大きくすることができます。今回は軸の向きを知りたいだけなので、${\rm w}$ の大きさに制約条件
|{\rm w}|^2 = {\rm w}{\rm w}^T = w_1^2 + w_2^2 = 1
を設けておきます。すると、これはラグランジュの未定乗数法で解くことができますね。ということで実際に解いていきます。
L(w_1, w_2, \lambda) = {\rm w}^T {\rm S} {\rm w} - \lambda ({\rm w}{\rm w}^T -1)
として、下記の連立方程式を解きます。
\begin{cases}
\frac{\partial L}{\partial {\rm w}} = 2 {\rm S}{\rm w} - 2 \lambda {\rm w} = 0 \\
\frac{\partial L}{\partial \lambda} = {\rm w}{\rm w}^T -1 = 0 \\
\end{cases}
から、${\rm S}{\rm w} = \lambda {\rm w}$ が得られます。
ちなみにラグランジュの未定乗数法では、「ある点 ${\rm w}={\rm w}_0$ で極値を取るなら
\begin{cases}
\frac{\partial L}{\partial {\rm w}}({\rm w}_0) = 0 \\
\frac{\partial L}{\partial \lambda}({\rm w}_0) = 0 \\
\end{cases}
が成り立つ」という命題は正しいですが、その逆が成り立つとは限りません。よって、${\rm w}={\rm w}_0$ が最大値を取るかどうかはチェックが必要になります。今回は最大値を取るとして進めてまいりましょう。
さて、得られた式 ${\rm S}{\rm w} = \lambda {\rm w}$ を確認しましょう。この式をどこかでご存じないでしょうか？この式は固有値・固有ベクトルを求めるときの式（固有方程式） になっていますね！実は、ばらつきが最も大きくなるとき（$s_y^2$ が最大になるとき）の軸（赤矢印）とは、分散共分散行列 ${\rm S}$ の固有ベクトルのことだったんですね！
さらに、
\begin{eqnarray}
{\rm S}{\rm w} &=& \lambda {\rm w} \\
{\rm w}^T{\rm S}{\rm w} &=& \lambda {\rm w}^T{\rm w} \\
s_y^2 &=& \lambda
\end{eqnarray}
から、$s_y^2$ の最大値は分散共分散行列 ${\rm S}$ の固有値 になっています！
ここまでのまとめ
ここまでの内容を少し整理しましょう。今求めようとしていたのは、画像の赤矢印、つまり
- 最もばらつき（分散）が大きくなる軸の方向を探すこと
でしたね。そして、
- その軸が「固有ベクトル」で、その時の最大値が「固有値」である
という事実を確認していきました。しかし、まだ疑問が残っていますね。固有値・固有ベクトルは2つあるはずです。もう1つは何を表しているのでしょうか？
もう1つの固有値・固有ベクトル
ここからはもう1つの固有ベクトル・固有値の計算を行っていきましょう。まず、分散共分散行列 ${\rm S}$ について振り返っていきましょう。
{\rm S} = \begin{pmatrix} s_{11} & s_{12} \\ s_{21} & s_{22} \end{pmatrix}
であって、
s_{jk} = \frac{1}{N-1}\sum_{i=1}^N (x_{ij} - {\bar x}_{*j})(x_{ik} - {\bar x}_{*k})
でした。$s_{jk}$ は $s_{kj}$ としても、同じ値になりますね。よって、${\rm S}$ は対称行列になります。また、明らかにこれらは実数ですね。そして、実対称行列では次のような性質を持っています。
定理
$A$ を実対称行列とする。このとき、次のことが成り立つ。
- 行列の固有値は全て実数である
- 行列の相異なる固有ベクトル同士は全て直交する
- $A$ は直交行列で対角化可能である
言い換えると、$A$ の固有値を $\lambda_1, \cdots, \lambda_p$ 、対応する固有ベクトルを ${\rm u}_1, \cdots, {\rm u}_p$ とすると、下記が成り立ちます。
- $\lambda_1, \cdots, \lambda_p \in {\mathbb R}, $
- $i \ne j$ のとき、${\rm u}_i \perp {\rm u}_j, $
- $ {\rm U}^T {\rm A} {\rm U} = {\operatorname {diag}(\lambda_1, \cdots, \lambda_p)} $ 、ただし ${\rm U}=({\rm u}_1, \cdots, {\rm u}_p).$
分散共分散行列の固有ベクトルをあらためて確認してみる
先ほどの例に戻って確認してみましょう。先ほどは、2次元のデータについて考えていました。ですので、先の1~3番をあらためて書き起こしてみます。
- $\lambda_1, \lambda_2 \in {\mathbb R}, $
- ${\rm u}_1 \perp {\rm u}_2, $
- $ {\rm U}^T {\rm S} {\rm U} = {\operatorname {diag}(\lambda_1, \lambda_2)} $ 、ただし ${\rm U}=({\rm u}_1, {\rm u}_2).$
さて、1番が成り立つので $\lambda_1 \ge \lambda_2$ とします。すると、$\lambda_1$ が一番大きいので対応する ${\rm u}_1$ が赤矢印ということになります（上画像の赤矢印は固有ベクトルだったことを思い出してください）。つまり、この ${\rm u}_1$ こそが赤矢印の正体になっています。
また、2番の式を見てみると、${\rm u}_1 \perp {\rm u}_2$ となっています。この式は ${\rm u}_1$ と ${\rm u}_2$ が直交するということを表しています。つまり、もう1つの固有ベクトルは赤矢印に直交する矢印（下画像の緑矢印） になるということが分かりますね！
また、その固有値は緑矢印にデータの点を投影したときの分散を表しています！
ここまでのまとめ
さて、この章では、もう1つの固有値・固有ベクトルの正体について説明してまいりました。
- 固有値は固有ベクトルの向きにデータを投影した時の分散
- 最も大きい固有値 $\lambda_1$ に対応する固有ベクトル ${\rm u}_1$ は分散が最も大きくなる軸の向き
- 2番目に大きい固有値 $\lambda_2$ に対応する固有ベクトル ${\rm u}_2$ は ${\rm u}_1$ に直交する軸のうち、分散が最も大きくなる軸の向き
ということが分かりましたね。
ここまででも主成分分析について概ねご理解いただけたかと思います。ここからは、統計検定準1級の範囲を念頭に置いて、主成分分析の寄与率というものについて説明していきます。主成分分析そのものの話からは少し逸れてきますので、主成分分析の意味だけを知りたい方はまとめの方にお進みいただければ幸いです。
寄与率とは？
主成分分析には寄与率というものがよく現れます。これは何を表しているのでしょうか？それを説明するために、分散共分散行列 ${\rm S}$ について少し深堀していきます。
${\rm S}$ の固有値を $\lambda_1, \lambda_2(\lambda_1 \ge \lambda_2)$ とすると、これらは散らばり具合（分散）を表すということを説明しました。しかし、ここで 「いつでも固有値が正の数になるとは限らないのでは？」 という疑問が浮かんできた方はかなり鋭いと思います。もし負の数になることがあれば、分散や共分散の値として相応しくないのではないのか…。
しかし、その点は実は問題ありません。分散共分散行列のような半正定値行列の場合、固有値は全て正の数または $0$ になることが知られているからです。
正定値行列とは
$N \times N$ の正方行列 $A$ が正定値行列であるとは、任意の ${\rm t} = (t_1, \cdots, t_N)^T \ne {\rm 0}$ に対して、
{\rm t}^T A {\rm t} \gt 0
が成り立つことをいいます。ちなみに $\ge 0$ におきかえることで半正定値行列といわれます。
分散共分散行列 ${\rm S}$ に戻ってみると、任意の ${\rm t} = (t_1, t_2)^T \ne (0,0)^T$ に対して、
{\rm t}^T {\rm S} {\rm t} \ge 0
が成り立てば、${\rm S}$ が半正定値行列であることが分かりますが、これは $s_y^2$ の計算を思い出して、
\begin{eqnarray}
{\rm t}^T {\rm S} {\rm t}
&=& \frac{1}{N-1}\sum_{i=1}^N \{
(x_{i1} - {\bar x}_{*1})t_1 + (x_{i2} - {\bar x}_{*2})t_2
\}^2
\ge 0
\end{eqnarray}
からお分かりだと思います。2乗の和を取っているだけですね。そして、正定値行列（半正定値行列）の固有値はいずれも正（正または $0$ ）になります。なぜなら、
\begin{eqnarray}
{\rm t}^T {\rm S} {\rm t}
&=& {\rm t}^T \lambda {\rm t} \\
&=& \lambda {\rm t}^T {\rm t} \\
&=& \lambda |{\rm t}|^2
\end{eqnarray}
と式変形し、${\rm t}^T {\rm S} {\rm t} \gt 0, \ |{\rm t}|^2 \ge 0$ から $\lambda \gt 0$ が導かれるからですね。
寄与率について
※以降はデータの点集合は平均が $0$ になるように正規化されていると仮定します。 つまり、
{\rm x}_i = (x_{i1}, x_{i2}), \ (i=1, \cdots, N).
に対し、
{\bar x}_{*1} = \frac{1}{N} \sum_{i=1}^N x_{i1}, \ {\bar x}_{*2} = \frac{1}{N} \sum_{i=1}^N x_{i2}
がいずれも $0$ になるように正規化されています。（注意終わり）
ということで、いよいよ寄与率について説明していきましょう。先ほど、
{\rm U}^T {\rm S} {\rm U} = {\operatorname {diag}(\lambda_1, \lambda_2)}
という式を見たかと思います。さて、両辺のトレースを取ってみます。すると、
{\operatorname {tr}} ({\rm U}^T {\rm S} {\rm U}) = {\operatorname {tr}}({\rm S}) = s_{11}+s_{22} = \lambda_1 + \lambda_2
が導かれます。この式は各固有値 $\lambda_i$ が各元データの散らばり具合（ここでは $s_{11}, s_{22}$ ）を反映しているかを表しているんですね！ということで、散らばり具合を反映している度合い、つまり、
\frac{\lambda_i}{\lambda_1 + \lambda_2}
を第$i$主成分の寄与率といいます！
主成分負荷量とは？
最後に主成分負荷量を紹介しましょう。主成分負荷量とは、主成分と元の変数の相関係数のことです。順を追って説明していきましょう。
{\rm S} = \begin{pmatrix} s_{11} & s_{12} \\ s_{21} & s_{22} \end{pmatrix}
を分散共分散行列、$\lambda_1 \ge \lambda_2$ を ${\rm S}$ の固有値、
{\rm u}_1=(u_{11}, u_{12})^T, {\rm u}_2=(u_{21}, u_{22})^T
を対応する固有ベクトルとします。このとき、
y_i = {\rm x}^T{{\rm u}_i} = (x_1, x_2) \begin{pmatrix} u_{i1} \\ u_{i2} \end{pmatrix} = x_1u_{i1} + x_2u_{i2}
を第$i$主成分といいます。これは、データの点を ${\rm u}_i$ の軸に投影した後の点の位置と解釈できるかと思います。
さて、主成分負荷量とは、主成分と元の変数の相関係数のことと述べました。主成分と元の変数の相関係数とは何を表しているのでしょうか？これは、元の変数が主成分にどれくらい影響を表しているのかを表しています。
例えば、身長、体重の2次元データを考えてみましょう。データを主成分分析すると、最も散らばりの大きい軸の向き（画像の赤矢印）を求めることができました。しかし、その矢印の向きがどれくらい体重（または身長）に影響を受けたのか分かりませんよね。
そこで、主成分負荷量です。主成分負荷量を求めることで、矢印の向きに対する体重（または身長）の影響度が分かるということなんですね。
ということで計算を続けていきましょう。実際に計算してみます。$x_i, y_j$ 間の主成分負荷量は、下記のように求められます。
\begin{eqnarray}
r_{x_i, y_j}
&=& \frac{\operatorname{Cov}[x_i, y_j]}{\sqrt{V[x_i]V[y_j]}} \\
&=& \frac{\operatorname{Cov}[x_i, y_j]}{\sqrt{s_{ii}\lambda_j}}.
\end{eqnarray}
ここで、$V[x_i] = s_{ii}, V[y_j]=\lambda_j$ を思い出してください。また、分子はこのように計算できます。
\begin{eqnarray}
\operatorname{Cov}[x_i, y_j]
&=& E[x_iy_j] - E[x_i]E[y_j] \\
&=& E[x_i(x_1u_{j1}+x_2u_{j2})] - E[x_i]E[y_j] \\
&=& u_{j1} E[x_ix_1]+u_{j2}E[x_ix_2] \\
&=& u_{j1} s_{i1}+u_{j2}s_{i2}. \\
\end{eqnarray}
計算は少し複雑かもしれませんが、
\operatorname{Cov}[x_i, x_1] = E[x_ix_1] - E[x_i]E[x_1]
を思い出せば、それほど難しくないかと思います。よって、下記のように計算することができます。
\begin{eqnarray}
\begin{pmatrix}
\operatorname{Cov}[x_1, y_j] \\
\operatorname{Cov}[x_2, y_j]
\end{pmatrix}
&=& \begin{pmatrix}
u_{j1} s_{11}+u_{j2}s_{12} \\
u_{j1} s_{21}+u_{j2}s_{22}
\end{pmatrix} \\
&=& \begin{pmatrix}
s_{11} & s_{12} \\
s_{21} & s_{22}
\end{pmatrix}
\begin{pmatrix}
u_{j1} \\
u_{j2}
\end{pmatrix} \\
&=& {\rm S}{\rm u}_j \\
&=& \lambda_j {\rm u}_j \\
&=& \begin{pmatrix}
\lambda_j u_{j1} \\
\lambda_j u_{j2}
\end{pmatrix} \\
\end{eqnarray}
ということで、$\operatorname{Cov}[x_i, y_j] = \lambda_j u_{ji}$ が分かりました。この式を最後に代入して、
\begin{eqnarray}
r_{x_i, y_j}
&=& \frac{\operatorname{Cov}[x_i, y_j]}{s_{ii}\sqrt{\lambda_j}} \\
&=& \frac{\lambda_j u_{ji}}{\sqrt{s_{ii}\lambda_j}} \\
&=& \frac{\sqrt{\lambda_j} u_{ji}}{\sqrt{s_{ii}}} \\
\end{eqnarray}
が得られました。これが主成分負荷量です。お疲れさまでした。
まとめ
ということで今回は主成分分析をできるだけ詳しくまとめてみました。主成分分析はデータの点の集まりから最も散らばっている軸の向きを探す計算手法ということがお分かりいただけたかと思います。そして、それを求めるには、分散共分散行列の固有値・固有ベクトルを求めればいいということもお分かりいただけたかと思います。
その後は、寄与率と主成分負荷量について説明していきました。この資料が皆さんの学習への参考になれば幸いです。
また、今回は簡単のため2次元を例にして計算を紹介しました。計算の練習として一般の次元で計算してみても面白いかもしれません。
それでは皆さんまたお会いしましょう！