Help us understand the problem. What is going on with this article?

PCAって数学的には何やってるんだっけ

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()

スクリーンショット 2020-11-30 15.38.37.png

正規化

単位が不揃いで値に大きなばらつきがあるため、正規化を行います。

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()

スクリーンショット 2020-11-30 16.40.23.png

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]])

参考文献

  1. 日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック
  2. 主成分分析の基礎 - inoccu.com
  3. sklearn.decomposition.PCA - scikit-learn.org
roki18d
某SIerでビッグデータ分析基盤のR&D業務に携わっています。Qiitaでは統計に関する記事をメインに書いていこうと思います。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away