【Python】主成分分析の理論と Python による実装の比較 (PCA, Kernel PCA, 2DPCA)

次元削減や特徴量抽出の手法として有名な主成分分析 (Principal Component Analysis, PCA) とその応用 (Kernel PCA, 2DPCA[^1]) の理論と Python での実装の比較を行なう．

## 主成分分析の理論
### 第 1 主成分
$M$ 個のデータ $a_1,\ldots,a_M\in\mathbb{R}^n$ の第 $1$ 主成分 $x_1\in\mathbb{R}^n$ は，中心化されたデータ 
$$a_i - \bar{a},\quad \bar{a} = \frac1M\sum_{i=1}^M a_i$$ 
の分散が最大となる方向で長さ $1$ のベクトルとして求められる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/2e320998-5b92-0bb2-4b07-fd073b66be03.png)
つまり，次の最大化問題の最適解と解釈できる．

```math
\begin{align}
\max_{x\in\mathbb{R}^n} &\quad \frac1M\sum_{i=1}^M\left((a_i - \bar{a})^\top x\right)^2\\
\text{s.t.} &\quad \|x\|_2 = 1
\end{align}
```

ここで，データの分散共分散行列を $C := \frac1M\sum_{i=1}^M(a_i - \bar{a})(a_i - \bar{a})^\top$ とすると，

```math
\begin{align}
\max_{x\in\mathbb{R}^n} \ \frac1M\sum_{i=1}^M\left((a_i - \bar{a})^\top x\right)^2 &= \max_{x\in\mathbb{R}^n} \ x^\top\left[\frac1M\sum_{i=1}^M(a_i - \bar{a})(a_i - \bar{a})^\top\right] x\\
&= \max_{x\in\mathbb{R}^n}\ x^\top C x
\end{align}
```

より，第 $1$ 主成分 $x_1$ は

```math
x_1 = \mathop{\mathrm{arg~max}}_{x\in\mathbb{R^n}}\left\{x^\top C x \mid \|x\|_2 = 1\right\}
```

となり，分散共分散行列 $C$ の最大固有値に対応する固有ベクトルと一致する．

### 第 2 主成分
中心化されたデータ $a_i - \bar{a}$ から第 1 主成分方向を引いた上で，第 1 主成分と同様の計算を行なう．つまり，

```math
\begin{align}
\max_{x\in\mathbb{R}^n} &\quad \frac1M\sum_{i=1}^M\left(\left(a_i - \bar{a} - \left(x_1^\top (a_i - \bar{a})\right)x_1\right)^\top x\right)^2\\
\text{s.t.} &\quad \|x\|_2 = 1
\end{align}
```

を扱う．ここで，

```math
\begin{align}
a_i - \bar{a} - \left(x_1^\top (a_i - \bar{a})\right)x_1 &= a_i - \bar{a} - \left(x_1x_1^\top\right)(a_i - \bar{a})\\
&= \left(I_n - x_1x_1^\top\right)(a_i - \bar{a})
\end{align}
```

と式変形できる．$P_1 = I_n - x_1x_1^\top$ とし，分散共分散行列 $C$ の固有値を $\lambda_1\geq\cdots\geq\lambda_n$，対応する固有ベクトルを $v_1, \ldots, v_n$ とすると，実対称行列の固有ベクトルの直交性より，

```math
\begin{cases}
P_1v_i = 0 & (i = 1)\\
P_1v_i = v_i & (i \neq 1)
\end{cases}
```

が成り立つ．また，

```math
\begin{align}
\max_{x\in\mathbb{R}^n} \ \frac1M\sum_{i=1}^M\left(\left(a_i - \bar{a} - \left(x_1^\top (a_i - \bar{a})\right)x_1\right)^\top x\right)^2 &= \max_{x\in\mathbb{R}^n} \ \frac1M\sum_{i=1}^M\left((a_i - \bar{a})^\top P_1 x\right)^2\\
&= \max_{x\in\mathbb{R}^n}\ (P_1x)^\top C (P_1x)
\end{align}
```

