TL;DR
特徴次元数削減アルゴリズムの代表格といえば主成分分析 (PCA : Principal Component Analysis) ですね。ただ、やってることは感覚的に何となく分かっていても、線形代数的に何やってるのかをよく覚えておらず、気になったので復習してみた記事です。
PCAを手計算で実装し、結果をscikit-learnのPCAサブパッケージのものと比較します。（手計算といっても sklearn.decomposition.PCA を使わないでという意味であってNumPyは使います）
主成分分析の概要
高次元の特徴からなるデータの次元削減・データ圧縮等を目的として、データをより少数の主成分と呼ばれる合成変数で記述します。高次元の特徴空間におけるデータのばらつき（分散）をより低次元な空間において出来るだけ再現するため、「主成分の分散最大化」と「主成分の直行化」が必要必要な手続きになります。これは固有値と固有ベクトルを求めることにより実現されます。
実装
まずは必要なライブラリをインポートします。
In
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.datasets import load_breast_cancer
from sklearn.decomposition import PCA # 結果比較用
データの読込み
データセットは何でも良いですが、そこそこ特徴次元数のある breast_cancer を使うことにします。(特徴次元数 : 30)
Bunch型のデータから使用する部分を抽出します。
In
breast_cancer = load_breast_cancer()
data = pd.DataFrame(breast_cancer.data, columns=breast_cancer.feature_names)
data.head()
正規化
単位が不揃いで値に大きなばらつきがあるため、正規化を行います。
In
data_stats = data.describe().transpose()
def norm(x):
return (x-data_stats['mean'])/data_stats['std']
data = norm(data)
手計算用にNumPy配列に変換します。
In
data_array = data.values
print(data_array)
print(data_array.shape)
Out
[[ 1.09609953 -2.0715123 1.26881726 ... 2.2940576 2.74820411 1.93531174]
[ 1.82821197 -0.35332152 1.68447255 ... 1.08612862 -0.24367526 0.28094279]
[ 1.5784992 0.45578591 1.56512598 ... 1.95328166 1.15124203 0.20121416]
...
[ 0.70166686 2.04377549 0.67208442 ... 0.41370467 -1.10357792 -0.31812924]
[ 1.83672491 2.33440316 1.98078127 ... 2.28797231 1.9173959 2.21768395]
[-1.80681144 1.22071793 -1.81279344 ... -1.7435287 -0.04809589 -0.75054629]]
(569, 30)
分散共分散行列の計算
変数$j$の分散$s_{j, j}$, および変数$j$と変数$k$の共分散$s_{j, k} (j \neq k)$はそれぞれ次のように表されます。
$$s_{j, j}=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i, j}-\bar{x}_{.j})^2$$
$$s_{j, k}=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i, j}-\bar{x_{.j}})(x_{i, k}-\bar{x_{.k}})$$
したがって、分散共分散行列$S$は次のように求められます。
In
n_instances = data_array.shape[0]
n_features = data_array.shape[1]
data_mean = np.mean(data_array, axis=0)
S = np.zeros((n_features, n_features), dtype=float)
for j in range(n_features):
for k in range(n_features):
for i in range(n_instances):
S[j][k] += (data_array[i][j] - data_mean[j])*(data_array[i][k] - data_mean[k]) / n_instances
固有値・固有ベクトルの計算
続いて$S$ の固有値・固有ベクトルを計算します。
固有値・固有ベクトルの計算には numpy.linalg.eig を利用します。
In
# 固有値
eigenvalues = np.linalg.eig(S)[0]
# 固有ベクトル
U = np.linalg.eig(S)[1]
print(eigenvalues)
Out
array([1.32582657e+01, 5.68135223e+00, 2.81299652e+00, 1.97715956e+00, 1.64583295e+00, 1.20523472e+00, 6.74033435e-01, 4.75779500e-01, 4.16162133e-01, 3.50077124e-01, 2.93399148e-01, 2.60702387e-01, 2.40933318e-01, 1.56733784e-01, 9.39695257e-02, 7.97224445e-02, 5.92946458e-02, 5.25263076e-02, 4.93906364e-02, 1.32811001e-04, 7.47487099e-04, 1.58654466e-03, 6.88833652e-03, 8.16326791e-03, 1.54540635e-02, 1.80232759e-02, 2.42980595e-02, 2.73911786e-02, 3.11046408e-02, 2.99202175e-02])
これで元データを混合変数へと変換する行列$U$を求めることができました。
ここで得られた$U$は、次で表されるようにように$S$を対角化する行列です。
$$U^TSU=diag(λ1, ..., λ_p)$$
In
U_inv = np.linalg.inv(U)
diag = np.dot(np.dot(U_inv, S), U)
diag
Out
array([[ 1.32582657e+01, 1.27416995e-15, 2.59655658e-16, -1.25261426e-15, 4.04168592e-16, -6.09310220e-16, -1.74696331e-15, 2.38945931e-15, -4.38301407e-16, -1.20268061e-15, 3.54961084e-16, -6.99035253e-16, 3.30506578e-16, -2.56258925e-16, 2.75352231e-15, 2.24139149e-16, -1.33023827e-15,
-6.37619774e-16, 3.22970991e-15, 4.21239268e-15, 3.08658922e-15, -2.76021803e-15, 2.97518279e-15, 9.45346327e-17, 1.06484694e-15, -7.14021028e-16, 2.37497868e-15, -1.50135500e-15, -5.45871698e-16, -1.40709743e-15],
[ 8.80527252e-16, 5.68135223e+00, -2.31827482e-16, 1.37268856e-16, -9.73260548e-16, -4.09033288e-16, -4.35578462e-16, 8.43022527e-17, -2.69433268e-16, -9.23904427e-16, 1.19330641e-15, -1.07633683e-15, 7.72295909e-16, 2.10526158e-16, -4.86324918e-16, 1.17279288e-15, 1.39606877e-16,
8.77355419e-16, 1.43290742e-17, -6.60369435e-16, -1.41242614e-15, 1.84603892e-15, 2.74993617e-16, -2.57557446e-16, 4.23073643e-16, -1.21349087e-16, 3.49565329e-16, 1.95562317e-17, 2.80092695e-16, 1.42138866e-16],
...
[ 2.11188838e-15, 5.39314735e-16, 1.63313324e-16, -6.34796943e-16, 3.03143545e-16, -1.45993795e-16, -1.09593888e-16, 1.35971780e-16, -1.47104768e-17, 5.63829035e-17, -2.31515095e-16, 4.54289262e-17, 4.46706707e-17, 1.42878079e-17, -6.69645307e-17, 2.12551787e-17, -9.76503610e-17,
-9.83490621e-17, 1.30410900e-16, -3.26122703e-17, -4.39348644e-17, 3.27142835e-17, 6.77437666e-18, 1.41538219e-17, -8.21726757e-17, 1.75046448e-17, -6.75116428e-17, -1.40021512e-17, 3.11046408e-02, -3.80548301e-17],
[ 1.81691436e-15, -4.19837541e-16, 1.69790524e-16, -1.32978080e-15, -2.92109118e-17, -9.51987993e-17, -3.02233535e-16, 1.12616692e-16, -7.39060669e-17, 5.25773657e-17, -3.85867511e-17, 1.17268833e-16, -9.59743723e-18, 1.11617116e-17, -1.12287607e-16, -8.42366298e-18, 9.72260224e-17,
1.10219780e-16, 1.01366564e-17, 2.55796481e-16, -2.68489853e-17, 1.20429197e-16, 1.10363419e-16, -1.23429551e-17, 6.66921317e-17, 3.95025490e-17, -3.87011501e-17, -8.63280591e-17, 2.40399949e-17, 2.99202175e-02]])
対角成分以外は0、対角成分は先に求めた固有値に一致していることが確認できます。
寄与率の計算
第$j$成分の寄与率 (contribution ratio) は次の式で求められます。各固有値の総和に対する割合と言えます。
$$c_j = \frac{λj}{λ1+ ... +λ_p}$$
In
contribution_ratios = []
sum_eigen = np.sum(eigenvalues)
for i in range(eigenvalues.shape[0]):
contribution_ratio = eigenvalues[i] / sum_eigen
contribution_ratios.append(contribution_ratio)
print(contribution_ratios)
Out
[0.44272025607526355, 0.1897118204403308, 0.09393163257431368, 0.06602134915470165, 0.054957684923462646,
0.040245220398833485, 0.02250733712982507, 0.015887238000213286, 0.013896493745591104, 0.0116897818941315,
0.009797189875980144, 0.008705379007378822, 0.008045249871967327, 0.005233657454926348, 0.0031378321676274047,
0.0026620933651523055, 0.0019799679253242703, 0.001753959450226365, 0.001649253059225152, 4.4348274273267065e-06,
2.4960103246860732e-05, 5.297792903810818e-05, 0.00023001546250597875, 0.00027258799547752054, 0.0005160423791651501,
0.0006018335666716387, 0.000811361258899119, 0.0009146467510543492, 0.0010386467483387103, 0.0009990964637002163]
累積寄与率とともに寄与率を可視化してみます。
In
plt.figure()
left = range(1, len(contribution_ratios)+1)
plt.bar(left, contribution_ratios, width=1.0)
plt.plot(left, cumulated_contribution_ratios, c='black')
plt.ylim([0, 1])
plt.show()
scikit-learnパッケージとの比較
scikit-learnパッケージを使えば上記の計算が2行です。
pca = PCA(n_components=30)
pca.fit(data)
fit後の explained_variance_ratio_ 属性が寄与率に相当し、手計算の結果と一致していることが分かります。
In
pca.explained_variance_ratio_
Out
array([4.42720256e-01, 1.89711820e-01, 9.39316326e-02, 6.60213492e-02, 5.49576849e-02, 4.02452204e-02, 2.25073371e-02, 1.58872380e-02, 1.38964937e-02, 1.16897819e-02, 9.79718988e-03, 8.70537901e-03, 8.04524987e-03, 5.23365745e-03, 3.13783217e-03, 2.66209337e-03, 1.97996793e-03, 1.75395945e-03,
1.64925306e-03, 1.03864675e-03, 9.99096464e-04, 9.14646751e-04, 8.11361259e-04, 6.01833567e-04, 5.16042379e-04, 2.72587995e-04, 2.30015463e-04, 5.29779290e-05, 2.49601032e-05, 4.43482743e-06])
固有ベクトルはというと、components_属性に逆行列として存在するようです。
こちらも手計算で求めたU_invに一致します。
In
pca.components_
Out
array([[ 2.18902444e-01, 1.03724578e-01, 2.27537293e-01, 2.20994985e-01, 1.42589694e-01, 2.39285354e-01, 2.58400481e-01, 2.60853758e-01, 1.38166959e-01, 6.43633464e-02, 2.05978776e-01, 1.74280281e-02, 2.11325916e-01, 2.02869635e-01, 1.45314521e-02, 1.70393451e-01, 1.53589790e-01,
1.83417397e-01, 4.24984216e-02, 1.02568322e-01, 2.27996634e-01, 1.04469325e-01, 2.36639681e-01, 2.24870533e-01, 1.27952561e-01, 2.10095880e-01, 2.28767533e-01, 2.50885971e-01, 1.22904556e-01, 1.31783943e-01],
[-2.33857132e-01, -5.97060883e-02, -2.15181361e-01, -2.31076711e-01, 1.86113023e-01, 1.51891610e-01, 6.01653628e-02, -3.47675005e-02, 1.90348770e-01, 3.66575471e-01, -1.05552152e-01, 8.99796818e-02, -8.94572342e-02, -1.52292628e-01, 2.04430453e-01, 2.32715896e-01, 1.97207283e-01,
1.30321560e-01, 1.83848000e-01, 2.80092027e-01, -2.19866379e-01, -4.54672983e-02, -1.99878428e-01, -2.19351858e-01, 1.72304352e-01, 1.43593173e-01, 9.79641143e-02, -8.25723507e-03, 1.41883349e-01, 2.75339469e-01],
...
[ 2.11460455e-01, -1.05339342e-02, 3.83826098e-01, -4.22794920e-01, -3.43466700e-03, -4.10167739e-02, -1.00147876e-02, -4.20694931e-03, -7.56986244e-03, 7.30143287e-03, 1.18442112e-01, -8.77627920e-03, -6.10021933e-03, -8.59259138e-02, 1.77638619e-03, 3.15813441e-03, 1.60785207e-02,
-2.39377870e-02, -5.22329189e-03, -8.34191154e-03, -6.35724917e-01, 1.72354925e-02, 2.29218029e-02, 4.44935933e-01, 7.38549171e-03, 3.56690392e-06, -1.26757226e-02, 3.52404543e-02, 1.34042283e-02, 1.14776603e-02],
[-7.02414091e-01, -2.73661018e-04, 6.89896968e-01, 3.29473482e-02, 4.84745766e-03, -4.46741863e-02, -2.51386661e-02, 1.07726530e-03, 1.28037941e-03, 4.75568480e-03, 8.71109373e-03, 1.07103919e-03, -1.37293906e-02, -1.10532603e-03, 1.60821086e-03, -1.91562235e-03, 8.92652653e-03,
2.16019727e-03, -3.29389752e-04, -1.79895682e-03, 1.35643056e-01, -1.02053601e-03, -7.97438536e-02, -3.97422838e-02, -4.58327731e-03, 1.28415624e-02, -4.02139168e-04, 2.28844179e-03, -3.95443454e-04, -1.89429245e-03]])