[Python] Vergleich der Theorie und Implementierung der Hauptkomponentenanalyse durch Python (PCA, Kernel PCA, 2DPCA)

Wir werden die Theorie der Hauptkomponentenanalyse (PCA), die als Methode zur Dimensionsreduktion und Merkmalsextraktion bekannt ist, mit ihrer Anwendung (Kernel PCA, 2DPCA [^ 1]) und ihrer Implementierung in Python vergleichen.

Theorie der Hauptkomponentenanalyse

1. Hauptkomponente

$ M $ data $ a_1, \ ldots, a_M \ in \ mathbb {R} ^ n $, die $ 1 $ Hauptkomponente $ x_1 \ in \ mathbb {R} ^ n $ sind die zentralisierten Daten $a_i - \bar{a},\quad \bar{a} = \frac1M\sum_{i=1}^M a_i$ Es wird als Vektor der Länge $ 1 $ in der Richtung berechnet, die die Streuung von maximiert. Mit anderen Worten kann es als die optimale Lösung für das folgende Maximierungsproblem interpretiert werden.

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

Wenn hier die verteilte, gemeinsam verteilte Matrix der Daten $ C ist: = \ 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}

Daher ist die erste $ 1 $ Hauptkomponente $ x_1 $

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

Und stimmt mit dem Eigenvektor überein, der dem maximalen Eigenwert der Kovarianzmatrix $ C $ entspricht.

Zweite Hauptkomponente

Nach dem Subtrahieren der Richtung der ersten Hauptkomponente von den zentralisierten Daten $ a_i- \ bar {a} $ wird dieselbe Berechnung wie für die erste Hauptkomponente durchgeführt. Mit anderen Worten

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

Wird behandelt. Hier,

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

Die Formel kann transformiert werden. Sei $ P_1 = I_n-x_1x_1 ^ \ top $, der Eigenwert der verteilten gemeinsam verteilten Matrix $ C $ sei $ \ lambda_1 \ geq \ cdots \ geq \ lambda_n $, und die entsprechenden Eigenvektoren seien $ v_1, \ ldots, v_n $. Aus der Orthogonalität der Eigenvektoren einer reellen symmetrischen Matrix

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

Hält. Ebenfalls,

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

Daher ist die zweite $ 2 $ Hauptkomponente $ 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\}

Und stimmt mit dem Eigenvektor überein, der dem zweitgrößten Eigenwert $ \ lambda_2 $ in der Kovarianzmatrix $ C $ entspricht.

Aus einer ähnlichen Diskussion ist ersichtlich, dass die $ k $ Hauptkomponente $ x_k $ mit dem Eigenvektor übereinstimmt, der $ \ lambda_k $ entspricht, und die Hauptkomponentenanalyse der Daten $ a_1, \ ldots, a_M $ eine verteilte, gemeinsam verteilte Matrix $ C ist. Es wird auf das Eigenwertproblem von $ reduziert.

Dimensionsreduzierung

Wenn $ d $ Hauptkomponenten $ x_1, \ ldots, x_d $ aus den Daten $ a_1, \ ldots, a_M $ erhalten werden, wird die Dimensionsreduktion der neuen Daten $ a \ in \ mathbb {R} ^ n $ durchgeführt. Überlegen. In Anbetracht der Tatsache, dass $ a $ durch eine lineare Kombination von $ d $ Hauptkomponenten angenähert wird, wird $ X_d = [x_1, \ ldots, x_d] $ verwendet.

a \approx X_d y + \bar{a}

Anruf. Daher ergibt sich aus der Orthogonalität der Hauptkomponenten $ \ left (X_d ^ \ top X_d = I_d \ right) $

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

Und $ d $ Dimensionen können reduziert werden.

Kernel-PCA-Theorie [^ 2]

Hauptbestandteil

$ \ Phi $ sei die Zuordnung vom Datenraum $ \ mathbb {R} ^ n $ zum Hilbert-Raum des reproduzierenden Kernels (RKHS) $ \ mathcal {H} $. des Weiteren,

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

Dann ist die Zielfunktion von 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}

Es wird. Da die Hauptkomponentenanalyse zu einem eindeutigen Wertproblem von $ \ frac1M \ tilde \ Phi \ tilde \ Phi ^ \ top $ führt,

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

Denk an. Wenn $ z = \ tilde \ Phi ^ \ top x $,

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

Anruf. Einsetzen in die linke bzw. rechte Seite der Eigengleichung

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

Und

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

Anruf. Wobei die Matrix $ \ tilde \ Phi ^ \ top \ tilde \ Phi \ in \ mathbb {R} ^ {M \ times M} $ $ k (a_i, a_j) = \ phi (a_i) ^ \ top \ ist Als 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}

Daher sei die Matrix mit allen Komponenten $ 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}

Es wird berechnet von. Aus dem Obigen ist der Eigenvektor $ x $, der dem Eigenwert $ \ lambda $ von $ \ frac1M \ tilde \ Phi \ tilde \ Phi ^ \ top $ entspricht, die $ M $ dimensionale quadratische Matrix $ \ frac1M \ tilde \ Phi ^ \ top \ Verwenden des Eigenvektors $ z \ in \ mathbb {R} ^ M $ entsprechend dem Eigenwert $ \ lambda $ von tilde \ Phi $

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

Ist berechnet.

