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

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

次元削減や特徴量抽出の手法として有名な主成分分析 (Principal Component Analysis, PCA) とその応用 (Kernel PCA, 2DPCA1) の理論と 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
つまり,次の最大化問題の最適解と解釈できる.

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

主成分

データ空間 $\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 の理論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$ 次元に写像することを考える.つまり,

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 での実行例

PCAKernel 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 個の固有ベクトルを学習した.KernelPCAfit_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

復元結果
image.png

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

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)

Mrrmm252
ポケモン好きなおじさん 論文を読んで Python に
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
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  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
ユーザーは見つかりませんでした