主成分分析(Principal Components Analysis, PCA)

主要思想

PCA将$D$维特征$\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N]\in {\mathbb{R}}^{D\times N}$（${\mathbf{x}}_i\in {\mathbb{R}}^D$）映射到$d(d\ll D)$维空间中（$\mathbf{Y}=[{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_N]={\mathbf{W}}^T\mathbf{X}\in {\mathbb{R}}^{d\times N}$, $\mathbf{W}=[{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d]\in {\mathbb{R}}^{D\times d}$，${\mathbf{w}}_i\in {\mathbb{R}}^D$），使得在降维后的空间中，特征的方差最大，即保留主成分。

推导方法

根据方差最大，可以确定优化目标为
$$\begin{array}{c}\underset{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}{\mathrm{argmax}}\frac{1}{N}\sum _{j=1}^d\sum _{i=1}^N({\mathbf{w}}_j^T{\mathbf{x}}_i-{\mathbf{w}}_j^T\overline{\mathbf{x}}{)}^2\end{array}$$
其中
$$\overline{\mathbf{x}}=\frac{1}{N}\sum _{i=1}^N{\mathbf{x}}_i$$
则，可进一步将优化目标化简为矩阵形式：
$$\begin{array}{c}\begin{array}{rl}\underset{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}{\mathrm{argmax}}& \frac{1}{N}\sum _{j=1}^d\sum _{i=1}^N({\mathbf{w}}_j^T{\mathbf{x}}_i-{\mathbf{w}}_j^T\overline{\mathbf{x}}{)}^2\\ \underset{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}{\mathrm{argmax}}& \sum _{j=1}^d{\mathbf{w}}_j^T\left(\frac{1}{N}\sum _{i=1}^N({\mathbf{x}}_i-\overline{\mathbf{x}})({\mathbf{x}}_i-\overline{\mathbf{x}}{)}^T\right){\mathbf{w}}_j\end{array}\end{array}$$
令$\tilde{\mathbf{X}}=\left[\begin{array}{c}{\mathbf{x}}_1-\overline{\mathbf{x}},{\mathbf{x}}_2-\overline{\mathbf{x}},\cdots ,{\mathbf{x}}_N-\overline{\mathbf{x}}\end{array}\right]\in {\mathbb{R}}^{D\times N}$，则
$$\begin{array}{c}\begin{array}{rl}\underset{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}{\mathrm{argmax}}& \sum _{j=1}^d{\mathbf{w}}_j^T\left(\frac{1}{N}\sum _{i=1}^N({\mathbf{x}}_i-\overline{\mathbf{x}})({\mathbf{x}}_i-\overline{\mathbf{x}}{)}^T\right){\mathbf{w}}_j\\ \underset{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}{\mathrm{argmax}}& \sum _{j=1}^d{\mathbf{w}}_j^T\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}{\mathbf{w}}_j\\ \underset{{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_d}{\mathrm{argmax}}& trace({\mathbf{W}}^T\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}\mathbf{W})\end{array}\end{array}$$
一般来说引入${\mathbf{W}}^T\mathbf{W}=\mathbb{I}$的约束，这样能确保只关注${\mathbf{w}}_i$的方向而不是尺度大小，并且对任意$i\ne j$${\mathbf{w}}_i$和${\mathbf{w}}_j$都是不相关的，从而能够保留更多的原数据信息，使特征映射不退化到更低维的空间。新的优化目标为
$$\begin{array}{c}\begin{array}{rl}\underset{\mathbf{W}}{\mathrm{argmax}}& trace({\mathbf{W}}^T\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}\mathbf{W})\\ s.t.& {\mathbf{W}}^T\mathbf{W}=\mathbb{I}\end{array}\end{array}$$
采用拉格朗日乘子的方法求解以上优化问题，引入$\mathbf{\Lambda }\in {\mathbb{R}}^{d\times d}$，由于约束条件是对称的，因此$\mathbf{\Lambda }$也满足${\mathbf{\Lambda }}^T=\mathbf{\Lambda }$，则
$$L\left(\mathbf{W},\mathbf{\Lambda }\right)=-trace\left({\mathbf{W}}^T\frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}\mathbf{W}\right)+trace\left(\mathbf{\Lambda }\left({\mathbf{W}}^T\mathbf{W}-\mathbb{I}\right)\right)$$
则
$$\begin{gathered} \frac{\partial L\left(\mathbf{W},\mathbf{\Lambda }\right)}{\partial \mathbf{W}}=-2\frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}\mathbf{W}+2\mathbf{W\Lambda }=0 \\ \frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}\mathbf{W}=\mathbf{W\Lambda } \\ \left[\frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}{\mathbf{w}}_1,\cdots ,\frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}{\mathbf{w}}_d\right]=\left[{\mathbf{w}}_1,\cdots ,{\mathbf{w}}_d\right]\left[\begin{array}{cccc}{\lambda }_{11}& {\lambda }_{12}& \cdots & {\lambda }_{1d}\\ {\lambda }_{21}& {\lambda }_{22}& \cdots & {\lambda }_{2d}\\ ⋮& ⋮& \ddots & ⋮\\ {\lambda }_{d1}& {\lambda }_{d2}& \cdots & {\lambda }_{dd}\end{array}\right] \end{gathered}$$
$$\begin{gathered} \frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}{\mathbf{w}}_i=\left[{\mathbf{w}}_1,\cdots ,{\mathbf{w}}_d\right]\left[\begin{array}{c}{\lambda }_{1i}\\ ⋮\\ {\lambda }_{di}\end{array}\right] \\ {\mathbf{w}}_i^T\frac{{\tilde{\mathbf{X}}\tilde{\mathbf{X}}}^T}{N}{\mathbf{w}}_i={\lambda }_{ii} \end{gathered}$$
以上优化问题的最终结果，实质上就是$\frac{\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^T}{N}$的最大的$d$特征值（${\lambda }_{11},{\lambda }_{22},\cdots ,{\lambda }_{dd}$）所对应的特征向量。

