主成分分析

主成分分析は、機械学習にとって重要な手法ですが、なかなかイメージしにくい部分です。今回はそんな主成分分析の全体像を明らかにした上で、線形代数を利用した詳細なロジックを説明し、実装したいと思います。

主成分分析とは

主成分分析は、教師なし学習により特徴量の数を減らす次元削減手法を指します。

例えば、トレーニングデータも、テストデータも、二次元モデルで十分にターゲットベクトルを予測できるとします。一方、データにはターゲットベクトルの予測に関係のない特徴量が何十個も存在し、トレーニングデータでこれら全ての特徴量を加味した場合、何十次元にわたる複雑なモデルが作られてしまいます。二次元モデルで十分なデータに対し、このようにトレーニングデータにオーバーフィットした複雑なモデルを使った場合、テストデータにおける予測精度は上がりません。

そんなとき、主成分分析はトレーニングデータへのオーバーフィットを解消し、より汎用的で予測精度の高いモデルを実現してくれるかもしれまん。おおよそ以下の3ステップで、主成分分析を行うことができます。

  1. データから分散共分散行列を作成する。
  2. 分散共分散行列から固有値・固有ベクトルを計算する。
  3. 任意の特徴量の数まで、固有値の大きい固有ベクトルを新たな軸(基底)としてデータを変換する。

分散共分散行列

主成分を特定するには、データのふるまいを理解する必要があります。

特徴量x1とx2における以下のようなデータを考えてみましょう。各特徴量ごとのふるまいを見ていくと、x1方向の広がりがx2方向の広がりよりも大きいことがわかります。x1とx2の広がりの関係はどうでしょうか。x1方向の増加に対し、x2方向へ少し増加していることから、x1とx2には緩やかな正の相関があることがわかります。

これらの情報をまとめたものが分散共分散行列です。

納得しにくい部分なので、具体的に計算してみましょう。座標(s, t)を用意し、行列Aにより変換します。すると、(s, t) = (1, 0)のとき(x1, x2) = (9, 4)になり、(s, t) = (0, -1)のとき(x1, x2) = (-4, -3)になります。これらの座標は上図の座標と一致することから、分散共分散行列がデータのふるまいを表していると言えそうです。

\[
\left[\begin{array}{rrr}
9 & 4 \\
4 & 3
\end{array}\right]
\left[\begin{array}{rrr}
s \\
t
\end{array}\right] =
\left[\begin{array}{rrr}
9s + 4t \\
4s + 3t
\end{array}\right]
\]

固有値・固有ベクトル

1.固有値・固有ベクトルとは

上記の通り、分散共分散行列はデータのふるまいを表しているのでした。そうしたデータのふるまいをベクトル(向きと大きさを表す矢印みたいなもの)で表せると嬉しいことがあります。標準座標系(普段使う普通の座標系)に比べ、そのベクトルはより良くデータのふるまいを表しているため、そのベクトルを中心にした座標系を再構築すると、各座標の値が圧縮されるのです。※主成分分析が固有値・固有ベクトルに帰着することを証明するにはラグランジュ未定乗法数を使うのですが、今回は直感的な説明にとどめました。

こうしたベクトルはテキトーに求められるわけではありません。通常ある行列Aでベクトルを変換すると、向きも大きさも異なるベクトルに変換されてしまうからです。例えば以下のように、ベクトルbを行列Aで変換した場合、向きも大きさも異なるベクトルb’に変換されます。ところが、ベクトルvの場合、向きは変わらず大きさだけ異なるベクトルv’に変換されています。

行列Aで変換しても、向きが変わらないベクトル(データのふるまいを保持するようなベクトル)でなければ、標準座標系からそのベクトルを中心にした座標系へ再構築した際、変換前座標と変換後座標の関係が失われてしまいます。従って、求めたいのは行列Aで変換しても、向きが変わらないベクトルです。

ある行列Aで変換しても向きは変わらず大きさだけ変わるベクトルを固有ベクトルvとし、変換の際に変わるベクトルの大きさを固有値λとすると、以下の方程式が成り立ちます。固有値λはデータの散らばり度合いを表しているため、大きいほどデータのふるまいをより良く表す固有ベクトルだと判断できます。

