主成分分析と固有値問題

主成分分析がいかにして固有値問題へと帰着されるのか、詳しめに解説してみる。
少し数式が多め。

目次

主成分分析とは

主成分分析とは、"与えられたデータ情報の次元を圧縮する手法"の一つである。
データの次元を圧縮するとはどういうことだろうか？実は日常生活でもデータを低次元化するということはよく行われている。例えば、肥満度を測るためのBMIという指標が存在する。BMIは体重÷(身長)^2という計算式で導出されるものであるが、よく考えてみると元のデータは(身長,体重)という2成分を持つデータであったのに対し、BMIはただ1成分の数値となっている。これが次元の圧縮であり、情報を削ぎ落としたにも関わらず、肥満度という特徴を表すのに十分な情報を持っている。このように、データに適切な処理を行えば、情報量の削減と特徴の抽出を同時に行う事ができる。また、人がイメージできないような高次元データも、3次元以下に圧縮すればグラフとして描画できるようになる、といったメリットもある。

しかし、現実問題ではBMIの場合のように特徴量があらかじめわかっている場合だけではない。そこで、与えられたデータの傾向から自動的に特徴量を見つけ出し、その特徴を良く表す低次元データへと次元圧縮を行うのが「主成分分析」である(これはある種の機械学習であり、特に自動的に特徴量を見出すという点において、「教師なし学習」と分類されるものである)。

具体的にどのように次元圧縮を行えばいいのだろうか？簡単に、二次元のデータを一次元へと圧縮する場合を考えてみよう。二成分の数値を一成分へと変換する計算式など無数に存在しそうだが、ここでは図のように、データを二次元平面にプロットした時に何らかの直線に向かって全データを射影するという方法を考えることにする。射影した結果、直線の垂直方向の情報が完全に失われ、平行方向の情報のみの一次元データが残る。この残された情報が出来るだけ元データの特徴を良く表していて欲しい。この時、残された一次元データがなるべくバラついている、つまり分散の値が大きいほどよいのである。なぜなら、バラついているということは、各データ点１つ１つの違いをより多く情報として保っていることになるからである。逆に言えば、バラツキの少ない方向というのは、各データが共通して持っている自明な情報なので削除してしまっても問題ないのである。

まとめると、主成分分析では、データの次元圧縮を行う。その際、圧縮後のデータの分散が大きくなるような射影をすることで、特徴量を自動的に抽出するのである。

共分散行列？固有値問題？

統計や数学の勉強を始めると、何に使うのかよくわからない概念が数多く出て来る。
例えば、統計学の教科書や授業などでは「分散」や「共分散」という概念を習う。
「分散」とは、ある一次元のデータ(x1,x2,x3)$(x_1,x_2,x_3)$に対して、平均値x¯¯¯$\overline{x}$からの差の二乗和の平均である。式で書けば次のようになる。

σxx=(x1−x¯¯¯)2+(x2−x¯¯¯)2+(x3−x¯¯¯)23$$\begin{array}{r}{\sigma }_{xx}=\frac{(x_1-\overline{x}{)}^2+(x_2-\overline{x}{)}^2+(x_3-\overline{x}{)}^2}{3}\end{array}$$

「共分散」とは、ある多次元のデータ((x1,y1),(x2,y2),(x3,y3))$((x_1,y_1),(x_2,y_2),(x_3,y_3))$に対して、平均x¯¯¯$\overline{x}$またはy¯¯¯$\overline{y}$からの差、の積平均にあたる。式で書けば次のようになる。

σxy=(x1−x¯¯¯)(y1−y¯¯¯)+(x2−x¯¯¯)(y2−y¯¯¯)+(x3−x¯¯¯)(y3−y¯¯¯)3$$\begin{array}{r}{\sigma }_{xy}=\frac{(x_1-\overline{x})(y_1-\overline{y})+(x_2-\overline{x})(y_2-\overline{y})+(x_3-\overline{x})(y_3-\overline{y})}{3}\end{array}$$