より，第 $2$ 主成分 $x_2$ は

```math
x_2 = \mathop{\mathrm{arg~max}}_{x\in\mathbb{R^n}}\left\{(P_1x)^\top C (P_1x) \mid \|x\|_2 = 1\right\}
```

となり，分散共分散行列 $C$ の 2 番目に大きい固有値 $\lambda_2$ に対応する固有ベクトルと一致する．

同様の議論より，第 $k$ 主成分 $x_k$ は $\lambda_k$ に対応する固有ベクトルと一致することがわかり，データ $a_1, \ldots, a_M$ の主成分分析は分散共分散行列 $C$ の固有値問題に帰着される．

### 次元削減
データ $a_1, \ldots, a_M$ から $d$ 個の主成分 $x_1, \ldots, x_d$ が得られたとき，新たなデータ $a\in\mathbb{R}^n$ の次元削減を考える．$a$ を$d$ 個の主成分の線形結合で近似することを考えると，$X_d = [x_1,\ldots,x_d]$ を用いて，

```math
a \approx X_d y + \bar{a}
``` 

とかける．よって，主成分の直交性 $\left(X_d^\top X_d = I_d\right)$ より，

```math
y = X_d^\top (a - \bar{a})
```

と $d$ 次元に次元削減することができる．

## Kernel PCA の理論[^2]
### 主成分
データ空間 $\mathbb{R}^n$ からある再生核ヒルベルト空間 (Reproducing Kernel Hilbert Space, RKHS) $\mathcal{H}$ への写像を $\phi$ とする．さらに，

```math
\bar{\phi} = \frac1M \sum_{i=1}^M \phi(a_i)\\
\Phi = [\phi(a_1),\ldots, \phi(a_M)]\\
\tilde{\Phi} = \left[\phi(a_1) - \bar{\phi},\ldots, \phi(a_M) - \bar{\phi}\right]
```

とすると，PCA の目的関数は

```math
\begin{align}
\frac1M\sum_{i=1}^M\left(\left(\phi(a_i) - \bar{\phi}\right)^\top x\right)^2 &= \frac1M x^\top\left[\sum_{i=1}^M\left(\phi(a_i) - \bar{\phi}\right)\left(\phi(a_i) - \bar{\phi}\right)^\top\right] x\\
&= \frac1M x^\top\tilde\Phi\tilde\Phi^\top x
\end{align}
```

となる．主成分分析は $\frac1M\tilde\Phi\tilde\Phi^\top$ の固有値問題に帰着されるので，

```math
\frac1M\tilde\Phi\tilde\Phi^\top x = \lambda x
```

を考える．$z = \tilde\Phi^\top x$ とすると，

```math
x = \frac1{M\lambda}\tilde\Phi z
```

とかける．これを固有方程式の左辺，右辺にそれぞれ代入すると，

```math
\begin{align}
\frac1M\tilde\Phi\tilde\Phi^\top x &= \frac1M\tilde\Phi\tilde\Phi^\top\left(\frac1{M\lambda}\tilde\Phi z\right)\\
&= \frac1{M^2\lambda}\tilde\Phi\left(\tilde\Phi^\top\tilde\Phi\right)z\\
\lambda x &= \frac1M\tilde\Phi z
\end{align}
```

となり，

```math
\begin{align}
\frac1{M^2\lambda}\tilde\Phi\left(\tilde\Phi^\top\tilde\Phi\right)z = \frac1M\tilde\Phi z \iff \tilde\Phi\left(\frac1M\tilde\Phi^\top\tilde\Phi\right)z = \tilde\Phi \lambda z
\end{align}
```

とかける．ここで，行列 $\tilde\Phi^\top\tilde\Phi\in\mathbb{R}^{M\times M}$ は，$k(a_i, a_j) = \phi(a_i)^\top\phi(a_j)$ として，

