[Python] Comparaison de la théorie de l'analyse des composants principaux et de l'implémentation par Python (PCA, Kernel PCA, 2DPCA)

Nous comparerons la théorie de l'Analyse en Composantes Principales (PCA), célèbre comme méthode de réduction de dimension et d'extraction de caractéristiques, et son application (Kernel PCA, 2DPCA [^ 1]) et son implémentation en Python.

Théorie de l'analyse en composantes principales

1er composant principal

$ M $ data $ a_1, \ ldots, a_M \ in \ mathbb {R} ^ n $, le premier composant principal $ 1 $ de $ x_1 \ in \ mathbb {R} ^ n $ est des données centralisées $a_i - \bar{a},\quad \bar{a} = \frac1M\sum_{i=1}^M a_i$ Il est calculé comme un vecteur de longueur $ 1 $ dans la direction qui maximise la dispersion de. En d'autres termes, il peut être interprété comme la solution optimale pour le problème de maximisation suivant.

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

Ici, si la matrice distribuée co-distribuée des données est $ 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}

Par conséquent, la première composante principale de 1 $ $ x_1 $ est

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

Et est d'accord avec le vecteur propre correspondant à la valeur propre maximale de la matrice de covariance $ C $.

Deuxième composant principal

Après avoir soustrait la direction du premier composant principal des données centralisées $ a_i- \ bar {a} $, le même calcul que pour le premier composant principal est effectué. En d'autres termes

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

Est traité. ici,

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

La formule peut être transformée. Soit $ P_1 = I_n-x_1x_1 ^ \ top $, la valeur propre de la matrice co-distribuée distribuée $ C $ soit $ \ lambda_1 \ geq \ cdots \ geq \ lambda_n $, et les vecteurs propres correspondants seront $ v_1, \ ldots, v_n $. A partir de l'orthogonalité des vecteurs propres d'une matrice symétrique réelle

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

Tient. Aussi,

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

Par conséquent, la deuxième composante principale $ 2 $ $ x_2 $ est

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

Et correspond au vecteur propre correspondant à la deuxième plus grande valeur propre $ \ lambda_2 $ dans la matrice de covariance $ C $.

A partir d'une discussion similaire, on peut voir que la composante principale $ k $ $ x_k $ correspond au vecteur propre correspondant à $ \ lambda_k $, et l'analyse en composante principale des données $ a_1, \ ldots, a_M $ est une matrice distribuée co-distribuée $ C. Il est réduit au problème des valeurs propres de $.

Réduction de dimension

Lorsque $ d $ composants principaux $ x_1, \ ldots, x_d $ sont obtenus à partir des données $ a_1, \ ldots, a_M $, la réduction de dimension des nouvelles données $ a \ in \ mathbb {R} ^ n $ est effectuée. Pense. Considérant que $ a $ est approximé par une combinaison linéaire de $ d $ composants principaux, $ X_d = [x_1, \ ldots, x_d] $ est utilisé.

a \approx X_d y + \bar{a}

Appel. Donc, à partir de l'orthogonalité des composantes principales $ \ left (X_d ^ \ top X_d = I_d \ right) $

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

Et les dimensions $ d $ peuvent être réduites.

Théorie PCA du noyau [^ 2]

Composant principal

Soit $ \ phi $ le mappage de l'espace de données $ \ mathbb {R} ^ n $ à l'espace de Hilbert du noyau reproducteur (RKHS) $ \ mathcal {H} $. plus loin,

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

Ensuite, la fonction objective de l'ACP est

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

Il devient. Puisque l'analyse en composantes principales aboutit à un problème de valeur unique de $ \ frac1M \ tilde \ Phi \ tilde \ Phi ^ \ top $,

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

penser à. Si $ z = \ tilde \ Phi ^ \ top x $,

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

Appel. En le substituant respectivement aux côtés gauche et droit de l'équation propre,

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

Et

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

Appel. Où la matrice $ \ tilde \ Phi ^ \ top \ tilde \ Phi \ in \ mathbb {R} ^ {M \ times M} $ est $ k (a_i, a_j) = \ phi (a_i) ^ \ top \ Comme 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}

Par conséquent, laissez la matrice avec tous les composants $ 1 $ être $ \ mathbf {1} \ in \ mathbb {R} ^ {M \ fois 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}

Il est calculé par. D'après ce qui précède, le vecteur propre $ x $ correspondant à la valeur propre $ \ lambda $ de $ \ frac1M \ tilde \ Phi \ tilde \ Phi ^ \ top $ est la matrice carrée dimensionnelle $ M $ $ \ frac1M \ tilde \ Phi ^ \ top \ En utilisant le vecteur propre $ z \ in \ mathbb {R} ^ M $ correspondant à la valeur propre $ \ lambda $ du tilde \ Phi $,

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

