PCAをもうちょっと理解した上で使ってみる
恥ずかしながら、パッケージユーザーなもので、PCAがきちっとwかっていない状態で使おうとしてハマっていた感があった。
そこで、「何をしているか?」を明確にして、扱ってみることにする。
ちなみに、この記事で扱ってるscikit-learnはバージョンが0.15です。
他のバージョンだと、メソッドが存在してなかったりするんで、ご注意を。
scikit-learnでやってみる
まず、scikit-learnの基本的な使い方にそって、やってみる。
In [55]: a = np.array([[ 0, 1, 2, 3, 4, 2, 4, 1],[ 2, 1, 4, 2, 4, 5, 6, 7],[ 8, 2, 1, 2, 1 ,9, 10, 11]])
In [55]: pca = PCA(a, n_components=8)
In [55]: pca.fit(a)
In [56]: pca.components_
Out[56]:
array([[ 0.50488641, 0.06499704, -0.08008693, -0.05745209, -0.19499111,
0.43234443, 0.37489233, 0.60470071],
[ 0.3453956 , 0.1116729 , -0.65966052, 0.16232091, -0.33501871,
-0.04027111, 0.1220498 , -0.52723383],
[ 0.40199019, 0.14317705, -0.06778525, -0.19504746, 0.18034548,
0.27022208, -0.81834129, -0.00624017]])
In [57]: pca.components_.shape
Out[57]: (3, 8)
In [58]: pca.transform(a)
Out[58]:
array([[ -7.57714864e+00, 1.42522380e+00, -2.22044605e-16],
[ -9.95075408e-01, -2.40574369e+00, 1.59594560e-15],
[ 8.57222405e+00, 9.80519892e-01, -6.21031004e-16]])
In [59]: pca.transform(a).shape
Out[59]: (3, 3)
In [272]: pca.explained_variance_
Out[272]: array([ 4.39621273e+01, 2.92676162e+00, 3.36224570e-31])
scikit-learnの一部だけ流用してやってみる
では、自分できちんとやってみよう。
PCAの手順はそんなに難しいことでない。手順的には
1 入力行列(n * m)の共分散行列(m * m)を作成する
2 共分散行列を固有値分解する。固有値(1 * m)と固有ベクトル(m * m)に分離できる
3 固有値の高い順に欲しい次元数だけ、固有値に対応する固有ベクトルを取る。
4 とった固有ベクトルで入力行列を写像する。(n * m と m * m の内積計算)
これをscikit-learnを使ってやってみよう。
In [61]: cov_matrix = pca.get_covariance() # 共分散行列の生成
In [13]: import numpy as np
In [14]: [E, V] = np.linalg.eig(cov_matrix) # 固有値分解
Out[16]:
array([[ 0.36158968, -0.65653988, -0.58099728, 0.31725455],
[-0.08226889, -0.72971237, 0.59641809, -0.32409435],
[ 0.85657211, 0.1757674 , 0.07252408, -0.47971899],
[ 0.35884393, 0.07470647, 0.54906091, 0.75112056]])
これで、共分散行列の固有値分解ができた。
まず、主成分得点に対応する固有値から確認してみる。
In [67]: E
Out[67]:
array([ 4.39621273e+01 +0.00000000e+00j,
2.92676162e+00 +0.00000000e+00j,
1.76593839e-15 +3.86799686e-15j,
1.76593839e-15 -3.86799686e-15j,
-2.29569145e-15 +0.00000000e+00j,
-1.29253764e-16 +3.91441638e-16j,
-1.29253764e-16 -3.91441638e-16j, 5.35881238e-18 +0.00000000e+00j])
In [102]: projection_matrix = V[0:2,] # 固有ベクトルから2次元までの空間をとってくる
In [103]: projection_matrix
array([[-0.50488641+0.j , 0.34539560+0.j ,
0.19411192+0.41056554j, 0.19411192-0.41056554j,
0.07143826+0.j , -0.01791936-0.00421049j,
-0.01791936+0.00421049j, -0.00091420+0.j ],
[-0.06499704+0.j , 0.11167290+0.j ,
-0.76012933+0.j , -0.76012933-0.j ,
0.20260054+0.j , 0.03205764+0.06085498j,
0.03205764-0.06085498j, 0.00976708+0.j ]])
ところで、思い出していただきたいのは、scikit-learnの演算結果である。
scikit-learnで演算した結果で、主成分(つまり固有値)を表示すると、
In [272]: pca.explained_variance_
Out[272]: array([ 4.39621273e+01, 2.92676162e+00, 3.36224570e-31])
おいおい、自分で算出してみたやつと違うで!
pca.transfrom()メソッドの結果をみても
In [59]: pca.transform(a).shape
Out[59]: (3, 3)
勝手に、3次元に落とされている。
どうやら、scikit-learnは累積寄与率ってので、勝手に切った結果を返しているらしい。
まぁ、ちょっと試すくらいならいいんだろうが、あまり自由度は高くないようだ。
自分でやった部分が正しいのか?確認してみる
一応、この固有値と固有ベクトルが正しいのか?の確認をしてみる。
と、その前に、固有値分解の基本をちょっとおさらいしてみる。
この時、Ax = λx
を満たすλとxを見つければよい。
xが0にならないためには、| λE - A | = 0
の固有方程式を解けばいい。
この結果、スカラ値のλとベクトルのxが算出される。
行列演算で行うときには、行列Aの固有ベクトルを束ねた行列をX、固有値を対角成分にもつ行列をDとした時、
AX = XD
となるので、A = XDX^(-1)
とすると、元の行列Aに復元できる。
ここまでの参考はこのサイト
なので、この検算をしてみよう。
In [332]: cov_matrix
Out[332]:
array([[ 11.55555556, 1.55555556, -2.44444444, -1.11111111,
-4.66666667, 9.55555556, 8.44444444, 12.88888889],
[ 1.55555556, 0.22222222, -0.44444444, -0.11111111,
-0.66666667, 1.22222222, 1.11111111, 1.55555556],
[ -2.44444444, -0.44444444, 1.55555556, -0.11111111,
1.33333333, -1.44444444, -1.55555556, -1.11111111],
[ -1.11111111, -0.11111111, -0.11111111, 0.22222222,
0.33333333, -1.11111111, -0.88888889, -1.77777778],
[ -4.66666667, -0.66666667, 1.33333333, 0.33333333,
2. , -3.66666667, -3.33333333, -4.66666667],
[ 9.55555556, 1.22222222, -1.44444444, -1.11111111,
-3.66666667, 8.22222222, 7.11111111, 11.55555556],
[ 8.44444444, 1.11111111, -1.55555556, -0.88888889,
-3.33333333, 7.11111111, 6.22222222, 9.77777778],
[ 12.88888889, 1.55555556, -1.11111111, -1.77777778,
-4.66666667, 11.55555556, 9.77777778, 16.88888889]])
In [333]: origin_matrix = np.dot(np.dot(V ,np.diag(E)), np.linalg.inv(V)) # A = XDX^(-1)の演算
In [334]: cov_matrix.astype("int32") == origin_matrix.astype("int32") # 型が違ったぽいので、統一して比較演算
Out[334]:
array([[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True],
[ True, True, True, True, True, True, True, True]], dtype=bool)
同じになった。
まとめ
scikit-learnの一部だけ使って、PCAをやってみたよ。
scikit-learnには、好きなだけ次元数を残す。っていう機能がないらしいので、自分でやるの大切。
でも、scikit-learn使ったとこは、共分散行列を算出しただけだから、実は無理してscikit-learnを使う必要もなかった・・・という話。