分散は式や定義から、データのバラつきを表していることが直感的にわかる。また、共分散はすぐにはわかりづらいが、xデータとyデータの相関度合いを示すものである。続いて、これらを要素として並べた「共分散行列」なるものが存在する。

Σ==(σxxσxyσxyσyy)13((x1−x¯¯¯)2+(x2−x¯¯¯)2+(x3−x¯¯¯)2(x1−x¯¯¯)(y1−y¯¯¯)+(x2−x¯¯¯)(y2−y¯¯¯)+(x3−x¯¯¯)(y3−y¯¯¯)(x1−x¯¯¯)(y1−y¯¯¯)+(x2−x¯¯¯)(y2−y¯¯¯)+(x3−x¯¯¯)(y3−y¯¯¯)(y1−y¯¯¯)2+(y2−y¯¯¯)2+(y3−y¯¯¯)2)$$\begin{array}{rcl}\mathrm{\Sigma }& =& \left(\begin{array}{cc}{\sigma }_{xx}& {\sigma }_{xy}\\ {\sigma }_{xy}& {\sigma }_{yy}\end{array}\right)\\ & =& \frac{1}{3}\left(\begin{array}{cc}(x_1-\overline{x}{)}^2+(x_2-\overline{x}{)}^2+(x_3-\overline{x}{)}^2& (x_1-\overline{x})(y_1-\overline{y})+(x_2-\overline{x})(y_2-\overline{y})+(x_3-\overline{x})(y_3-\overline{y})\\ (x_1-\overline{x})(y_1-\overline{y})+(x_2-\overline{x})(y_2-\overline{y})+(x_3-\overline{x})(y_3-\overline{y})& (y_1-\overline{y}{)}^2+(y_2-\overline{y}{)}^2+(y_3-\overline{y}{)}^2\end{array}\right)\end{array}$$

といった感じである。行列の形に並べて何の意味があるのか、筆者自身も最近までその使い道を知らなかった。

他には、数学の線型代数において、何らかの行列を与えられた時に固有値と固有ベクトルを求めるという「固有値問題」がよく登場する。
突然なんらかの行列、例えば(1324)$\left(\begin{array}{cc}1& 2\\ 3& 4\end{array}\right)$が与えられた時、次の固有方程式

(1324)(xy)=λ(xy)$$\begin{array}{r}\left(\begin{array}{cc}1& 2\\ 3& 4\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)=\lambda \left(\begin{array}{c}x\\ y\end{array}\right)\end{array}$$

を満たす固有値λ$\lambda$と固有ベクトル(xy)$\left(\begin{array}{c}x\\ y\end{array}\right)$を求めよ、といった具合である。しかし、数学として固有方程式を解く方法だけ学んだとしても、結局何のためにあるのかわからないまま終わってしまう。

これら二つの謎概念(共分散行列、固有値問題)を応用する例の１つが「主成分分析」なのである。

共分散行列についての固有方程式

問題設定

ここからは具体的なデータ群に対して主成分分析を行ってみよう。
少し数学的にはなるが、流れはいたって単純で、

1. ある方向に射影した時のデータの分散を計算する
2. 分散が最小になるような方向を見つける

を行うだけである。結論としては、1を行うと共分散行列が出現し、2を行うとその共分散行列に対する固有方程式が得られる。

シンプルな例として、３つの二次元データ

a⃗ 1=(a1xa1y),a⃗ 2=(a2xa2y),a⃗ 3=(a3xa3y)$$\begin{array}{r}{\overrightarrow{a}}_1=\left(\begin{array}{c}{a}_{1x}\\ a_{1y}\end{array}\right),{\overrightarrow{a}}_2=\left(\begin{array}{c}{a}_{2x}\\ a_{2y}\end{array}\right),{\overrightarrow{a}}_3=\left(\begin{array}{c}{a}_{3x}\\ a_{3y}\end{array}\right)\end{array}$$

