- 1. はじめに
- 2. 前回記事
- 3. サンプルデータについて
- 4. 実装の方針
- 5. ソースコード
- 6. ソースコードの説明
- 7. 出力結果
- 8. 考察
- 9. 第3主成分を求めるコード
- 10. 出力結果(第3主成分を求めるコード)
- 11. まとめ
1. はじめに
前回記事で、第2主成分は「第1主成分と直交する＆分散が最大となる結合係数を求めるということだと理解しています」と書きました。本記事は、その理解で正しいか検証したときの記録です。検証した結果、その理解で正しいということがわかりました。
2. 前回記事
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 の情報が必要です。このように依存関係がありますのでコーディング時に注意します。
5. ソースコード
ソースコードについてですが、前回は行列のシンボルとして sympy.MatrixSymbol を使用していたのですが、今回は sympy.Symbol と np.array を使うことにしました。理由は、前回記事のソースコードでは、Matrixによるキャストが大量発生してしまい、期待していたほど可読性が上がるわけではなかったからです。(MatrixSymbol を使用すると行列の各要素の変数を自動的に作成してくれて便利だったのですが残念です)そのため本記事では MatrixSymbol の代わりに、 Symbol と np.array を使うことにしました。
# -*- coding:utf-8 -*- # MatrixSymbol を使わないで np.array で統一した & Second component 追加したバージョン import numpy as np from sympy import symbols, diff, Eq, solve from sklearn.decomposition import PCA from itertools import chain # 1行目が観測データ1、2行目が観測データ2、3行目が観測データ3 # 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3 とする data = np.array([[3,5,7],[11,13,17],[19,23,29]], dtype=np.float64) ################################################################ # 独自手法で # first component (第一主成分)を求めるための処理 ################################################################ # 観測データ入力用の 3x3 の行列のシンボルを作成 # 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3 Data = np.array([[symbols('Data%d%d'%(i, j)) for j in range(1, 4)] for i in range(1, 4)]) # first component 用のシンボルを作成(サイズ=3x1) Wa = np.array([[symbols('Wa%d1'%i)] for i in range(1, 4)]) # second component 用のシンボルを作成(サイズ=3x1) Wb = np.array([[symbols('Wb%d1'%i)] for i in range(1, 4)]) # second component は first component と直行する(内積=0となる)条件を使用し、1変数を置き換える inner_prod = Wa.T.dot(Wb)[0][0] # solve 関数を使用して Wa と 直行する Wb の Wb[0][0] を求める # solve を使用するとリストが戻るため、要素を取り出すため インデックスを指定して取り出す new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0] # 全要素の数式に含まれる Wb11 を new_Wb11 と置き換える関数を定義 replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11})) # Data を以下のように変換し、次元を1次元に削減する Y = Data.dot(Wa) # shape=(3, 1) Y2 = Data.dot(Wb) # Y の特徴量の分散を求める Ymean = np.mean(Y, axis=0, keepdims=True) Yhat = (Y - Ymean) # 重心を (0, 0, 0) にする、shape=(3, 1) V = Yhat.T.dot(Yhat) # 分散を求める、 shape=(1, 1) # Y2 特徴量の分散を求める Y2mean = np.mean(Y2, axis=0, keepdims=True) Y2hat = (Y2 - Y2mean) # 重心を (0, 0, 0) にする、shape=(3, 1) V2 = Y2hat.T.dot(Y2hat) # 分散を求める、 shape=(1, 1) # 分散 V が最大になるようにする # ここでは、数値解析的に勾配降下法を使用しVの最大値を求めるため微分する 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) #shape=(3, 1) dV2dWb = diff_V2(Wb) #shape=(3, 1) # first components の値を保持する変数を用意しランダムな数値で初期化 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('===========================================') # 第1主成分を学習する for i in range(101): # 勾配を求める dVdWa_val = dVdWa_func(Wa_val) # ハイパーパラメータ b をかけて 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('===========================================') # 第2主成分を学習する 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 を Wa と直交するようにする Wb_val = calc_Wb(replace_Wb11(Wb)).astype(np.float64) # 勾配を求める dV2dW_val = dV2dWb_func(Wb_val, Wa_val) # ハイパーパラメータ b をかけて Wb_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) ################################################################ # sklearn.decomposition.PCA で # first component (第一主成分)を求めるための処理 ################################################################ # solving by sklearn.decomposition.PCA pca = PCA(n_components=3) # 1-3 components を求める pca.fit(data) # 学習する ################################################################ # 独自手法 と sklearn.decomposition.PCA で # 求めた first component (第一主成分)を表示するための処理 ################################################################ print('\n================ 検証 ================') # 独自手法 で求めた first component を表示(比較しやすいようにflattendを使用) print('--- 独自手法 ---') print('第1主成分 = ', Wa_val.flatten()) print('第2主成分 = ', Wb_val.flatten()) # sklearn.decomposition.PCA で求めた frist component を表示 w1_pca, w2_pca, w3_pca = pca.components_ print('--- sklearn.decomposition.PCA ---') print('第1主成分 = ', w1_pca) print('第2主成分 = ', w2_pca) # sklearn.decomposition.PCA で求めた frist component と # コサイン類似度を求める、w1_pca と Wa_val はノルムが１なので内積を # とるだけで cosθ が求まる 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)))) # 第1主成分と第2主成分における分散の比率 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)
6. ソースコードの説明
ソースコードの説明を以下に記載します。
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) #shape=(3, 1) dV2dWb = diff_V2(Wb) #shape=(3, 1)
上記は、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 とします)を求めるコードを以下に記載し終わりとします。
# 連立方程式を解いて 最後の主成分である第3主成分を求める 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) # eq1 かつ eq2 かつ eq3 の条件を満たす Wc の各要素の変数についてシンボリックに解く kais = np.array(solve([eq1, eq2, eq3], *Wc.flatten())) # 求めたシンボリックな解に Wa と Wb の値を代入して実際の値を求める関数を定義 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]) # 連立方程式で求めた第3主成分と sklearn.decomposition.PCA で求めたものとの内積を求めて値が 1 に近いこと(方向が同じであること)を確認する 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の理解を深めました。
以上です。おやすみなさい。