```math
\begin{align}
\left(\tilde\Phi^\top\tilde\Phi\right)_{i,j} &= (\phi(a_i) - \bar{\phi})^\top(\phi(a_j) - \bar{\phi})\\
&= k(a_i, a_j) - \frac1M\sum_{\ell=1}^Mk(a_i, a_\ell) - \frac1M\sum_{k=1}^Mk(a_k, a_j) + \frac1{M^2}\sum_{k, \ell=1}^Mk(a_k, a_\ell)
\end{align}
```

より，全成分が $1$ の行列を $\mathbf{1}\in\mathbb{R}^{M\times M}$，$\Phi^\top\Phi = K$ として，

```math
\tilde\Phi^\top\tilde\Phi = K - \frac1M K\mathbf{1} - \frac1M\mathbf{1} K+ \frac1{M^2}\mathbf{1} K \mathbf{1}
```

で計算される．以上より，$\frac1M\tilde\Phi\tilde\Phi^\top$ の固有値 $\lambda$ に対応する固有ベクトル $x$ は，$M$ 次元正方行列 $\frac1M\tilde\Phi^\top\tilde\Phi$ の固有値 $\lambda$ に対応する固有ベクトル $z\in\mathbb{R}^M$ を用いて，

```math
x = \frac1{M\lambda}\tilde\Phi z
```

と計算される．

### 次元削減
データ $a_1, \ldots, a_M$ から $d$ 個の主成分 $x_1, \ldots, x_d$ が得られたとき，新たなデータ $a\in\mathbb{R}^n$ の次元削減を考える．$\phi(a)$ を$d$ 個の主成分の線形結合で近似することを考えると，$X_d = [x_1,\ldots,x_d]$ を用いて，

```math
\phi(a) \approx X_d y + \bar{\phi}
``` 

とかける．よって，主成分の直交性 $\left(X_d^\top X_d = I_d\right)$ より，

```math
y = X_d^\top \left(\phi(a) - \bar{\phi}\right)
```

と $d$ 次元に次元削減することができる．

(補足)

```math
\left<x_i, x_j\right>_{\mathcal H} = \delta_{i, j}
```

と直交化する際，

```math
\left<x_i, x_i\right>_{\mathcal H} = \frac1{M\lambda_i^2}z_i^\top \left(\frac1M\Phi^\top\Phi \right)z_i = \frac{1}{M\lambda_i}\|z_i\|_2^2 = 1
```

をみたすように $z_i$ は規格化される．

## 2DPCA の理論[^1]

画像のような 2 次元データ $A_1, \ldots, A_M\in\mathbb{R}^{m\times n}$ が与えられたとき，PCA では$\mathbb{R}^{mn}$ 上の 1 次元ベクトルに変換して扱うのに対し，Two-Dimensional PCA (2DPCA) では 2 次元のまま扱う．

### 主成分

確率変数 $A\in\mathbb{R}^{m\times n}$ を $n$ 次元ベクトル $x$ により $m$ 次元に写像することを考える．つまり，

```math
y = A x \tag{1}
```

により特徴量ベクトル $y\in\mathbb{R}^m$ を得ることを考える．このとき，$x$ の指標として $y$ の分散共分散行列

```math
\begin{align}
S_x &= \mathbb{E}\left[\left(y - \mathbb{E}[y]\right)\left(y - \mathbb{E}[y]\right)^\top\right]\\
&= \mathbb{E}\left[\left(Ax - \mathbb{E}[Ax]\right)\left(Ax - \mathbb{E}[Ax]\right)^\top\right]\\
&= \mathbb{E}\left[\left(A - \mathbb{E}[A]\right)x(\left(A - \mathbb{E}[A]\right)x)^\top\right]\\
\end{align}
```

のトレース $\mathrm{tr}\,S_x$，すなわち $y$ の分散を最大化する．