についての主成分分析を考えてみよう。簡単のため、重心が原点(0,0)となるようなデータであるとしている。(このように仮定しても問題ないことは最後まで読めばわかるはず)

準備

今後のために、次のようにデータ成分を並べた行列を定義しておく。

X≡=⎛⎝⎜⎜a⃗ T1a⃗ T2a⃗ T3⎞⎠⎟⎟⎛⎝⎜a1xa2xa3xa1ya2ya3y⎞⎠⎟$$\begin{array}{rcl}X& \equiv & \left(\begin{array}{c}{\overrightarrow{a}}_1^{\mathrm{T}}\\ {\overrightarrow{a}}_2^{\mathrm{T}}\\ {\overrightarrow{a}}_3^{\mathrm{T}}\end{array}\right)\\ & =& \left(\begin{array}{cc}{a}_{1x}& a_{1y}\\ a_{2x}& a_{2y}\\ a_{3x}& a_{3y}\end{array}\right)\end{array}$$

転置行列XT$X^{\mathrm{T}}$との積XTX$X^{\mathrm{T}}X$を計算してみると、

XTX==(a1xa1ya2xa2ya3xa3y)⎛⎝⎜a1xa2xa3xa1ya2ya3y⎞⎠⎟(a21x+a22x+a23xa1xa1y+a2xa2y+a3xa3ya1xa1y+a2xa2y+a3xa3ya21y+a22y+a23y)$$\begin{array}{rcl}{X}^{\mathrm{T}}X& =& \left(\begin{array}{ccc}{a}_{1x}& a_{2x}& a_{3x}\\ a_{1y}& a_{2y}& a_{3y}\end{array}\right)\left(\begin{array}{cc}{a}_{1x}& a_{1y}\\ a_{2x}& a_{2y}\\ a_{3x}& a_{3y}\end{array}\right)\\ & =& \left(\begin{array}{cc}{a}_{1x}^2+a_{2x}^2+a_{3x}^2& a_{1x}{a}_{1y}+a_{2x}{a}_{2y}+a_{3x}{a}_{3y}\\ a_{1x}{a}_{1y}+a_{2x}{a}_{2y}+a_{3x}{a}_{3y}& a_{1y}^2+a_{2y}^2+a_{3y}^2\end{array}\right)\end{array}$$

最終的な表式が、上で見た共分散行列の形になっていることがわかる("1/データ数"の因子を除いて)。共分散行列は後で出てくるので、次のようにΣ$\mathrm{\Sigma }$として定義しておこう。

Σ≡13XTX$$\begin{array}{r}\mathrm{\Sigma }\equiv \frac{1}{3}{X}^{\mathrm{T}}X\end{array}$$

1. 単位ベクトルにデータを射影した時の分散を計算

ある方向の単位ベクトルをe⃗ =(exey)$\overrightarrow{e}=\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)$と表す。データをe⃗ $\overrightarrow{e}$の方向に射影して、1次元データとなった後の分散値の表式を求める。例えば、a⃗ 1${\overrightarrow{a}}_1$をe⃗ $\overrightarrow{e}$に射影した値はa⃗ 1⋅e⃗ ${\overrightarrow{a}}_1\cdot \overrightarrow{e}$のようにベクトルの内積を取ることで簡単に計算でき、これを用いて分散(var)を定義通りに計算すると、

var=(a⃗ 1⋅e⃗ )2+(a⃗ 2⋅e⃗ )2+(a⃗ 3⋅e⃗ )23$$\begin{array}{r}\mathrm{v}\mathrm{a}\mathrm{r}=\frac{({\overrightarrow{a}}_1\cdot \overrightarrow{e}{)}^2+({\overrightarrow{a}}_2\cdot \overrightarrow{e}{)}^2+({\overrightarrow{a}}_3\cdot \overrightarrow{e}{)}^2}{3}\end{array}$$