\[
Av = \lambda v
\]

2.固有値の計算

行列Aを以下の分散共分散行列とします。

\[
A =
\left[\begin{array}{rrr}
9 & 4 \\
4 & 3
\end{array}\right]
\]

固有値・固有ベクトルの方程式を、移項してvでくくると、

\[
Av = \lambda v \\
|A-\lambda| v = 0
\]

行列Aのサイズは2×2である一方、固有値λはスカラー(サイズ1×1の値)なのでこのままでは計算できません。単位行列Eは左右いずれから掛けてもベクトルや行列を変えない(サイズが変わるだけ)ため、行列Aのサイズに合うような2×2の単位行列Eを固有値λへ掛けると、

\[
|A-\lambda E| v = 0 \verb| ・・・①|
\]

固有ベクトルvが零ベクトルとなるv=0の場合には関心がないため、

\[
|A-\lambda E| = 0
\]

行列Aと単位行列Eを代入すると、

\[
|A-\lambda E| =
\left[\begin{array}{rrr}
9 & 4 \\
4 & 3
\end{array}\right] –
\lambda\left[\begin{array}{rrr}
1 & 0 \\
0 & 1
\end{array}\right] =
\left[\begin{array}{rrr}
9-\lambda & 4 \\
4 & 3-\lambda
\end{array}\right]
\]

この行列式が0になる時の固有値λを求めるため、因数分解すると、固有値λが求まります。

\[
(9-\lambda)(3-\lambda) -4 \times 4 = (\lambda-11)(\lambda-1) = 0 \\
\lambda = 11, 1
\]

3.固有ベクトルの計算
(i)λ=11のとき

固有値λを①に代入し、固有ベクトルv1の要素をw1とw2で列ベクトル表記すると、