```math
\begin{align}
\mathrm{tr}\,S_x &= \mathrm{tr}\,\mathbb{E}\left[\left(A - \mathbb{E}[A]\right)x(\left(A - \mathbb{E}[A]\right)x)^\top\right]\\
&= \mathrm{tr}\,\mathbb{E}\left[x^\top\left(A - \mathbb{E}[A]\right)^\top\left(A - \mathbb{E}[A]\right)x\right]\\
&= x^\top\mathbb{E}\left[\left(A - \mathbb{E}[A]\right)^\top\left(A - \mathbb{E}[A]\right)\right]x
\end{align}
```

より，データ $A_1, \ldots, A_M\in\mathbb{R}^{m\times n}$ が観測されたとき，$\mathrm{tr}\,S_x$ は，

```math
\bar{A} = \frac1M\sum_{i=1}^MA_i\\
G_t = \frac1M\sum_{i=1}^M\left(A_i - \bar{A}\right)^\top\left(A_i - \bar{A}\right)
```

を用いて

```math
\mathrm{tr}\,S_x = x^\top G_t x
```

と計算でき，このとき，

```math
x_1 = \mathop{\mathrm{arg~max}}_{x\in\mathbb{R}^n}\left\{x^\top G_t x \mid \|x\|_2 = 1\right\}
```

と $G_t$ の最大固有値問題に帰着する．

### 次元削減

データ $A_1, \ldots, A_M$ から $G_t$ の $d$ 個の固有ベクトル $x_1, \ldots, x_d$ が得られたとき，新たなデータ $A\in\mathbb{R}^{m\times n}$ の次元削減を考える．$X_d = [x_1, \ldots, x_d] \in\mathbb{R}^{n\times d}$ とすると (1) 式より

```math
Y_d = AX_d
```

と $m\times d$ 次元に次元削減することができる．また，

```math
A \approx Y_dX_d^\top
``` 

と近似的に復元される．

### PCA と 2DPCA の比較表

|| PCA | 2DPCA|
|:-:|:-:|:-:|
| 入力次元 | $\mathbb{R}^{mn}$ | $\mathbb{R}^{m\times n}$ |
| 圧縮次元 | $\mathbb{R}^{d}$ | $\mathbb{R}^{m\times d}$ |
| 分散共分散行列の次元 | $\mathbb{R}^{mn\times mn}$ | $\mathbb{R}^{n\times n}$ |
| $X_d$ の次元 | $\mathbb{R}^{mn\times d}$ | $\mathbb{R}^{n\times d}$ |

## Python での実行例

[PCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html) と [Kernel PCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.KernelPCA.html) は scikit-learn の実装を，2DPCA は自分で実装したものを使用した．実装したコードは最後に記載する．

### データセットの準備

実行例として，28x28 の手書き数字画像データセット MNIST を使用した．次のように 0 の画像 6902 枚を学習に，0 と 1 の画像を 1 枚ずつ検証に用いた．

```python
from sklearn import datasets
mnist = datasets.fetch_openml('mnist_784', data_home='../../data/')
mnist_data = mnist.data / 255.
mnist_0 = mnist_data[mnist.target == '0']

mnist_shape = (28, 28)
train = mnist_0[:-1]
test_0 = mnist_0[-1].reshape(1, -1)
test_1 = mnist_data[mnist.target == '1'][-1].reshape(1, -1)
```

### PCA の学習

分散共分散行列から 5 個の固有ベクトルを学習した．`KernelPCA` は `fit_inverse_transform=True` としないと次元削減から復元ができないことに注意する．

```python
from sklearn.decomposition import PCA, KernelPCA

n_components = 5
%time pca = PCA(n_components=n_components).fit(train)
%time kpca = KernelPCA(n_components=n_components, fit_inverse_transform=True).fit(train)
%time tdpca = TwoDimensionalPCA(n_components=n_components).fit(train.reshape(-1, *mnist_shape))
```

学習時間は PCA, Kernel PCA, 2DPCA の順に以下の通り．Kernel PCA が特に重いことがわかる．

```
CPU times: user 545 ms, sys: 93.9 ms, total: 639 ms
Wall time: 310 ms
CPU times: user 10.9 s, sys: 1.22 s, total: 12.1 s
Wall time: 7.8 s
CPU times: user 97.6 ms, sys: 59.3 ms, total: 157 ms
Wall time: 195 ms
```