import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets
from sklearn.decomposition import PCA

def plot_2d(X_r, y, target_names, name):
    plt.subplot(1,2,1 if 'Sklearn' in name else 2)
    colors = ["navy", "turquoise", "darkorange"]
    lw = 2

    for color, i, target_name in zip(colors, [0, 1, 2], target_names):
        plt.scatter(
            X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=0.8, lw=lw, label=target_name
        )
    plt.legend(loc="best", shadow=False, scatterpoints=1)
    plt.xlabel(f"{name} of IRIS dataset")

def sklearn_pca(X, y, target_names):
    pca = PCA(n_components=2)
    X_r = pca.fit(X).transform(X)

    print(pca.explained_variance_)

    plot_2d(X_r, y, target_names, 'Sklearn PCA')

    return X_r

def my_pca(X, y, target_names):
    n_components=2
    # 去中心化
    N = X.shape[0]
    X_mean = np.mean(X, axis=0)
    X_ = X - X_mean
    # 构建协方差矩阵
    XX = 1. / N * np.matmul(X_.T, X_)
    # 求特征向量
    values, vectors = np.linalg.eig(XX)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,n_components+1)]]
    
    X_r = np.matmul(X_, W)

    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    plot_2d(X_r, y, target_names, 'My Imple. PCA')

    return X_r



if __name__ == '__main__':
    iris = datasets.load_iris()

    X = iris.data
    y = iris.target
    target_names = iris.target_names

    plt.figure()

    X_r1 = sklearn_pca(X, y, target_names)

    X_r2 = my_pca(X, y, target_names)

    print(np.sum(np.abs(X_r1 - X_r2)))

    plt.show()

说明自己实现的方法与sklearn内核一致。

核技巧

