12.2 基于核函数的主成分分析#
在12.1节内容中我们详细介绍了PCA算法的思想原理，在本节内容中将会介绍另外一种基于PCA算法改进的核主成分分析算法（Kernel Principal Component Analysis, KPCA）[1]。如图12-12便是基于核函数的主成分分析算法的学习路线图，首先仍旧是先介绍KPCA背后的思想和在sklearn中的用法，然后介绍KPCA的原理与求解过程，最后介绍如何从零实现KPCA算法。
12.2.1 KPCA算法思想#
从名字来看，KPCA算法是在原有的PCA降维算法中引入了核技巧（Kernel Trick），而所谓核技巧则是指能够快速计算得到两个低维向量映射到高维（甚至是无穷维）空间后内积的方法，详细内容可以阅读本书10.8节内容。之所以需要在原始PCA算法中引入特征映射，是因为PCA算法只能够处理线性分布的数据样本，而对于非线性的场景来说PCA算法将会显得无能为力。
如图12-13所示，在二维空间中有一个非线性可分的数据集，其一共包含有两个类别。假设现在通过PCA算法来对其进行降维处理将会产生什么样的结果呢？
根据PCA算法原理可知，图12-13中的数据样本无论是朝着哪个方向进行投影降维，在子空间中都会将两个类别的样本点完全混合在一起，而这样的降维操作显然是无效的。因为降维算法的目的是在尽可能多的保留原有数据结构信息的前提条件下降低特征维度，而不仅仅只是为了单方面的降低维度而过多（或全部）丢失了数据原有的结构信息。如图12-14所示便是图12-13中的样本通过PCA算法降维后的结果。
如图12-14所示，可以看到经过PCA算法降维后的结果已经完全不能够体现原有的样本结构信息（这里特指类别）。因此，这就需要借助特征映射的思想先将原始样本点映射到一个高维的线性可分的特征空间中，然后再通过PCA算法对齐进行降维处理。如图12-15所示便是经过KPCA算法降维后的结果。
如图12-15所示，可以看到经过KPCA算法降维后的结果依旧能够体现出两个完全不同的样本类别。
根据上面的介绍可知，KPCA算法的核心思想主要可以分为两步：第1步，先将原始非线性可分的样本点映射到线性可分的高维或无穷维（高斯核）的特征空间中；第2步，在线性可分的高维空间中通过PCA算法进行降维处理。
如图12-6所示，可以选择某个合适的核函数将原始非线性可分的样本点映射到线性可分的高维空间中。
然后再通过PCA算法对其进行降维处理便可以得到类似图12-15中的结果，如图12-17所示。
如图12-17所示，左边为PCA降维后使用2个主成分表示的可视化结果，右边为使用1个主成分表示的结果。上述可视化代码可参见 AllBooKCode/Chapter12/C05_KPCA_idea.py 文件。
12.2.2 KPCA示例代码#
在清楚了KPCA算法的基本思想后，我们便可以通过sklearn中的KernelPCA模块来进行使用，示例代码如下：
1 import numpy as np
2 from sklearn.decomposition import KernelPCA
3 from sklearn.datasets import make_circles
4 import matplotlib.pyplot as plt
5
6 def make_nonlinear_cla_data(num_points = 500):
7 x, y = make_circles(n_samples=num_points, factor=0.2, noise=0.1,
8 random_state=np.random.seed(10))
9 x = x.reshape(-1, 2)
10 return x, y.reshape(-1, 1)在上述代码中，第1~4行是导入相关模块。第6~10行是构造一个线性不可分的二分类数据集。
进一步，可以通过KernelPCA模块来对其进行降维处理并可视化，示例代码如下：
1 def visualization():
2 x, y = make_nonlinear_cla_data()
3 pca = KernelPCA(n_components=2, kernel='rbf', gamma=2.)
4 plt.figure(figsize=(12, 4))
5 plt.subplot(1, 3, 1)
6 plt.title('Original', fontsize=15)
7 plt.scatter(x[:, 0], x[:, 1], c=y)
8 x = pca.fit_transform(x)
9 plt.subplot(1, 3, 2)
10 plt.title('Projection with two components', fontsize=15)
11 plt.scatter(x[:, 0], x[:, 1], c=y)
12 plt.subplot(1, 3, 3)
13 plt.title('Projection with one component', fontsize=15)
14 plt.scatter(x[:, 0], [0] * len(x), c=y)
15 plt.tight_layout()
16 plt.show()
17
18 if __name__ == '__main__':
19 visualization()在上述代码中，第3行是实例化一个KernelPCA类对象，其中n_components用于指定主成分的数量，kernel用于指定对应的核函数，gamma为高斯核函数中的参数。第4~14行是分别对未降维、用两个主成分和一个主成分表示的结果进行可视化。
最后，上述代码运行结束后将会得到如图12-18所示的结果。
上述可视化及示例代码可参见 AllBooKCode/Chapter12/C06_KPCA_train.py 文件。
12.2.3 KPCA算法原理#
在清楚了KPCA算法的思想之后，其背后的数学原理就变得十分清晰了。根据PCA算法的原理可知，对于任意数据集$X=\{{x^{(1)},x^{(2)},...,x^{(m)}}\}$来说需要求解的方程为
$$ \Sigma u = \lambda u\tag{12-16} $$其中
$$ \Sigma=\frac{1}{m}\sum_{i=1}^mx^{(i)}x^{(i)^T}=\frac{1}{m}X^T_{m\times n}X_{m\times n}\tag{12-17} $$最后根据协方差矩阵$\Sigma$求得其对应的特征值和特征向量，并选择前$k$个最大的特征值所对应的特征向量作为主成分进行降维即可。
对于KPCA算法来说，首先需要定义一个映射函数 $\phi(x)$ 将原始的低维向量 $x^{(i)}\in R^n$ 映射为高维向量 $\phi(x^{(i)})\in R^N$ ，其中 $n\ll N$ 。
假设在映射后的高维样本空间中，所有样本点的均值为0，即
$$ \mu=\frac{1}{m}\sum_{i=1}^m\phi(x^{(i)})=0\tag{12-18} $$根据式(12-17)可知，此时新的协方差矩阵为
$$ \overline{\Sigma}=\frac{1}{m}\sum_{i=1}^m\phi(x^{(i)})\phi(x^{(i)})^T\tag{12-19} $$进一步，根据式(12-11)可得
$$ \overline{\Sigma}v=\lambda v\tag{12-20} $$最后，只需要根据协方差矩阵$\overline{\Sigma}$求得对应的特征值和特征向量，并同样选择前$k$个最大的特征值所对应的特征向量作为主成分进行降维即可。同时，这里需要注意的是，当高维空间$N$的取值非常大时，直接求解协方差矩阵$\overline{\Sigma}$的特征值与特征向量将是一件非常困难的事情，因此就需要借助核函数来完成整个计算过程，这部分内容将在12.2.5内容中节进行介绍。不过核函数只是一个计算上的技巧，它并不影响对于KPCA算法的整体把握。
12.2.4 KPCA算法求解#
在介绍完KPCA算法的基本原理和在sklearn中的使用示例后，下面再来介绍KPCA算法的求解问题。
根据式(12-19)和式(12-20)可知
$$ \overline{\Sigma}v=\frac{1}{m}\sum_{i=1}^m\phi(x^{(i)})\phi(x^{(i)})^Tv=\lambda v\tag{12-21} $$进一步由式(12-21)可得
$$ v=\frac{1}{\lambda m}\sum_{i=1}^m\phi(x^{(i)})\phi(x^{(i)})^Tv=\frac{1}{\lambda m}\sum_{i=1}^m\left(\phi(x^{(i)})^T v\right)\phi(x^{(i)})^T\tag{12-22} $$其中对于式(12-22)右侧的变形过程如下
设$x=(x_1,x_2,...,x_n)^T$，$v=(v_1,v_2,...,v_n)^T$，则此时有
$$ \begin{aligned} xx^Tv &=\begin{pmatrix}x_1\\x_2 \\ \vdots \\ x_n \end{pmatrix}_{n\times1} \begin{pmatrix} x_1&x_2&\cdots&x_n\end{pmatrix}_{1\times n} \begin{pmatrix}v_1\\v_2 \\ \vdots \\ v_n \end{pmatrix}_{n\times 1} = \begin{pmatrix}x_1x_1v_1+x_1x_2v_2+&\cdots&+x_1x_nv_n\\x_2x_1v_1+x_2x_2v_2+&\cdots&+x_2x_nv_n \\ \vdots \\x_nx_1v_1+x_nx_2v_2+&\cdots&+x_nx_nv_n\end{pmatrix}_{n\times 1} \\[3ex] &= \begin{pmatrix}(x_1v_1+x_2v_2+&\cdots&+x_nv_n)x_1\\(x_1v_1+x_2v_2+&\cdots&+x_nv_n)x_2 \\ \vdots \\(x_1v_1+x_2v_2+&\cdots&+x_nv_n)x_n\end{pmatrix}_{n\times 1}=x^Tvx^T \end{aligned}\tag{12-23} $$