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は累積寄与率ってので、勝手に切った結果を返しているらしい。

まぁ、ちょっと試すくらいならいいんだろうが、あまり自由度は高くないようだ。

自分でやった部分が正しいのか?確認してみる

一応、この固有値固有ベクトルが正しいのか?の確認をしてみる。

と、その前に、固有値分解の基本をちょっとおさらいしてみる。

Aを元の行列、xを固有ベクトル、λを固有値とする。

この時、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を使う必要もなかった・・・という話。