主成分分析（PCA：Principal Component Analysis）とはデータの散らばりが最も大きい方向へ次元圧縮（射影）する機械学習手法である.
例えば３次元空間上のデータを可視化するために2次元平面へ射影したデータを取得したいときや, 多くの特徴量をもつデータからより特徴を捉えた成分を分析したいときなどに役立つ.
数学的な理論
以下ではまず, １次元空間への射影を考えたときに, 最も散らばり（分散）が大きくなる方向を考える.
データ \({\boldsymbol x}_i\; (i=1,\ldots, n)\) の 最も分散の大きい方向を表す単位ベクトルを \({\boldsymbol w}\) とする. このとき, \(\boldsymbol{w}\) は共分散行列 \[S=\frac{1}{n} \sum_{i=1}^n (\boldsymbol{x}_i-\bar{\boldsymbol{x}}) (\boldsymbol{x}_i-\bar{\boldsymbol{x}})^T\] の固有値最大の固有ベクトルである. （ \(\bar{\boldsymbol{x}}=\frac{1}{n} \sum_{i=1}^n \boldsymbol{x}_i\) はデータの平均）
これは次の筋道を辿れば分かる.
ここでは最も分散の大きい1つの方向だけを考えたが, 一般に固有値の大きい方から m個の固有ベクトルをとれば、その方向へ射影してm次元空間へ次元圧縮することができる. このとき最も固有値の大きい固有ベクトルから順に, 第一主成分, 第二主成分などという. \(S\) は対称行列だから, これらは全て直交する.
Pythonによる実装
実際に実装するには, 主に2つの方法が考えられるが, 2つ目の方法を用いることが多い.
- 分散行列 \(S\) は対称行列なので, 直交行列による対角化を考える.
- \(A= \begin{pmatrix} (\boldsymbol{x}_1-\bar{\boldsymbol{x}})^T \\ \vdots \\ (\boldsymbol{x}_n-\bar{\boldsymbol{x}})^T\end{pmatrix}\) とおき, \(A\) の特異値分解を考える.
特異値分解や特異ベクトルなどの用語については以下の記事を参照.
1つ目の方法については上の考え方通りである.
2つ目の方法については \(S=\frac{1}{n} A^T A\) であり, \(A\) の右特異ベクトルは \(S\) の固有ベクトルになるので, 1と同じ結果を得ることができる. 以下において, この2つ目の方法で実装する際の手順を示す.
- \(A= \begin{pmatrix} (\boldsymbol{x}_1-\bar{\boldsymbol{x}})^T \\ \vdots \\ (\boldsymbol{x}_n-\bar{\boldsymbol{x}})^T\end{pmatrix}\) とおく.
- \(A\) の特異値分解 \(A= U \Sigma V^T\) を計算する.
- m次元空間へ射影したいときは, \(V\) の1列目からm列目までの列ベクトルへデータを射影する
実際に実装すると以下のようになる.
import numpy as np
from scipy.linalg import svd
def mypca(X, n_components):
A = X-X.mean(axis=0)
U, S, VT = svd(A)
m = n_components
return np.dot(A,VT[:m].T) #np.dot(X,VT[:m].T)でもよいが, Aを射影すれば平均を0にできる
#Xは形状が(データ数,特徴量の数)のnumpy配列
#数式で書いた場合と転置の関係になっているので注意
#Aは数式での説明と一致する例として数学、理科、国語、社会における10人分のテストの点数を主成分分析して2次元に次元削減してみる.
math=np.array([83,95,32,58,47,75,61,20,66,49])
science=np.array([88,90,48,48,61,78,53,33,51,62])
japanese=np.array([82,51,91,51,78,44,57,41,63,55])
social_studies=np.array([90,43,98,63,85,33,67,21,70,62])
X=np.vstack((math,science,japanese,social_studies)).T
X_new=mypca(X,n_components=2)
print(X_new)
"""出力
[[-3.80385501e+01 3.14644645e+01]
[ 1.68539450e+01 4.87078431e+01]
[-4.19439364e+01 -3.43079247e+01]
[ 7.39416289e+00 -7.74348465e+00]
[-2.64705716e+01 -1.24466633e+01]
[ 3.15134096e+01 2.71713856e+01]
[ 4.14490372e-02 -3.16003795e+00]
[ 5.15254482e+01 -4.20553957e+01]
[-5.77822534e+00 -1.20971508e+00]
[ 4.90286860e+00 -6.42047185e+00]]
"""散布図を描いてみると次のようになる. 全体的に散らばっているように見える.
import matplotlib.pyplot as plt
plt.scatter(X_new[:,0],X_new[:,1])
plt.show()
第一主成分, 第二主成分を見てみると, それぞれ文系科目と理系科目の値（絶対値）が大きいので, 「文系学力」と「理系学力」の2成分でデータを捉えることができる.
A = X-X.mean(axis=0)
U, S, VT = svd(A)
print(VT)
"""出力
[[-0.05761705 -0.11743878 -0.54797673 -0.8262021 ]
[ 0.78156567 0.61240935 -0.07680933 -0.09061027]
[ 0.53568137 -0.68517325 -0.38152664 0.31308229]
[ 0.31445326 -0.37643133 0.74044531 -0.45952078]]
"""
コメント