利用$\varphi$将${\mathbf{x}}_i\in {\mathbb{R}}^D$映射到高维度空间，得到${\mathbf{X}}_{\varphi }=[\varphi ({\mathbf{x}}_1),\varphi ({\mathbf{x}}_2),\cdots ,\varphi ({\mathbf{x}}_N)]$，在该高维空间中进行PCA操作。首先，去中心化
$$\widetilde{{\mathbf{X}}_{\varphi }}=\left[\varphi ({\mathbf{x}}_1)-\overline{\varphi (\mathbf{x})},\varphi ({\mathbf{x}}_2)-\overline{\varphi (\mathbf{x})},\cdots ,\varphi ({\mathbf{x}}_N)-\overline{\varphi (\mathbf{x})}\right]$$
其中
$$\overline{\varphi (\mathbf{x})}=\frac{1}{N}\sum _{i=1}^N\varphi ({\mathbf{x}}_i)$$
则
$$\begin{array}{rl}\frac{\widetilde{{\mathbf{X}}_{\varphi }}{\widetilde{{\mathbf{X}}_{\varphi }}}^T}{N}& =\frac{1}{N}\sum _{i=1}^N(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j))(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j){)}^T\end{array}$$
需计算以上矩阵的特征向量${\mathbf{w}}_k$（对应第$k$大的特征值为${\lambda }_k$所对应的特征向量）
$$\begin{array}{c}\begin{array}{rl}\frac{1}{N}\sum _{i=1}^N(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j))(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j){)}^T{\mathbf{w}}_k& ={\lambda }_k{\mathbf{w}}_k\end{array}\end{array}$$
所以
$$\begin{array}{rl}{\mathbf{w}}_k& =\left[\sum _{i=1}^N\frac{\varphi ({\mathbf{x}}_i{)}^T{\mathbf{w}}_k-\frac{1}{N}{\sum }_{j=1}^N\varphi ({\mathbf{x}}_j{)}^T{\mathbf{w}}_k}{N{\lambda }_k}\varphi ({\mathbf{x}}_i)\right]\\ & =\sum _{i=1}^N{\alpha }_{ki}\varphi ({\mathbf{x}}_i)\end{array}$$
其中${\alpha }_{ki}=\left(\varphi ({\mathbf{x}}_i{)}^T{\mathbf{w}}_k-\frac{1}{N}{\sum }_{j=1}^N\varphi ({\mathbf{x}}_j{)}^T{\mathbf{w}}_k\right)/(N{\lambda }_k)$，则下面将求出相应的${\alpha }_{ki}$，将上式带入公式（5），可得
$$\begin{array}{c}\begin{array}{rl}\frac{1}{N}\sum _{i=1}^N(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j))(\varphi ({\mathbf{x}}_i)-\frac{1}{N}\sum _{j=1}^N\varphi ({\mathbf{x}}_j){)}^T\sum _{i=1}^N{\alpha }_{ki}\varphi ({\mathbf{x}}_i)& ={\lambda }_k\sum _{i=1}^N{\alpha }_{ki}\varphi ({\mathbf{x}}_i)\\ \frac{1}{N}({\mathbf{X}}_{\varphi }{\mathbf{X}}_{\varphi }^T-\frac{1}{N}{\mathbf{X}}_{\varphi }{1}_{N\times 1}{1}_{N\times 1}^T{\mathbf{X}}_{\varphi }^T){\mathbf{X}}_{\varphi }{\alpha }_k& ={\lambda }_k{\mathbf{X}}_{\varphi }{\alpha }_k\\ \frac{1}{N}({\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k-{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }\frac{1_{N\times 1}{1}_{N\times 1}^T}{N}{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k)& ={\lambda }_k{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{\alpha }_k\end{array}\end{array}$$
而核矩阵$\mathbf{K}={\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }$，所以
$$\begin{array}{rl}\frac{1}{N}(\mathbf{KK}{\alpha }_k-\mathbf{K}\frac{1_{N\times 1}{1}_{N\times 1}^T}{N}\mathbf{K}{\alpha }_k)& ={\lambda }_k\mathbf{K}{\alpha }_k\\ \frac{1}{N}\left(\mathbf{I}-\frac{1_{N\times N}}{N}\right)\mathbf{K}{\alpha }_k& ={\lambda }_k{\alpha }_k\\ \tilde{\mathbf{K}}{\alpha }_k& ={\lambda }_k{\alpha }_k\end{array}$$
那么可以求出$\tilde{\mathbf{K}}$的第$k$大的特征值${\lambda }_k$所对应的特征向量${\alpha }_k$，从而得到
$${\mathbf{w}}_k={\mathbf{X}}_{\varphi }{\alpha }_k$$
对一个新的数据${\mathbf{x}}_{new}$，降维到${\mathbf{w}}_k$的维度
$$\begin{array}{rl}{\mathbf{w}}_k^T\left(\varphi \left({\mathbf{x}}_{new}\right)-\frac{1}{N}{\mathbf{X}}_{\varphi }{1}_{N\times 1}\right)& ={\alpha }_k^T{\mathbf{X}}_{\varphi }^T\varphi \left({\mathbf{x}}_{new}\right)-\frac{1}{N}{\alpha }_k^T{\mathbf{X}}_{\varphi }^T{\mathbf{X}}_{\varphi }{1}_{N\times 1}\\ & ={\alpha }_k^T\mathbf{K}\left(\cdot ,{\mathbf{x}}_{new}\right)-\frac{1}{N}{\alpha }_k^T\mathbf{K}{1}_{N\times 1}\end{array}$$

import matplotlib.pyplot as plt
import numpy as np

from sklearn import datasets
from sklearn.decomposition import PCA, KernelPCA
from sklearn.datasets import make_circles
from sklearn.model_selection import train_test_split

def plot_2d(X, y, X_pca, X_kpca, X_my_kpca):
    fig, (orig_data_ax, pca_proj_ax, kernel_pca_proj_ax, kernel_pca_proj_ax2) = plt.subplots(
    ncols=4, figsize=(17, 4)
    )

    orig_data_ax.scatter(X[:, 0], X[:, 1], c=y)
    orig_data_ax.set_ylabel("Feature #1")
    orig_data_ax.set_xlabel("Feature #0")
    orig_data_ax.set_title("Testing data")

    pca_proj_ax.scatter(X_pca[:, 0], X_pca[:, 1], c=y)
    pca_proj_ax.set_ylabel("Principal component #1")
    pca_proj_ax.set_xlabel("Principal component #0")
    pca_proj_ax.set_title("my Imple. PCA")

    kernel_pca_proj_ax.scatter(X_kpca[:, 0], X_kpca[:, 1], c=y)
    kernel_pca_proj_ax.set_ylabel("Principal component #1")
    kernel_pca_proj_ax.set_xlabel("Principal component #0")
    _ = kernel_pca_proj_ax.set_title("Sklearn KernelPCA")

    kernel_pca_proj_ax2.scatter(X_my_kpca[:, 0], X_my_kpca[:, 1], c=y)
    kernel_pca_proj_ax2.set_ylabel("Principal component #1")
    kernel_pca_proj_ax2.set_xlabel("Principal component #0")
    _ = kernel_pca_proj_ax2.set_title("my Imple. KernelPCA")

    plt.show()

def my_pca(X_tr, X_te):
    n_components=2
    # 去中心化
    N = X_tr.shape[0]
    X_mean = np.mean(X_tr, axis=0)
    X_ = X_tr - X_mean
    # 构建协方差矩阵
    XX = 1. / N * np.matmul(X_.T, X_)
    # 求特征向量
    values, vectors = np.linalg.eig(XX)
    idxs = np.argsort(values)
    
    W = vectors[:,[idxs[-i] for i in range(1,n_components+1)]]
    
    X_r = np.matmul(X_te - X_mean, W)

    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    return X_r

def my_kpca(X_tr, X_te, gamma=10):
    # exp(-||x_i-x_j||^2*gamma)
    N = X_tr.shape[0]
    K_tr0 = np.sum(X_tr**2, axis=1, keepdims=True)
    K_tr0 = np.exp(-(K_tr0 + K_tr0.T - 2 * np.matmul(X_tr, X_tr.T))*gamma)
    K_ = 1. / N * (K_tr0 - np.ones([N, N]) @ K_tr0 /N)
    # 求特征向量
    values, vectors = np.linalg.eig(K_)
    idxs = np.argsort(values)
    
    alphas = vectors[:,[idxs[-i] for i in range(1,3)]] # N_tr x 2
    K_new = np.sum(X_tr**2, axis=1, keepdims=True) + \
            np.sum(X_te**2, axis=1, keepdims=True).T - \
            2 * X_tr @ X_te.T
    K_new = np.exp(-K_new * gamma) # N_tr x N_te
    X_r = K_new.T @ alphas - 1. / N * np.ones([1, N]) @ K_tr0 @ alphas

    X_r[:,0] = - X_r[:,0] #确保与sklearn得到的结果一致
    X_r[:,1] = - X_r[:,1] #确保与sklearn得到的结果一致

    return X_r

def sklearn_kpca(X_tr, X_te, gamma=10):
    kernel_pca = KernelPCA(
        n_components=None, kernel="rbf", gamma=gamma, fit_inverse_transform=False
    )
    X_test_kernel_pca = kernel_pca.fit(X_train).transform(X_test)
    return X_test_kernel_pca



if __name__ == '__main__':
    X, y = make_circles(n_samples=1000, factor=0.3, noise=0.05, random_state=0)
    X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=0)

    X_pca = my_pca(X_train, X_test)
    X_kpca = sklearn_kpca(X_train, X_test, gamma=10.0)
    X_my_kpca = my_kpca(X_train, X_test, gamma=10.0)

    plot_2d(X_test, y_test, X_pca, X_kpca, X_my_kpca)

虽然与sklearn的KPCA在倍数上有点出入，但是不想深究了，就这样吧！