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

# TL;DR

特徴次元数削減アルゴリズムの代表格といえば主成分分析 (PCA : Principal Component Analysis) ですね。ただ、やってることは感覚的に何となく分かっていても、線形代数的に何やってるのかをよく覚えておらず、気になったので復習してみた記事です。

PCAを手計算で実装し、結果をscikit-learnのPCAサブパッケージのものと比較します。（手計算といっても `sklearn.decomposition.PCA` を使わないでという意味であってNumPyは使います）

# 主成分分析の概要

高次元の特徴からなるデータの次元削減・データ圧縮等を目的として、データをより少数の主成分と呼ばれる合成変数で記述します。高次元の特徴空間におけるデータのばらつき（分散）をより低次元な空間において出来るだけ再現するため、「**主成分の分散最大化**」と「**主成分の直行化**」が必要必要な手続きになります。これは**固有値**と**固有ベクトル**を求めることにより実現されます。

# 実装

まずは必要なライブラリをインポートします。

```python: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型のデータから使用する部分を抽出します。

```python: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/fa1d3572-3365-9e47-3fd5-08edc62e2b66.png)

## 正規化

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

```python:In
data_stats = data.describe().transpose()

def norm(x):
    return (x-data_stats['mean'])/data_stats['std']

data = norm(data)
```

手計算用にNumPy配列に変換します。

```python:In
data_array = data.values
print(data_array)
print(data_array.shape)
```

```r: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$は次のように求められます。

```python: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` を利用します。

```python:In
# 固有値
eigenvalues = np.linalg.eig(S)[0]
# 固有ベクトル
U = np.linalg.eig(S)[1]

print(eigenvalues)
```

```r: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)$$

```python:In
U_inv = np.linalg.inv(U)
diag = np.dot(np.dot(U_inv, S), U)
diag
```

```r: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}$$

```python: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)
```

```r: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]
```

累積寄与率とともに寄与率を可視化してみます。

```python: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/a10ef00d-b0c9-5fff-df18-97a44586c82c.png)

## scikit-learnパッケージとの比較

scikit-learnパッケージを使えば上記の計算が2行です。

```
pca = PCA(n_components=30)
pca.fit(data)
```

fit後の `explained_variance_ratio_` 属性が寄与率に相当し、手計算の結果と一致していることが分かります。

```python:In
pca.explained_variance_ratio_
```
```r: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`に一致します。

```python:In
pca.components_
```
```r: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級対応 統計学実践ワークブック](https://www.amazon.co.jp/dp/478060852X/)
2. [主成分分析の基礎 - inoccu.com](https://inoccu.com/docs/datascience/pca.html)
3. [sklearn.decomposition.PCA - scikit-learn.org](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html)