次元削減や特徴量抽出の手法として有名な主成分分析 (Principal Component Analysis, PCA) とその応用 (Kernel PCA, 2DPCA¹) の理論と 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$ のベクトルとして求められる．

つまり，次の最大化問題の最適解と解釈できる．

\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$ とすると，

\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$ は

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 主成分と同様の計算を行なう．つまり，

\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}

を扱う．ここで，

\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$ とすると，実対称行列の固有ベクトルの直交性より，

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

が成り立つ．また，

\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$ は

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]$ を用いて，

a \approx X_d y + \bar{a}

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

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

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

Kernel PCA の理論²

主成分

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

\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 の目的関数は

\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$ の固有値問題に帰着されるので，

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

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

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

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

\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}

となり，

\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)$ として，

\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$ として，

\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$ を用いて，

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]$ を用いて，

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

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

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

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

(補足)

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

と直交化する際，

\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 の理論¹

画像のような 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$ 次元に写像することを考える．つまり，

y = A x \tag{1}

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

\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$ の分散を最大化する．

\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$ は，

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

を用いて

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

と計算でき，このとき，

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) 式より

Y_d = AX_d

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

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 と Kernel PCA は scikit-learn の実装を，2DPCA は自分で実装したものを使用した．実装したコードは最後に記載する．

データセットの準備

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

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 としないと次元削減から復元ができないことに注意する．

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 の画像の復元を行なった．

%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

復元結果

1 の画像に対しても同様に次元削減・復元を行った．

2DPCA のソースコード

実行したノートブックは Github に．

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)

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. ↩ ↩²
福水健次，『カーネル法入門―正定値カーネルによるデータ解析―』，朝倉書店，2010. ↩