となる(平均値は0となっていることに注意)。これを行列の積の形に分解して表現してみると、次に示す式の一行目のようになる。実は、この分解した行列がそれぞれ(Xe⃗ $X\overrightarrow{e}$)とその転置行列(Xe⃗ )T$(X\overrightarrow{e}{)}^{\mathrm{T}}$になっている(実際に計算してみれば確認できる)。二行目から三行目では行列の性質を用いている。こうするとXTX$X^{\mathrm{T}}X$という並びが出現する。これは上で計算した共分散行列Σ$\mathrm{\Sigma }$である。

var====13(a⃗ 1⋅e⃗ a⃗ 2⋅e⃗ a⃗ 3⋅e⃗ )⎛⎝⎜a⃗ 1⋅e⃗ a⃗ 2⋅e⃗ a⃗ 3⋅e⃗ ⎞⎠⎟13(Xe⃗ )T(Xe⃗ )13e⃗ TXTXe⃗ e⃗ TΣe⃗ $$\begin{array}{rcl}\mathrm{v}\mathrm{a}\mathrm{r}& =& \frac{1}{3}\left(\begin{array}{ccc}{\overrightarrow{a}}_1\cdot \overrightarrow{e}& {\overrightarrow{a}}_2\cdot \overrightarrow{e}& {\overrightarrow{a}}_3\cdot \overrightarrow{e}\end{array}\right)\left(\begin{array}{c}{\overrightarrow{a}}_1\cdot \overrightarrow{e}\\ {\overrightarrow{a}}_2\cdot \overrightarrow{e}\\ {\overrightarrow{a}}_3\cdot \overrightarrow{e}\end{array}\right)\\ & =& \frac{1}{3}(X\overrightarrow{e}{)}^{\mathrm{T}}(X\overrightarrow{e})\\ & =& \frac{1}{3}{\overrightarrow{e}}^{\mathrm{T}}{X}^{\mathrm{T}}X\overrightarrow{e}\\ & =& {\overrightarrow{e}}^{\mathrm{T}}\mathrm{\Sigma }\overrightarrow{e}\end{array}$$

このvarの表式はデータが多成分(つまり高次元)、多データ化しても一般的に成り立つはずである。なぜなら、行列Xやベクトルe⃗ $\overrightarrow{e}$の縦成分と横成分数をデータ次元数とデータ数に併せて拡張すれば同じような計算が行えるからである。

2. varが最小になる方向を求める

さて、主成分分析ではこのvarを最小化するようなe⃗ $\overrightarrow{e}$を求めることが目的であった。また、e⃗ $\overrightarrow{e}$には単位ベクトル、つまりノルムが1であるという制約条件も付いている。このように、ある制約条件のもとで、関数の最大(最小)を決定する時には、「ラグランジュの未定乗数法」というものがよく用いられる。

　ラグランジュの未定乗数法:
　　変数x,yについて、g(x,y)=0$g(x,y)=0$の制約条件下でf(x,y)$f(x,y)$という関数を最小化するのは、

L(x,y,λ)≡f(x,y)−λg(x,y)$L(x,y,\lambda )\equiv f(x,y)-\lambda g(x,y)$

　　という関数を定義した時に、

∂L∂x=0,∂L∂y=0,∂L∂λ=0$\frac{\mathrm{\partial }L}{\mathrm{\partial }x}=0,\frac{\mathrm{\partial }L}{\mathrm{\partial }y}=0,\frac{\mathrm{\partial }L}{\mathrm{\partial }\lambda }=0$

　　を満たすような解(x,y)$(x,y)$である。

ラグランジュの未定乗数法を用いると、

