2. はじめに
本記事は、主成分分析(以下PCA)について自分の理解があっているか確かめた時の記録です。PCAの第1主成分(First component)について、自分の理解があっているか検証を行いました。Sympy で独自実装した数値解析のアルゴリズムによって求めた第1主成分と sklearn.decomposition.PCA によって求めたものを比較し、誤差が少ないことを確認しました。
2019/07/04 更新:
続編は以下です。
(第2主成分、第3主成分についても検証しました)
keimina.hatenablog.jp
4. PCAの第一主成分について自分が理解していること
私の理解している、PCAの第一主成分の求め方を書きます。例えば3次元の特徴量をもつ複数のデータがあったとして、それらのデータは 点(x, y, z) のように３次元で表現できます。(説明を簡略化するため、これらの点の集合の重心は (0, 0, 0) とします。)点の集合の要素のi番目の特徴量を、x[i], y[i], z[i] としその線形結合を Y[i](次元削減後の特徴量) とし、線形結合の係数を w1, w2, w3 とすると
Y[i] = w1*x[i] + w2*y[i] + w3*z[i]
とかけます。
ここで、例えば、 データが２つあったとして、
data[0] = (1, 2, 3) data[1] = (4, 5, 6)
とします。すると Y は、
Y[0] = w1*1 + w2*2 + w3*3 Y[1] = w1*4 + w2*5 + w3*6
となります。当たり前ですが、以下のような変換が行われることがわかります。
data[0] (3次元) → Y[0] (スカラ) data[1] (3次元) → Y[1] (スカラ)
この変換後の Y を使用し以下の V を定義します。
V = Y[0]*Y[0] + Y[1]*Y[1]
つまり、 V は Y[i] の分散です。
そして V が最大になるような w1, w2, w3 が第1主成分です。
5. 検証の方針
上記の理解で本当に正しいのか検証するため、以下の方針に従い実装しました。
- V を定式化する
- 勾配降下法により V が最大になるような第1主成分(線形結合の係数)を求める
- sklearn.decomposition.PCA により第1主成分を求める
- 2. と 3. で求めた第1主成分を比較し差分が少ないことを確認する
6. ソースコード
方針に従い実装したソースコードは以下です。
※実装には主に Sympy を使用しています。(TensorFlowでもできますが、Sympy だとprint文で数式の表示ができデバッグが楽なため Sympy を使用しました。)
import numpy as np from sympy import MatrixSymbol, Matrix from sklearn.decomposition import PCA # 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 (第一主成分)を求めるための処理 ################################################################ # 特徴量を定義する # 1行目が観測データ1、2行目が観測データ2、3行目が観測データ3 # 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3 Data = MatrixSymbol('Data', 3, 3) # 観測データ入力用の 3x3 の行列のシンボルを作成 Data_a = np.array(Data) # numpy 表現にする # 1行3列の行列(縦ベクトル)を作成(first components の変数) W = MatrixSymbol('W', 3, 1) # first components 入力用の 3x1 の行列のシンボルを作成 W_a = np.array(Matrix(W)) # numpy 表現にする # Data_a を以下のように変換し、次元を1次元に削減する Y = Data_a.dot(W_a) # shape=(3, 1) # 変換後の特徴量の分散を求める 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) # 分散 V が最大になるようにする # ここでは、数値解析的に勾配降下法を使用しVの最大値を求めるため微分する # Matrix(V).diff(W[0, 0]) # Matrix(V).diff(W[1, 0]) # Matrix(V).diff(W[2, 0]) # 上記を求めるのに jacobian() というものが便利だったため使用しました。 dVdW = Matrix(V).jacobian(W).T # shape=(3, 1) dVdW_a = np.array(dVdW) # numpy 表現にする # first components の値を保持する変数を用意しランダムな数値で初期化 np.random.seed(51) # 他の人がやっても同じ結果になるようにする W_a_val = np.random.random((3,1)) # ハイパーパラメータ(学習係数)を定義 b = 0.00003 # 勾配降下法により first component の最大値を求める。 for i in range(101): # 勾配を求める dVdW_a_val = np.array(dVdW.subs({W:Matrix(W_a_val), Data:Matrix(data)})).astype(np.float64) # ハイパーパラメータ b をかけて W_a_val を更新する W_a_val += b*dVdW_a_val # 正規化する W_a_val /= np.linalg.norm(W_a_val) #np.sum(W_a_val.T.dot(W_a_val)) でも同じ if i%10==0: print('学習 i = ', i) print("分散=", Matrix(V).subs({W:Matrix(W_a_val), Data:Matrix(data)})) print("W=", Matrix(W_a_val)) ################################################################ # 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 ================') # 独自手法 で求めた first component を表示(比較しやすいようにflattendを使用) print('--- 独自手法 ---') print(W_a_val.flatten()) # sklearn.decomposition.PCA で求めた frist component を表示 w1_pca = pca.components_[0] print('--- sklearn.decomposition.PCA ---') print(w1_pca) # 差分 を表示 print('--- 差分 ---') print(w1_pca - W_a_val.flatten()) # sklearn.decomposition.PCA で求めた frist component と # コサイン類似度を求める、w1_pca と W_a_val はノルムが１なので内積を # とるだけで cosθ が求まる print('--- コサイン類似度 ---') cos_theta = w1_pca.dot(W_a_val) print('cosθ = {} rad ({}°)'.format(*cos_theta, *np.rad2deg(np.arccos(cos_theta))))
ソースコードの解説についてはソースコードのコメントをみてください。最後に誤差が少ないか確認するため、コサイン類似度を求めています。コサイン類似度については以下を参考にしてください。
7. 出力結果
学習 i = 0 分散= Matrix([[326.893552407606]]) W= Matrix([[0.885069441686711], [0.0711731365114243], [0.459985291103423]]) 学習 i = 10 分散= Matrix([[398.915653114400]]) W= Matrix([[0.828195710513823], [0.186379419417903], [0.528540043046834]]) 学習 i = 20 分散= Matrix([[452.012116765482]]) W= Matrix([[0.765356021528304], [0.282142685524890], [0.578468378837745]]) 学習 i = 30 分散= Matrix([[486.534391735230]]) W= Matrix([[0.705988103517652], [0.356151889455055], [0.612157356673220]]) 学習 i = 40 分散= Matrix([[507.182771621072]]) W= Matrix([[0.655167743479454], [0.410912766518736], [0.633960508404080]]) 学習 i = 50 分散= Matrix([[518.921766869242]]) W= Matrix([[0.614225060024774], [0.450555959412570], [0.647863336727261]]) 学習 i = 60 分散= Matrix([[525.403750244901]]) W= Matrix([[0.582431990459451], [0.479009309986766], [0.656751899453244]]) 学習 i = 70 分散= Matrix([[528.925385131233]]) W= Matrix([[0.558289211467764], [0.499398761609065], [0.662505874133985]]) 学習 i = 80 分散= Matrix([[530.821835476805]]) W= Matrix([[0.540205689228539], [0.514032980713747], [0.666294160310340]]) 学習 i = 90 分散= Matrix([[531.838240751787]]) W= Matrix([[0.526776505606956], [0.524565081220116], [0.668833304123724]]) 学習 i = 100 分散= Matrix([[532.381592136768]]) W= Matrix([[0.516858344804962], [0.532166525145605], [0.670564121408182]]) ================ first component ================ --- 独自手法 --- [0.51685834 0.53216653 0.67056412] --- sklearn.decomposition.PCA --- [0.48958803 0.55232128 0.67471829] --- 差分 --- [-0.02727032 0.02015476 0.00415417] --- コサイン類似度 --- cosθ = 0.9994164292655139 (1.9575168306292177°)
出力結果を見ると 「学習 i = 0」の時、分散が326だったのが、「学習 i = 100」で 532 に増えていることがわかります。最大化の方向に学習がうまく進んでいるようです。最後に「--- 独自手法 ---」に独自手法で求めた第1主成分と「--- sklearn.decomposition.PCA ---」にsklearn.decomposition.PCA で求めた第1主成分を表示しています。若干の誤差はあるもののほぼ等しいことがわかります。等しさの指標としてベクトルの差分とコサイン類似度も表示しました。コサイン類似度を見るとベクトルの方向がほとんど一致していることがわかります。
※学習の経過を見るため、学習係数bはあえて小さくしています。学習係数bを大きくするとすぐに sklearn.decomposition.PCAで求めた値と同じ値へ収束します。
8. まとめ・考察
PCAの第1主成分(First component)について、自分の理解があっているか検証しました。PCA ではデータの特徴量を線形変換し次元削減した後それらの分散が最大となるような線形変換の結合係数を求めているということがわかりました。ただ、PCAは特徴量に線形変換を施すことが前提としてありますので、特徴量に非線形の関係性がある場合は、たとえ分散を最大化しても、元の特徴をうまく表せるとは限らないということも言えそうです。PCA の理解が深まりました。後は第2主成分以降の成分なのですが、これについては、第1主成分と直交する＆分散が最大となる結合係数を求めるということだと理解しています。余裕があれば、こちらも記事にしたいと思います。
以上です。おやすみなさい。