### PCA による次元削減と画像の復元

`transform` で次元削減，`inverse_transform` で 0 の画像の復元を行なった．

```python
%time pca_result_0 = pca.inverse_transform(pca.transform(test_0))
%time kpca_result_0 = kpca.inverse_transform(kpca.transform(test_0))
%time tdpca_result_0 = tdpca.inverse_transform(tdpca.transform(test_0.reshape(1, *mnist_shape)))
```

計算時間は PCA, Kernel PCA, 2DPCA の順に以下の通り．

```
CPU times: user 1.17 ms, sys: 2.72 ms, total: 3.89 ms
Wall time: 4.26 ms
CPU times: user 14.3 ms, sys: 2.24 ms, total: 16.5 ms
Wall time: 11.5 ms
CPU times: user 90 µs, sys: 3 µs, total: 93 µs
Wall time: 103 µs
```

復元結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/735d6131-0f4e-d267-2211-2405da1bf141.png)


1 の画像に対しても同様に次元削減・復元を行った．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394475/54c34461-2c6d-ffb7-2951-3267e010aa11.png)


### 2DPCA のソースコード

実行したノートブックは [Github](https://github.com/Mrrmm510/ManifoldOptimization/blob/master/notebook/pca/pca-comparison.ipynb) に．

```python
from typing import Union

import numpy as np
from scipy.linalg import eigh
from sklearn.exceptions import NotFittedError


class TwoDimensionalPCA:
    def __init__(self, n_components: Union[int, None] = None):
        self.n_components_ = n_components

        self.components_ = None
        self.mean_ = None

    def fit(self, X: np.ndarray):
        """
        Parameters
        ----------
        X : np.ndarray, shape (n_samples, height, width)
        """
        if X.ndim != 3:
            raise ValueError(f"Expected 3D array, got {X.ndim}D array instead")

        self.mean_ = np.mean(X, axis=0)
        cov = np.mean([x.T @ x for x in X - self.mean_], axis=0)
        n_features = cov.shape[0]

        if self.n_components_ is None or self.n_components_ > n_features:
            self.n_components_ = n_features
        self.components_ = eigh(cov, eigvals=(n_features - self.n_components_, n_features - 1))[1][:, ::-1]
        return self

    def transform(self, X: np.ndarray) -> np.ndarray:
        """
        Parameters
        ----------
        X : np.ndarray, shape (n_samples, height, width)

        Returns
        -------
        : np.ndarray, shape (n_samples, height, n_components)
        """
        if self.components_ is None:
            raise NotFittedError(
                "This PCA instance is not fitted yet. "
                "Call 'fit' with appropriate arguments before using this estimator.")

        if X.ndim != 3:
            raise ValueError(f"Expected 3D array, got {X.ndim}D array instead")

        return np.array([x @ self.components_ for x in X])

    def inverse_transform(self, X: np.ndarray) -> np.ndarray:
        """
        Parameters
        ----------
        X : np.ndarray, shape (n_samples, height, n_components)

        Returns
        -------
        : np.ndarray, shape (n_samples, height, width)
        """
        return np.array([x @ self.components_.T for x in X])

    def fit_transform(self, X: np.ndarray) -> np.ndarray:
        """
        Parameters
        ----------
        X : np.ndarray, shape (n_samples, height, width)

        Returns
        -------
        : np.ndarray, shape (n_samples, height, n_components)
        """
        self.fit(X)
        return self.transform(X)

```

[^1]: [Yang, J., Zhang, D., Frangi, A. F., & Yang, J. Y. (2004). Two-dimensional PCA: a new approach to appearance-based face representation and recognition. IEEE transactions on pattern analysis and machine intelligence, 26(1), 131-137.](http://158.132.160.118/bitstream/10397/190/1/137.pdf)

[^2]: 福水健次，『カーネル法入門―正定値カーネルによるデータ解析―』，朝倉書店，2010.