1. はじめに
前回記事で、第2主成分は「第1主成分と直交する＆分散が最大となる結合係数を求めるということだと理解しています」と書きました。本記事は、その理解で正しいか検証したときの記録です。検証した結果、その理解で正しいということがわかりました。
3. サンプルデータについて
前回同様、使用するサンプルデータは 3次元の特徴量をもつ3つのデータとします。
現実にこのサイズのデータで PCA を行うことは稀かと思われますが、PCA の理解には十分だと思います。
4. 実装の方針
第2主成分が第1主成分と直交するというのは第1主成分をベクトル Wa 、 第2主成分をベクトル Wb とすると Wa と Wb の内積が 0 になるということです。 Wa と Wb が以下の時
Wa = [Wa1, Wa2, Wa3]
Wb = [Wb1, Wb2, Wb3]
内積は以下のように定義されます。
inner_product = Wa1*Wa2 + Wa2*Wb2 + Wa3*Wb3
Wa と Wb が直交しているとき、内積が 0 なので
Wa1*Wb1 + Wa2*Wb2 + Wa3*Wb3 = 0
を Wa, Wb は満たします。ちなみに前回記事で第1主成分 Wa は求めることができましたので、 Wa1, Wa2, Wa3 には数値が入ることがわかります。上記の式を変形して Wb1 について解くと
Wb1 = -(Wa2*Wb2 + Wa3*Wb3) / Wa1
となります。 これを Wb の Wb1 に代入すると
Wb = [-(Wa2*Wb2 + Wa3*Wb3) / Wa1, Wb2, Wb3]
と書くことができます(細かいですが Wa1が 0 でないときに限ります)。Wa と Wb が直交するという条件があるので、それを反映させたときの Wb がこの式になります。この Wb をサンプルデータに適用し適用後のデータ間の分散が最大になる Wb が第2主成分であると理解していますので、その Wb と sklearn.decomposition.PCA を使用して求めた第2主成分を比較し、同じ方向を向いているか検証します。
※当たり前ですが Wb を求めるときに Wa2, Wa3 すなわち Wa の情報が必要です。このように依存関係がありますのでコーディング時に注意します。
ソースコードについてですが、前回は行列のシンボルとして sympy.MatrixSymbol を使用していたのですが、今回は sympy.Symbol と np.array を使うことにしました。理由は、前回記事のソースコードでは、Matrixによるキャストが大量発生してしまい、期待していたほど可読性が上がるわけではなかったからです。(MatrixSymbol を使用すると行列の各要素の変数を自動的に作成してくれて便利だったのですが残念です)そのため本記事では MatrixSymbol の代わりに、 Symbol と np.array を使うことにしました。
import numpy as np
from sympy import symbols, diff, Eq, solve
from sklearn.decomposition import PCA
from itertools import chain
data = np.array([[3,5,7],[11,13,17],[19,23,29]], dtype=np.float64)
Data = np.array([[symbols('Data%d%d'%(i, j)) for j in range(1, 4)] for i in range(1, 4)])
Wa = np.array([[symbols('Wa%d1'%i)] for i in range(1, 4)])
Wb = np.array([[symbols('Wb%d1'%i)] for i in range(1, 4)])
inner_prod = Wa.T.dot(Wb)[0][0]
new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0]
replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11}))
Y = Data.dot(Wa)
Y2 = Data.dot(Wb)
Ymean = np.mean(Y, axis=0, keepdims=True)
Yhat = (Y - Ymean)
V = Yhat.T.dot(Yhat)
Y2mean = np.mean(Y2, axis=0, keepdims=True)
Y2hat = (Y2 - Y2mean)
V2 = Y2hat.T.dot(Y2hat)
diff_V = np.vectorize(lambda x: diff(V[0, 0], x))
diff_V2 = np.vectorize(lambda x: diff(V2[0, 0], x))
dVdWa = diff_V(Wa)
dV2dWb = diff_V2(Wb)
np.random.seed(51)
Wa_val = np.random.random((3,1))
Wb_val = np.random.random((3,1))
b = 0.05
dVdWa_func = lambda Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
chain(data.flatten(), Wa_val.flatten()))))
for w in dVdWa], dtype=np.float64).reshape((3, 1))
dV2dWb_func = lambda Wb_val, Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
for w in replace_Wb11(dV2dWb)], dtype=np.float64).reshape((3, 1))
print('===========================================')
print('================ 第1主成分 ================')
print('===========================================')
for i in range(101):
dVdWa_val = dVdWa_func(Wa_val)
Wa_val += b*dVdWa_val
Wa_val /= np.linalg.norm(Wa_val)
if i%10==0:
print('================ 第1主成分 学習 i = ', i, "================")
v = V[0, 0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
chain(data.flatten(), Wa_val.flatten()))))
print("分散V = ", v)
print("Wa = ", Wa_val.T)
print('===========================================')
print('================ 第2主成分 ================')
print('===========================================')
calc_Wb = np.vectorize(lambda x: x.subs(dict(zip(chain(Wb.flatten(), Wa.flatten()),
chain(Wb_val.flatten(), Wa_val.flatten())))))
for i in range(101):
Wb_val = calc_Wb(replace_Wb11(Wb)).astype(np.float64)
dV2dW_val = dV2dWb_func(Wb_val, Wa_val)
Wb_val += b*dV2dW_val
Wb_val /= np.linalg.norm(Wb_val)
if i%10==0:
print('================ 第2主成分 学習 i = ', i, "================")
v2 = V2[0, 0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
print("分散V2 = ", v2)
print("Wb = ", Wb_val.T)
pca = PCA(n_components=3)
pca.fit(data)
print('\n================ 検証 ================')
print('--- 独自手法 ---')
print('第1主成分 = ', Wa_val.flatten())
print('第2主成分 = ', Wb_val.flatten())
w1_pca, w2_pca, w3_pca = pca.components_
print('--- sklearn.decomposition.PCA ---')
print('第1主成分 = ', w1_pca)
print('第2主成分 = ', w2_pca)
print('--- コサイン類似度 ---')
cos_theta1 = w1_pca.dot(Wa_val)
cos_theta2 = w2_pca.dot(Wb_val)
print(' 第1主成分 cosθ = {} ({}°)'.format(*cos_theta1, *np.rad2deg(np.arccos(cos_theta1))))
print(' 第2主成分 cosθ = {} ({}°)'.format(*cos_theta2, *np.rad2deg(np.arccos(cos_theta2))))
print('--- 第1主成分と第2主成分における分散の比率(sklearn.decomposition.PCA) ---')
print('pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1] \n =', pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1])
print('--- 第1主成分と第2主成分における分散の比率(独自手法) ---')
print('V / V2 \n = ', v / v2)
ソースコードの説明を以下に記載します。
diff_V = np.vectorize(lambda x: diff(V[0, 0], x))
diff_V2 = np.vectorize(lambda x: diff(V2[0, 0], x))
dVdWa = diff_V(Wa)
dV2dWb = diff_V2(Wb)
上記は、V の Wa に関する微分と V2 の Wb に関する微分を求めています。微分を行う関数である diff 関数をブロードキャストできるようにするため np.vectorize 関数を使用しています。
new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0]
上記は、 Eq(inner_prod) で inner_prod = 0 となる数式を構築し、 solve(数式, Wb[0, 0]) で 数式をWb[0, 0] について解いています。solve 関数は解が複数ある場合に備えてリストで解を戻します。今回は、数式をみれば解が一つであることは自明なので、インデックスに 0 を指定してリストから要素をを取り出しています。ここで求めた new_Wb11 を使用して Wb の要素 Wb_11 を置き換えます。
replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11}))
上記は、全要素の数式に含まれる Wb11 を new_Wb11 と置き換える関数を定義しています。
dVdWa_func = lambda Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
chain(data.flatten(), Wa_val.flatten()))))
for w in dVdWa], dtype=np.float64).reshape((3, 1))
上記は Wa の勾配を求める関数を定義しています。シンボルに値を代入するために、key にシンボルを value に シンボルに対応する値を入れた辞書を渡しています。(TensorFlow がわかる人には placeholder に値を代入しているというと伝わるかもしれません。)また、前回記事では、微分をするのにjacobian を使用していましたが今回は for文で回して、各要素を微分しています。
dV2dWb_func = lambda Wb_val, Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
for w in replace_Wb11(dV2dWb)], dtype=np.float64).reshape((3, 1))
上記は Wb の勾配を求める関数を定義しています。 Wb を求めるときに Wa の情報が必要でコーディング時に注意と冒頭で説明しましたが、 上記がそのコードになります。 dict の中身をみると Wb だけでなく Wa の情報も使用していることがわかります。 また、 dV2dWb に直交条件を反映させるため replace_Wb11 関数を適用しています。これにより、 Wa に直交する Wb における分散の勾配を求めることができます。
7. 出力結果
出力結果は以下になります。
===========================================
================ 第1主成分 ================
===========================================
================ 第1主成分 学習 i = 0 ================
分散V = 532.873871148815
Wa = [[0.50212335 0.54326168 0.67285874]]
================ 第1主成分 学習 i = 10 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 20 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 30 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 40 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 50 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 60 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 70 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 80 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 90 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i = 100 ================
分散V = 533.003128892760
Wa = [[0.48958803 0.55232128 0.67471829]]
===========================================
================ 第2主成分 ================
===========================================
================ 第2主成分 学習 i = 0 ================
分散V2 = 0.282286132046230
Wb = [[-0.85629139 0.15854369 0.49155768]]
================ 第2主成分 学習 i = 10 ================
分散V2 = 0.303318086521561
Wb = [[-0.86855083 0.24061133 0.43327317]]
================ 第2主成分 学習 i = 20 ================
分散V2 = 0.315596549516175
Wb = [[-0.87190187 0.30165001 0.38573876]]
================ 第2主成分 学習 i = 30 ================
分散V2 = 0.322411743461676
Wb = [[-0.87100083 0.34594949 0.3488216 ]]
================ 第2主成分 学習 i = 40 ================
分散V2 = 0.326088827550761
Wb = [[-0.86849373 0.37775612 0.32096564]]
================ 第2主成分 学習 i = 50 ================
分散V2 = 0.328042467143946
Wb = [[-0.8656847 0.40051315 0.30029855]]
================ 第2主成分 学習 i = 60 ================
分散V2 = 0.329071954623274
Wb = [[-0.86312883 0.41679385 0.28511667]]
================ 第2主成分 学習 i = 70 ================
分散V2 = 0.329612106418076
Wb = [[-0.86100686 0.42845606 0.27403029]]
================ 第2主成分 学習 i = 80 ================
分散V2 = 0.329894869159100
Wb = [[-0.85933015 0.43682358 0.26596403]]
================ 第2主成分 学習 i = 90 ================
分散V2 = 0.330042715591538
Wb = [[-0.85804391 0.44283655 0.26010852]]
================ 第2主成分 学習 i = 100 ================
分散V2 = 0.330119970975593
Wb = [[-0.85707561 0.44716326 0.25586408]]
================ 検証 ================
--- 独自手法 ---
第1主成分 = [0.48958803 0.55232128 0.67471829]
第2主成分 = [-0.85707561 0.44716326 0.25586408]
--- sklearn.decomposition.PCA ---
第1主成分 = [0.48958803 0.55232128 0.67471829]
第2主成分 = [ 0.85440069 -0.45835872 -0.24475855]
--- コサイン類似度 ---
第1主成分 cosθ = 1.0 (0.0°)
第2主成分 cosθ = -0.9998720868190456 (179.0835687095191°)
--- 第1主成分と第2主成分における分散の比率(sklearn.decomposition.PCA) ---
pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1]
= 1614.1609966447247
--- 第1主成分と第2主成分における分散の比率(独自手法) ---
V / V2
= 1614.57402082519
8. 考察
出力結果を見ると、第1主成分は前回と同じで、問題なく動作していることが確認できます。第2主成分については i = 0 のとき分散が 0.28 で i = 100 のとき 0.33 と増加していることがわかりますので、学習はうまく進んだといえます。検証の箇所で sklearn.decomposition.PCA で求めた第2主成分とのコサイン類似度を出していますがおおよそ cosθ = -1 となりました。これは、 ベクトルが反対方向を向いているということで、 分散の計算でこのベクトルの符号は打ち消されますので、この方向でも反対方向でも分散が最大なのは変わりません。ちなみにどちらの方向になるかは Wb_val の初期状態できまります。最後に、 第2主成分の分散の値が 0.33 と少なく思えるため、分散の計算があっているか確かめるのに、 第1主成分の分散と第2主成分の分散の比率を sklearn.decomposition.PCA で導出したものと比較しています。これにより、分散の比率は 1614 とほとんど同じですので、分散の計算はあっていることがわかります。このサンプルデータにおいては第2主成分は分散が少ないですので、ほとんど意味がない情報といえそうです。以上にて、第2主成分について私の理解について問題ないことが確認できました。最後に、 3x3 のサンプルデータにおける最後の主成分の第3主成分ですが、これは第1主成分と第2主成分に直行するノルム1のベクトルですので連立方程式をとくと解が求まります。連立方程式を解いて第3主成分( Wc とします)を求めるコードを以下に記載し終わりとします。
9. 第3主成分を求めるコード
連立方程式を解いて第3主成分( Wc とします)を求めるコードを以下に記載し終わりとします。
Wc = np.array([[symbols('Wc%d1'%i)] for i in range(1, 4)])
wc_wa_inner = Wc.T.dot(Wa)[0, 0]
wc_wb_inner = Wc.T.dot(Wb)[0, 0]
wc_norm = Wc.T.dot(Wc)[0, 0]
eq1 = Eq(wc_wa_inner)
eq2 = Eq(wc_wb_inner)
eq3 = Eq(wc_norm - 1)
kais = np.array(solve([eq1, eq2, eq3], *Wc.flatten()))
subs = np.vectorize(lambda x: x.subs(dict(zip(chain(Wa.flatten(), Wb.flatten()),chain(Wa_val.flatten(), Wb_val.flatten())))))
print('第3主成分 = ', subs(kais)[0])
cos_theta3 = float(np.inner(subs(kais)[0], w3_pca))
print('--- コサイン類似度 ---')
print(' 第3主成分 cosθ = {} ({}°)'.format(cos_theta3, np.rad2deg(np.arccos(cos_theta3))))
10. 出力結果(第3主成分を求めるコード)
第3主成分 = [0.160390052085319 0.703552574067647 -0.692306873225203]
--- コサイン類似度 ---
第3主成分 cosθ = 0.9998720868190454 (0.9164312904816843°)
第3主成分についても、コサイン類似度をみると、少し誤差はあるものの、sklearn.decomposition.PCA のものとほとんど同じであることが確認できました。誤差を少なくするには、勾配降下法で第1主成分と第2主成分を求めるときに学習率を大きくするなどの工夫が必要です。
11. まとめ
3x3 のサンプルデータについて第1主成分と第2主成分を数値解析的に求め、第3主成分を代数的に求めました。それぞれ sklearn.decomposition.PCA で求めた値ととほぼ同じ値になることを確認しました。PCAの理解を深めました。
以上です。おやすみなさい。
