PCA(Principal Component Analysis)を具体的に計算する


このエントリーをはてなブックマークに追加

(2013/6/18修正)

PCAの計算方法を天下り的に説明します。専門家ではないので誤りが含まれる可能性があることをあらかじめご了承ください。データの例はExample of Principal Component Analysis (PCA).mp4 - YouTubeで使われているものをそのまま流用しています。というか、このあやしい説明を読むよりかこの動画を見た方がいいかもしれません。

PCAの目的

ベクトル基底を変換することにより、できるだけ情報を失うことなく次元を削減すること。

データの例

2つの独立変数X, Yを持つ10組のデータがあるとする。

X Y
1.4 1.65
1.6 1.975
-1.4 -1.775
-2.0 -2.525
-3.0 -3.95
2.4 3.075
1.5 2.025
2.3 2.75
-3.2 -4.05
-4.1 -4.85

プロットするとこんな感じ。いかにもXとYに相関がありそうなデータである。

ステップ1: 共分散行列(Covariance Matrix)を求める

共分散行列の前に、共分散(Covariance)とは何なのか。Wikipediaによると、共分散とは、「2組の対応するデータ間での、平均からの偏差の積の平均値」である。多分、共分散の定義は2種類ある。1つ目の定義は標本共分散(??)で、
Cov(X, Y) = \frac{1}{N}\Sigma_i (x_i-E[X])(y_i-E[Y])
2つ目の定義は不偏共分散(??)で、
Cov(X, Y) = \frac{1}{N-1}\Sigma_i (x_i-E[X])(y_i-E[Y])
上で使った「標本共分散」と「不偏共分散」、ほんとにこのような用語の使い方で合ってるか全然自信がない。あてはまる英語の訳もよくわからない。この稿で参考にしているYouTubeでは、不偏共分散の方を使って計算を進めているようなので、ここでもそれに従うことにする。

変数X, Yに対して共分散行列は以下のように定義される。
\begin{bmatrix} Cov(X, X) & Cov(X, Y) \\ Cov(Y, X) & Cov(Y, Y) \end{bmatrix}
なお、Wikipediaでは共分散行列ではなく分散共分散行列 - Wikipediaという語で項目が立っているが、これらは同じ物をさしているようだ。分散共分散行列は、行列の対角成分が分散になっていることから来た名前のようだ。

実際に共分散行列を求めてみる。まずE[X], E[Y]は単にそれぞれX, Yの平均値なので
E[X] = (1.4 + 1.6 + ... + (-4.1)) / 10 = -0.45
E[Y] = (1.65 + 1.975 + ... + (-4.85)) / 10 = -0.5675


次にCov(X, Y), Cov(Y, X)を求めると
Cov(X, Y) = Cov(Y, X) = \{ (1.4 - (-0.45)) * (1.65 - (-0.5675)) + ... + (-4.1 - (-0.45)) * (-3.65 - (-0.5675)) \} / (10 - 1)\\= 7.987638889

同様にCov(X, X), Cov(Y, Y)を求めると
Cov(X, X) = \{(1.4 - (-0.45))^2 + ... + (-4.1 - (-0.45))^2 \} / (10-1) = 6.422777778\\Cov(Y, Y) = \{(1.65 - (-0.5675))^2 + ... + (-3.65 - (-0.5675))^2 \} / (10-1) = 9.952784722
ちなみに、これらの値は定義よりそれぞれXの不偏分散、Yの不偏分散と等しい。

以上から、共分散行列は以下のような対称行列となる。
\begin{bmatrix} 6.422777778 & 7.987638889 \\ 7.987638889 & 9.952784722 \end{bmatrix}

ステップ2, 3: 共分散行列の固有値固有ベクトルを求める

この共分散行列の固有値を求める。固有値の求め方は省略。今回はさぼってWolfram Alpha先生に聞いてみると、
二つの固有値\lambda_1 =16.3681, \lambda_2 = 0.00746266
それぞれに対応する固有ベクトル(長さ1とする)はv_1 = \begin{bmatrix} 0.626194 \\ 0.779667 \end{bmatrix}, v_2 = \begin{bmatrix} -0.779667 \\ 0.626194 \end{bmatrix}
と出た。v_1, v_2は互いに直行しているのに注意。この固有ベクトルv_1, v_2が、PCAによって求められた新しい軸に相当する。

ステップ4: 座標を変換する

さっき求めた固有ベクトル固有値の大きい順に並べると、座標の変換行列
\begin{bmatrix} 0.626194 & -0.779667 \\ 0.779667 & 0.626194 \end{bmatrix}
ができる。変換行列を介した旧座標の点(x,y)と新座標の点(x',y')の関係は以下のようになっている。

\begin{bmatrix} x' & y' \end{bmatrix} = \begin{bmatrix} x - E(X) & y - E(Y) \end{bmatrix} \begin{bmatrix} 0.626194 & -0.779667 \\ 0.779667 & 0.626194 \end{bmatrix}

例えば、点(x, y) = (1.4, 1.65)が新座標ではどうなるかを計算すると、
\begin{bmatrix} x' & y' \end{bmatrix} = \begin{bmatrix} 1.4 - (-0.45) & 1.65 - (-0.5675) \end{bmatrix} \begin{bmatrix} 0.626194 & -0.779667 \\ 0.779667 & 0.626194 \end{bmatrix} \\ = \begin{bmatrix} 2.89 & -0.054 \end{bmatrix}
となる。

同様に、他の点も新座標に変換すると以下のようになる。

X Y X' Y'
1.4 1.65 2.887 -0.054
1.6 1.975 3.266 -0.006
-1.4 -1.775 -1.536 -0.015
-2 -2.525 -2.497 -0.017
-3 -3.95 -4.234 -0.130
2.4 3.075 4.625 0.059
1.5 2.025 3.242 0.103
2.3 2.75 4.309 -0.067
-3.2 -4.05 -4.437 -0.037
-4.1 -4.85 -5.625 0.164

これをX'軸-Y'軸上にプロットすると以下のようになる。


この表とグラフを見ると、新座標(X', Y')では、Y'の値がほとんど0に近くなっていることが分かる。つまり、X'の値だけ教えてもらえば、Y'の値を教えてもらわなくてもほとんどその点がどんな点なのかわかってしまうということになる。ここで、Y'の値を捨ててX'の値だけ使うことを、次元削減と呼ぶ。次元削減を行うには、変換行列を一つ目の固有値のみから作る、つまり
\begin{bmatrix} 0.626194 \\ 0.779667 \end{bmatrix}ととればよい。すると旧座標では2次元だった点が、新座標では1次元に縮退することになる。

参考

もともとの旧座標において、Xの不偏分散 + Yの不偏分散 = Cov(X, X) + Cov(Y, Y) = 16.3755625であった。実は共分散行列の固有値の合計\lambda_1 + \lambda_2も同じ値になっている。