Est calculé.

Réduction de dimension

Lorsque $ d $ composants principaux $ x_1, \ ldots, x_d $ sont obtenus à partir des données $ a_1, \ ldots, a_M $, la réduction de dimension des nouvelles données $ a \ in \ mathbb {R} ^ n $ est effectuée. Pense. Considérant que $ \ phi (a) $ est approximé par une combinaison linéaire de $ d $ composants principaux, $ X_d = [x_1, \ ldots, x_d] $ est utilisé.

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

Appel. Donc, à partir de l'orthogonalité des composantes principales $ \ left (X_d ^ \ top X_d = I_d \ right) $

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

Et les dimensions $ d $ peuvent être réduites.

(Supplément)

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

Lors de l'orthogonalisation avec

\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 $ est standardisé pour répondre.

2 Théorie DPCA [^ 1]

Étant donné les données bidimensionnelles de type image $ A_1, \ ldots, A_M \ in \ mathbb {R} ^ {m \ times n} $, PCA sur $ \ mathbb {R} ^ {mn} $ Alors qu'il est converti en un vecteur unidimensionnel et manipulé, PCA bidimensionnel (2DPCA) le gère tel quel en deux dimensions.

Composant principal

Envisagez de mapper la variable stochastique $ A \ in \ mathbb {R} ^ {m \ times n} $ à la dimension $ m $ avec le vecteur de dimension $ n $ $ x $. En d'autres termes

y = A x \tag{1}

Pensez à obtenir le vecteur de caractéristiques $ y \ in \ mathbb {R} ^ m $ by. À ce stade, la matrice distribuée co-distribuée de $ y $ comme un indice de $ 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}

Tracez $ \ mathrm {tr} , S_x $, c'est-à-dire maximisez la distribution de $ 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}

Par conséquent, lorsque les données $ A_1, \ ldots, A_M \ in \ mathbb {R} ^ {m \ times n} $ sont observées, $ \ mathrm {tr} , S_x $ est

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

En utilisant

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

En ce moment,

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

Et $ G_t $ aboutit au problème de valeur propre maximale.

Réduction de dimension

Lorsque $ d $ vecteurs propres $ x_1, \ ldots, x_d $ de $ G_t $ sont obtenus à partir des données $ A_1, \ ldots, A_M $, nouvelles données $ A \ in \ mathbb {R} ^ {m \ Considérons la réduction de dimension des temps n} $. Si $ X_d = [x_1, \ ldots, x_d] \ in \ mathbb {R} ^ {n \ fois d} $, alors de l'équation (1)

Y_d = AX_d

Et $ m \ times d $ peuvent être réduits en dimensions. Aussi,

A \approx Y_dX_d^\top

Est approximativement restauré.

Tableau de comparaison PCA et 2DPCA

Exemple d'exécution en Python

PCA et [Kernel PCA](https://scikit-learn.org/stable/modules/generated/ sklearn.decomposition.KernelPCA.html) a utilisé l'implémentation de scikit-learn, et 2DPCA a utilisé l'implémentation par lui-même. Le code implémenté est décrit à la fin.

Préparation du jeu de données

Comme exemple d'exécution, nous avons utilisé l'ensemble de données d'images numériques manuscrites 28x28 MNIST. Comme indiqué ci-dessous, 6902 images de 0 ont été utilisées pour l'apprentissage, et une image de 0 et une image de 1 ont été utilisées pour la vérification.

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)

Apprentissage de l'APC

Cinq vecteurs propres ont été tirés de la matrice de variance-co-distribution. Notez que KernelPCA ne peut pas être restauré à partir de la réduction de dimension à moins que 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))

Le temps d'apprentissage est le suivant dans l'ordre PCA, Kernel PCA, 2DPCA. Il s'avère que Kernel PCA est particulièrement lourd.

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

Réduction dimensionnelle et restauration d'image par PCA

La réduction de dimension a été effectuée avec transform, et 0 image a été restaurée avec ʻinverse_transform`.

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

Le temps de calcul est le suivant dans l'ordre 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

Résultat de la restauration

La dimension de l'image 1 a été réduite et restaurée de la même manière.

Code source 2DPCA

Le bloc-notes exécuté se trouve sur 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, "Introduction à l'analyse des données de la méthode du noyau avec un noyau à valeur fixe positive-", Asakura Shoten, 2010.