Dimensionsreduzierung

Wenn $ d $ Hauptkomponenten $ x_1, \ ldots, x_d $ aus den Daten $ a_1, \ ldots, a_M $ erhalten werden, wird die Dimensionsreduktion der neuen Daten $ a \ in \ mathbb {R} ^ n $ durchgeführt. Überlegen. In Anbetracht der Tatsache, dass $ \ phi (a) $ durch eine lineare Kombination von $ d $ Hauptkomponenten angenähert wird, wird $ X_d = [x_1, \ ldots, x_d] $ verwendet.

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

Anruf. Daher ergibt sich aus der Orthogonalität der Hauptkomponenten $ \ left (X_d ^ \ top X_d = I_d \ right) $

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

Und $ d $ Dimensionen können reduziert werden.

(Ergänzung)

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

Bei der Orthogonalisierung mit

\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 $ ist standardisiert, um zu erfüllen.

2 DPCA-Theorie [^ 1]

Angesichts der bildartigen zweidimensionalen Daten $ A_1, \ ldots, A_M \ in \ mathbb {R} ^ {m \ times n} $, PCA auf $ \ mathbb {R} ^ {mn} $ Während es in einen eindimensionalen Vektor konvertiert und behandelt wird, behandelt es die zweidimensionale PCA (2DPCA) wie in zwei Dimensionen.

Hauptbestandteil

Ziehen Sie in Betracht, die stochastische Variable $ A \ in \ mathbb {R} ^ {m \ times n} $ durch den $ n $ -Dimensionsvektor $ x $ der $ m $ -Dimension zuzuordnen. Mit anderen Worten

y = A x \tag{1}

Erwägen Sie, den Merkmalsvektor $ y \ in \ mathbb {R} ^ m $ von zu erhalten. Zu diesem Zeitpunkt ist die verteilte mitverteilte Matrix von $ y $ ein Index von $ x $

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

Verfolgen Sie $ \ mathrm {tr} , S_x $, dh maximieren Sie die Verteilung von $ 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}

Wenn daher die Daten $ A_1, \ ldots, A_M \ in \ mathbb {R} ^ {m \ times n} $ beobachtet werden, ist $ \ 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)

Verwenden von

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

In diesem Moment,

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

Und $ G_t $ ergibt das maximale Eigenwertproblem.

Dimensionsreduzierung

Wenn $ d $ Eigenvektoren $ x_1, \ ldots, x_d $ von $ G_t $ aus den Daten $ A_1, \ ldots, A_M $, neuen Daten $ A \ in \ mathbb {R} ^ {m \ erhalten werden Betrachten Sie die Dimensionsreduktion der Zeiten n} $. Wenn $ X_d = [x_1, \ ldots, x_d] \ in \ mathbb {R} ^ {n \ times d} $, dann aus Gleichung (1)

Y_d = AX_d

Und $ m \ times d $ kann auf Dimensionen reduziert werden. Ebenfalls,

A \approx Y_dX_d^\top

Ist ungefähr restauriert.

Vergleichstabelle von PCA und 2DPCA

	PCA	2DPCA
Eingabedimension	\mathbb{R}^{mn}	\mathbb{R}^{m\times n}
Komprimierte Abmessung	\mathbb{R}^{d}	\mathbb{R}^{m\times d}
Dimension der verteilten co-verteilten Matrix	\mathbb{R}^{mn\times mn}	\mathbb{R}^{n\times n}
X_dAbmessungen	\mathbb{R}^{mn\times d}	\mathbb{R}^{n\times d}

Beispiel für die Ausführung in Python

PCA und [Kernel PCA](https://scikit-learn.org/stable/modules/generated/ sklearn.decomposition.KernelPCA.html) verwendete die Implementierung von scikit-learn und 2DPCA verwendete die Implementierung selbst. Der implementierte Code wird am Ende beschrieben.

Datensatzvorbereitung

Als Ausführungsbeispiel haben wir den handgeschriebenen numerischen Bilddatensatz MNIST 28x28 verwendet. Wie unten gezeigt, wurden 6902 Bilder von 0 zum Lernen verwendet, und ein Bild von 0 und ein Bild von 1 wurden zur Verifizierung verwendet.

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 lernen

Aus der Varianz-Co-Verteilungsmatrix wurden fünf Eigenvektoren gelernt. Beachten Sie, dass "KernelPCA" nicht aus der Dimensionsreduktion wiederhergestellt werden kann, es sei denn "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))

Die Lernzeit ist wie folgt in der Reihenfolge PCA, Kernel PCA, 2DPCA. Es stellt sich heraus, dass Kernel PCA besonders schwer ist.

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

Dimensionsreduktion und Bildwiederherstellung durch PCA

Die Dimensionsreduktion wurde mit "transform" durchgeführt und 0 Bilder wurden mit "inverse_transform" wiederhergestellt.

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

Die Berechnungszeit ist wie folgt in der Reihenfolge 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

Restaurierungsergebnis

Die Dimension von Bild 1 wurde auf die gleiche Weise reduziert und wiederhergestellt.

2DPCA-Quellcode

Das ausgeführte Notizbuch befindet sich auf 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)

[^ 2]: Kenji Fukumizu, "Einführung in die Kernel-Methodendatenanalyse mit einem positiven Kernel mit festem Wert", Asakura Shoten, 2010.