L(ex,ey,λ)≡=var−λ(e2x+e2y−1)(exey)Σ(exey)−λ(e2x+e2y−1)$$\begin{array}{rcl}L(e_x,e_y,\lambda )& \equiv & \mathrm{v}\mathrm{a}\mathrm{r}-\lambda (e_x^2+e_y^2-1)\\ & =& \left(\begin{array}{cc}{e}_x& e_y\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)-\lambda (e_x^2+e_y^2-1)\end{array}$$

と定義した関数L(ex,ey,λ)$L(e_x,e_y,\lambda )$に対して、

∂L(ex,ey,λ)∂ex=0∂L(ex,ey,λ)∂ey=0$$\begin{array}{r}\frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }{e}_x}=0\\ \frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }{e}_y}=0\end{array}$$

が成立するようなex,ey$e_x,e_y$を見つければよい(λ$\lambda$についての偏微分は一旦無視する)。
∂L(ex,ey,λ)∂ex$\frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }{e}_x}$,∂L(ex,ey,λ)∂ey$\frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }{e}_y}$を計算すると、

∂L(ex,ey,λ)∂ex∂L(ex,ey,λ)∂ey=====(10)Σ(exey)+(exey)Σ(10)−2λex(10)Σ(exey)+[(exey)Σ(10)]T−2λex(10)Σ(exey)+(10)ΣT(exey)−2λex2(10)Σ(exey)−2λex2(01)Σ(exey)−2λey$$\begin{array}{rcl}\frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }{e}_x}& =& \left(\begin{array}{cc}1& 0\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)+\left(\begin{array}{cc}{e}_x& e_y\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}1\\ 0\end{array}\right)-2\lambda e_x\\ & =& \left(\begin{array}{cc}1& 0\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)+{\left[\left(\begin{array}{cc}{e}_x& e_y\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}1\\ 0\end{array}\right)\right]}^{\mathrm{T}}-2\lambda e_x\\ & =& \left(\begin{array}{cc}1& 0\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)+\left(\begin{array}{cc}1& 0\end{array}\right){\mathrm{\Sigma }}^{\mathrm{T}}\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)-2\lambda e_x\\ & =& 2\left(\begin{array}{cc}1& 0\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)-2\lambda e_x\\ \frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }{e}_y}& =& 2\left(\begin{array}{cc}0& 1\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)-2\lambda e_y\end{array}$$

であり、これらが0になるので

(10)Σ(exey)(01)Σ(exey)==λexλey$$\begin{array}{rcl}\left(\begin{array}{cc}1& 0\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)& =& \lambda e_x\\ \left(\begin{array}{cc}0& 1\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)& =& \lambda e_y\end{array}$$

と式変形出来る。2つの方程式をまとめると、

(1001)Σ(exey)=λ(exey)$$\begin{array}{rcl}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)& =& \lambda \left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)\end{array}$$

一番左の行列はただの単位行列なので省略して、

Σ(exey)=λ(exey)$$\begin{array}{rcl}\mathrm{\Sigma }\left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)& =& \lambda \left(\begin{array}{c}{e}_x\\ e_y\end{array}\right)\end{array}$$

これは共分散行列Σ$\mathrm{\Sigma }$に対する固有方程式となっている！！
(上で無視していた∂L(ex,ey,λ)∂λ=0$\frac{\mathrm{\partial }L(e_x,e_y,\lambda )}{\mathrm{\partial }\lambda }=0$、つまりe2x+e2y=1$e_x^2+e_y^2=1$というベクトルの正規化条件は、固有値問題を解いた後の固有ベクトルを拡大縮小することで調節出来る。)
実際にデータを主成分分析する際には、データから共分散行列を生成し、固有ベクトルを計算、全データを射影する、といった作業を行えば良い。プログラミング言語によってはこれらの作業を簡略化出来るパッケージが存在したりするので、とても簡単に処理出来る。

参考リンク

主成分分析に関しては、以下の解説およびブログを参照した。
http://manabukano.brilliant-future.net/document/text-PCA.pdf
yusuke-ujitoko.hatenablog.com