\[
\left[\begin{array}{rrr}
9-\lambda & 4 \\
4 & 3-\lambda
\end{array}\right] v_1 =
\left[\begin{array}{rrr}
-2 & 4 \\
4 & -8
\end{array}\right]
\left[\begin{array}{rrr}
w_1 \\
w_2
\end{array}\right] = 0\\
{\begin{equation}
\left\{\begin{aligned}
-2w_1 + 4w_2 = 0\\
4w_1 – 8w_2 = 0
\end{aligned}\right.
\end{equation}
\verb| ⇔ |
w_1 -2w_2 = 0
}
\]

二つの方程式が一つで表せてしまうことから、解は不定(固有ベクトルv1の要素w1かw2のいずれかを固定しなければもう一方が決まらない)となります。言い換えれば、任意の定数s1(s1 ≠ 0)でw1かw2のいずれか一方を固定すると、もう一方も自動的に決まります。最後にs1でくくると、一つ目の固有ベクトルv1が求まります。s1は任意の定数なので(ベクトルを同じ方向へ伸縮するだけなので)、固有ベクトルv1としては[2 1]が該当します。

\[
w_1 = 2s_1 \verb|とすると、| w_2 = s_1 \\
v_1 = \left[\begin{array}{rrr}
2s_1 \\
s_1
\end{array}\right] =
s_1\left[\begin{array}{rrr}
2 \\
1
\end{array}\right]
\]

(ii)λ=1のとき

固有値λを①に代入し、固有ベクトルv2をw1とw2で列ベクトル表記すると、

\[
\left[\begin{array}{rrr}
9-\lambda & 4 \\
4 & 3-\lambda
\end{array}\right] v_2 =
\left[\begin{array}{rrr}
8 & 4 \\
4 & 2
\end{array}\right]
\left[\begin{array}{rrr}
w_1 \\
w_2
\end{array}\right] = 0 \\
{\begin{equation}
\left\{\begin{aligned}
8w_1 + 4w_2 = 0\\
4w_1 + 2w_2 = 0
\end{aligned}\right.
\end{equation}
\verb| ⇔ |
2w_1 + w_2 = 0
}
\]

任意の定数s2(s2 ≠ 0)でw1かw2のいずれか一方を固定すると、もう一方も自動的に決まります。最後にs2でくくると、二つ目の固有ベクトルv2が求まります。s2は任意の定数なので(ベクトルを同じ方向へ伸縮するだけなので)、固有ベクトルv2としては[1 -2]が該当します。

\[
w_1 = s_2 \verb|とすると、| w_2 = -2s_2 \\
v_2 =
\left[\begin{array}{rrr}
s_2 \\
-2s_2
\end{array}\right] =
s_2\left[\begin{array}{rrr}
1 \\
-2
\end{array}\right]
\]

座標変換

固有ベクトルv1とv2を求めることで、主成分の方向が明らかになりました。ところが、元データの座標はx1とx2軸において作られたものであり(標準基底i = [1 0]とj = [0 1]を基底として作られた標準座標系における座標)、関心があるのはw1とw2軸における座標(固有ベクトルv1とv2を基底として作られた座標系における座標)です。w1とw2軸に焦点を当てて変換後座標を求めるには、以下の方程式をwについて解く必要があります。この式が意味するところは「行列P×変換後座標=変換前座標」であり、少しイメージがつかみにくい部分です。最後に検算をして式の正しさを確認しようと思います。

\[
Pw = x \verb| ・・・②| \\
\]

逆行列Pインバースを両辺へ掛けると、

\[
P^{-1}Pw = P^{-1}x
\]

行列P×逆行列Pインバースは単位行列Eなので、

\[
Ew = P^{-1} x
\]

単位行列Eはベクトルや行列を変えないため、

\[
w = P^{-1} x \verb| ・・・③|
\]

ここで、行列Pは固有ベクトルv1とv2を行列にしたものであるため、上記で求めた固有ベクトルv1とv2を列ベクトルとして行列表記すると、

\[
P = \left[\begin{array}{rrr}
2 & 1\\
1 & -2
\end{array}\right] \\
\]

逆行列Pインバースを求めると、

\[
P^{-1} =
\frac{1}{|P|}\tilde{P} =
\frac{1}{2\times(-2) – 1\times1}
\left[\begin{array}{rrr}
(-1)^{1+1}\times(-2) & (-1)^{1+2}\times1 \\
(-1)^{2+1}\times1 & (-1)^{2+2}\times2
\end{array}\right] =
-\frac{1}{5}
\left[\begin{array}{rrr}
-2 & -1 \\
-1 & 2
\end{array}\right] =
\left[\begin{array}{rrr}
0.4 & 0.2 \\
0.2 & -0.4
\end{array}\right] \verb| ・・・④|
\]

④を③に代入し、xとwを列ベクトル表記すると、x1とx2座標からw1とw2座標を求める式が得られます。

\[
\left[\begin{array}{rrr}
w_1\\
w_2
\end{array}\right] =
\left[\begin{array}{rrr}
0.4 & 0.2 \\
0.2 & -0.4
\end{array}\right]
\left[\begin{array}{rrr}
x_1\\
x_2
\end{array}\right] =
\left[\begin{array}{rrr}
0.4x_1+0.2x_2\\
0.2x_1-0.4x_2
\end{array}\right] \verb| ・・・⑤|
\]

⑤式はグラフ上何を意味しているのでしょう。例えば、標準基底i = [1 0]とj = [0 1]を基底として作られた標準座標系における座標(x1, x2) = (-4, -3)は、標準基底iとjを使い-4i + (-3)jと表せます。また、新基底v1 = [2 1]とv2 = [1 -2]を使い-2.2v1 + 0.4v2とも表せます。つまり、同じ座標を別のモノサシで表しているだけです。

関心があるのはw1とw2軸における座標(固有ベクトルv1とv2を基底として作られた座標系における座標)でした。それなら、標準基底iとjが構成する標準座標系における座標が、新基底v1とv2が構成する新座標系においてどのように表現されるかを考えた方がわかりやすくなります。これが座標変換の意義です。

⑤を試してみると、(x1, x2) = (9, 4)のとき(w1, w2) = (4.4, 0.2)になり、(x1, x2) = (-4, -3)のとき(w1, w2) = (-2.2, 0.4)になります。あわせて②を検算してみると、(w1, w2) = (4.4, 0.2)のとき(x1, x2) = (9, 4)になり、(w1, w2) = (-2.2, 0.4)のとき(x1, x2) = (-4, -3)になるため、正しいことが確認できました。

結果が正しいようなので、x1とx2座標をw1とw2軸におけるw1とw2座標に変換しました。第一主成分であるw1方向へは値が散らばっているため、保持されている情報が多いことがわかります。一方、第二主成分であるw2方向へは値が0へ近づいているため、保持されている情報が少ないことがわかります。今回の二次元データでは、第一主成分であるw1がおおよそのデータのふるまいを表しているようです。

結局主成分分析は何なのかというと、分散共分散行列からデータのふるまいを上手く表すような固有ベクトルを見つけて、その固有ベクトルを新たなモノサシとしてデータを表現し直す処理だと言えます。

実装

主成分分析のロジックは複雑ですが、Pythonでの実装は簡単です。まずライブラリ一式をインポートします。

# In
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

irisデータセットを利用します。

# In
iris = datasets.load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)

データのスケールが異なると主成分分析の精度に影響するため、スケーリングします。

# In
scaler = StandardScaler()
df_scaled = scaler.fit_transform(df)

スケーリング後のデータを主成分分析します。

# In
pca = PCA()
pca.fit(df_scaled)

各主成分がどれくらいデータを表現しているかは、累積寄与率で確認するこができます。第二主成分まででデータの95%が説明できるようです。従って、第三主成分と第四主成分は分類モデルのトレーニングデータに考慮する必要はないかもしれません。※特徴量の数が5個だけなので、トレーニングデータにオーバーフィットする可能性は低く、本当は主成分分析をする意味はあまり無いです。

# In
plt.figure()
plt.plot(np.hstack([0, pca.explained_variance_ratio_]).cumsum(), 'o-')
plt.xticks(range(5))
plt.xlabel('Components')
plt.ylabel('Explained variance ratio')
plt.grid()

主成分と特徴量の相関係数のヒートマップを作りました。主成分分析により座標変換されてしまうため、主成分から特徴量を一意に特定することはできません。ただし、主成分と特徴量の相関係数から、各主成分でどの特徴量がどの程度考慮されたかがわかります。

# In
fig = plt.figure(figsize=(10, 5))
sns.heatmap(pca.components_,
           cmap='coolwarm',
           annot=True,
           xticklabels=[i for i in df.columns],
           yticklabels=["PC{}".format(x+1) for x in range(4)])
plt.show()

第二主成分まででデータのほとんどが説明できるため、例として第一主成分と第二主成分における特徴量の寄与率をグラフ化しました。

# In
plt.figure(figsize=(6, 6))
for x, y, name in zip(pca.components_[0], pca.components_[1], df.columns):
    plt.text(x, y, name)
plt.scatter(pca.components_[0], pca.components_[1], alpha=0.8)
plt.grid()
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()

最後に、学習したモデルで座標変換し、変換前後で比較しました。値が0に近づくほど、データのふるまいを表現するにあたり、その特徴量の持っている情報が少ないことを意味しています。

変換前の特徴量x1~x4では、0から遠い値にも多く散らばっていることから、それぞれが重要な情報を持っていることがわかります。一方、変換後の特徴量w1~w4では、w1とw2においては0から遠い値にも多く散らばっているものの、w3とw4においては0付近に集中しています。従って、w1とw2が重要な情報を持っており、w3とw4はあまり重要な情報を持っていないようです。これは、累積寄与率から得られた知見とも一致しています。

# In
# 学習したモデルでデータを変換
df_transformed = pca.transform(df_scaled)

# 座標変換前後のデータからデータフレームを作成
df_concat = pd.DataFrame(data=np.concatenate((df_scaled, df_transformed), axis=1), 
                         columns=['x1 (sepal length)', 'x2 (sepal width)', 'x3 (petal length)', 'x4 (petal width)', 
                                  'w1 (PC1)', 'w2 (PC2)', 'w3 (PC3)', 'w4 (PC4)'])

# データフレームからグラフを作成
plt.figure(figsize=(16, 8))
for i in range(len(df_concat.columns)):
    plt.subplot(2, 4, i+1)
    plt.xlim([-5, 5])
    plt.xlabel(df_concat.columns[i])
    sns.distplot(df_concat.iloc[:, [i]])

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です