第5章：通过降维压缩数据

GitHub Notebook 地址: http://nbviewer.ipython.org/github/rasbt/python-machine-learning-book/blob/master/code/ch05/ch05.ipynb

在第4章中，我们介绍了通过不同的特征选择技术对数据集进行降维的方法。另一种常用于降维的特征选择方法就是特征抽取。在本章中，读者将学习三种可以帮助我们归纳总结数据集内所蕴含信息的技术，它们都可以将原始数据集变换到一个维度更低的新的特征子空间。数据压缩也是机器学习领域中的一个重要内容，随着现代技术的发展，将会产生越来越多的数据，数据压缩技术可以帮助我们队数据进行存储和分析。本章将涵盖如下主题:

• 无监督数据压缩——主成分分析(Principal Component Analysis, PCA)
• 基于类别可分最大化的监督降维技术——线性判别分析(Linear Discriminant Analysis, LDA)
• 通过核主成分分析(kernel principal component analysis) 进行非线性降维

1. 无监督数据降维技术——主成分分析

与特征选择类似，我们可以使用特征抽取来减少数据集中特征的数量。不过，当使用序列后向选择等特征选择算法时，能够保持数据的原始特征，而特征抽取算法则会将数据转换或者映射到一个新的特征空间。基于降维在数据预处理领域的含义，特征抽取可以理解为: 在尽可能多地保持相关信息的情况下，对数据进行压缩的一种方法。特征抽取通常用于提高计算效率，同样也可以帮助我们降低“维度灾难”——尤其当模型不适于正则化处理时。

主成分分析(principal component analysis, PCA) 是一种广泛应用于不同领域的无监督线性数据转换技术，其突出作用是降维。PCA 的其他常用领域包括: 股票交易市场的探索性分析和信号去噪，以及生物信息学领域的基因组和基因表达水平数据分析等。PCA 可以帮助我们基于特征之间的关系识别出数据内在的模式。简而言之，PCA 的目标是在高维数据中找到最大方差的方向，并将数据映射到一个维度不大于原始数据的新的子空间上。如下图所示，以新特征的坐标是相互正交为约束条件，新的子空间上正交的坐标轴(主成分)可被解释为方差最大的方向。在此，\(x_1\) 和 \(x_2\) 为原始特征的坐标轴，而 PC1 和 PC2 即为主成分。

如果使用 PCA 降维，我们将构建一个 \(d \times k\)的转换矩阵 \(W\)。这样就可以将一个样本向量 \(x\) 映射到一个新的 \(k\) 维特征子空间上去，此空间的维度小于原始的 \(d\) 维特征空间:

\[ \begin{align} x &= [x_1, x_2, ..., x_d], \quad x \in R^d \\ & \downarrow x^W, \quad W \in R^{d \times k} \\ z &= [z_1, z_2, ..., z_k], \quad z \in R^k \end{align} \]

完成从原始 \(d\) 维数据到新的 \(k\) 维子空间(一般情况下 \(k << d\)) 的转换后，第一主成分的方差应该是最大的，由于各大主成分之间是不相关的(正交的)，后续各主成分也具备尽可能大方差。需注意的是，主成分的方向对数据值的范围高度敏感，如果特征的值使用不同的度量标准，我们需要先对特征进行标准化处理，以让各特征具有相同的重要性。

在详细讨论使用 PCA 算法降维之前，我们先通过以下几个步骤来概括以下算法的流程:

1. 对原始 \(d\) 维数据集做标准化处理。
2. 构造样本的协方差矩阵。
3. 计算协方差矩阵的特征值和相应的特征向量。
4. 选择与前 \(k\) 个最大特征值对应的特征向量，其中 \(k\) 为新特征空间的维度(\(k \leq d\))。
5. 通过前 \(k\) 个特征向量构造映射矩阵 \(W\)。
6. 通过映射矩阵 \(W\) 将 d 维的输入数据集 \(X\) 转换到新的 \(k\) 维特征子空间。

1.1 总体方差与贡献方差

在本小节中，我们将学习主成分分析算法的前四个步骤: 数据标准化、构造协方差矩阵、获得协方差举着你的特征值和特征向量，以及按降序排列特征值对所对应的特征向量，以及按降序排列特征值所对应的特征向量:

首先，我们加载第4章已经使用过的葡萄酒数据集:

import pandas as pd
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)

接下来，我们将葡萄酒数据集划分为训练集和测试集——分别占数据集的 70% 和 30%，并使用单位方差将其标准化。

from sklearn.cross_validation import train_test_split
from sklearn.preprocessing import StandardScaler

X, y = df_wine.loc[:, 1:].values, df_wind.loc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
sc = StandardScaler()

X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit_transform(X_test)

通过上述代码完成数据预处理后，我们进入第二步: 构造协方差矩阵。此 \(d \times d\) 维协方差矩阵是沿主对角线对称的，其中 \(d\) 为数据集的维度，此谨遵成对地存储了不同特征之间的协方差。例如，对群体进行描述的两个特征 \(x_j\) 和 \(x_k\) 可通过如下公式计算它们之间的协方差:

\[\sigma_{jk} = \frac {1}{n} \sum_{i=1}^n (x_j^{(i)} - \mu_j)(x_k^{(i)} - \mu_k)\]

在此，\(\mu_j\) 和 \(\mu_k\) 分别为特征 \(j\) 和 \(k\) 的均值。请注意，我们对数据集做了标准化处理，样本的均值将为零。两个特征之间的协方差如果为正，说明它们会同时增加增减，而一个负的协方差值则表示两个特征会朝相反的方向变动。例如，一个包含三个特征的协方差矩阵可记为(注意: 此处 \(\Sigma\) 代表希腊字母 simga，在此请勿求和符号混为一谈):

\[ \Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_2^2 & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_3^2 \end{bmatrix} \]

协方差矩阵的特征向量代表主成分(最大方差方向)，而对应的特征值大小就决定了特征向量的重要性。就葡萄酒数据集看来说，我们可以得到 \(13 \times 13\) 维协方差矩阵的13个特征向量及其对应的特征值。

现在，我们来计算协方差矩阵的特征对。通过线性代数或微积分的相关知识我们知道，特征值 V 需满足如下条件:

$$ \Delta \nu = \lambda \nu $$ 此处的特征值 \(\lambda\) 是一个标量。由于手动计算特征向量和特征值从某种程序上来说，是一项繁琐且复杂的工作，我们将使用 NumPy 中的 linalg.eig 函数来计算葡萄酒数据集协方差矩阵的特征对:

>>> import numpy as np
>>> cov_mat = np.cov(X_train_std.T)
>>> eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
>>> print('\nEigenvalues \n%s' % eigen_vals)
Eigenvalues 
[4.8923083  2.46635032 1.42809973 1.01233462 0.84906459 0.60181514
 0.52251546 0.08414846 0.33051429 0.29595018 0.16831254 0.21432212
 0.2399553 ]

应用 numpy.cov 函数，我们计算得到了经标准化处理的训练数据集的协方差矩阵。使用 linalg.eig 函数，通过特征分解，得到一个包含 13 个特征值的向量(eigen_vals)，及其对应的特征值，特征向量以列的方式存储一个 \(13 \times 13\) 维的矩阵中(eigen_vecs)。

因为要将数据压缩到一个新的特征子空间上来实现数据降维，所以我们只选择那些包含最多信息(方差最大)的特征向量(主成分)组成子集。由于特征值的大小决定了特征向量的重要性，因此需要将特征值按降序排列，我们感兴趣的是排序在前 \(k\) 个的特征值所对应的特征向量。在整理包含信息最大的前 \(k\) 个特征向量前，我们先回执特征值的方差贡献率(variable explained ratios) 图像。

特征值 \(\lambda_j\) 的方差贡献率是指，特征值 \(\lambda_j\) 与所有特征值和的比值:

$$ \frac {\lambda_j} {\sum_{j=1}^d \lambda_j} $$ 使用 NumPy 的 cumsum 函数，我们可以计算出累计方差，其图像可通过 matplotlib 的 step 函数绘制:

tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
import matplotlib.pyplot as plt
plt.bar(range(1, 14), var_exp, alpha=0.5, align='center', 
        label='individual explained variabce')
plt.step(range(1, 14), cum_var_exp, where='mid',
    label='cumulative explained variable')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.show()

由下图可以看到，第一主成分占方差总和的 40% 左右；此外，还可以看出前两个主成分占总体方差的近 60%:

虽然方差贡献率图像可以让我们联想到第4章中通过随机森林计算出的关于特征的重要程度，但我们应该注意: PCA 是一种无监督方法，这意味着我们可以忽略类标信息。相对而言，随机森林通过类标信息来计算及节点不纯度，而方差度量的是特征值在轴线上的分布。

1.2 特征转换

在将方差矩阵分解为特征对后，我们继续执行 PCA 方法的最后三个步骤，将葡萄酒数据集中的信息转换到新的主成分轴上。在本小节中，我们将对特征值按降序进行排列，并通过挑选出对应的特征向量构造出映射矩阵，然后使用映射矩阵转换到低维的子空间上。

首先，按特征值的降序排列特征对:

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) 
    for i in range(len(eigen_vals))]
eigen_pairs.sort(reverse=True)

接下来，我们选取两个对应的特征值最大的特征向量，这两个值之和占据了数据集总体方差的 60%。请注意，出于演示的需要，我们只选择了两个特征向量，因为在本小节中我们将以二维散点图的方式回执相关图像。在实际应用中，确定主成分的数量时，我们需要根据实际情况在计算效率与分类器性能之间做出权衡。

>>> w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
            eigen_pairs[1][1][:, np.newaxis]))
>>> print('Matrix W: \n', w)
Matrix W: 
 [[ 0.14669811  0.50417079]
 [-0.24224554  0.24216889]
 [-0.02993442  0.28698484]
 [-0.25519002 -0.06468718]
 [ 0.12079772  0.22995385]
 [ 0.38934455  0.09363991]
 [ 0.42326486  0.01088622]
 [-0.30634956  0.01870216]
 [ 0.30572219  0.03040352]
 [-0.09869191  0.54527081]
 [ 0.30032535 -0.27924322]
 [ 0.36821154 -0.174365  ]
 [ 0.29259713  0.36315461]]

执行上述代码，通过选取的两个特征向量，我们得到了一个 \(13 \times 12\) 维的映射矩阵 \(W\)。通过映射矩阵，我们可以将一个样本 \(x\)(以13维的行向量表示) 转换到 PCA 的子空间上得到 \(x'\)，样本转换为一个包含两个新特征的二维向量:

\[ x'=xW \]

>>> X_train_std[0].dot(w)
array([2.59891628, 0.00484089])

类似地，通过计算矩阵的点积，我们可以将 \(124 \times 13\) 维的训练数据集转换到包含两个主成分的子空间上:

\[ X'=XW \]

>>> X_train_pca = X_train_std.dot(w)

最后，转换后的葡萄酒数据集将以 \(124 \times 2\) 维矩阵的方式存储，我们以二维散点图的方式来对其进行可视化展示:

colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train == l, 0],
                X_train_pca[y_train == l, 1],
                c=c, label=l, marker=m)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()

在下图中可以看到，相较于第二主成分(y轴)，数据更多地沿着x轴(第一主成分)方向分布，这与我们在上一小节中绘制的方差贡献率图是一致的。而且，可以直观地看到，线性分类器能够很好地对其进行划分：

为了演示散点图，我们使用了类标信息，但请注意 PCA 是无监督方法，无需类标信息。

1.3 使用 scikit-learn 进行主成分分析

通过上一小节中详尽的介绍，我们了解了 PCA 内部的工作原理，接下来讨论如何使用 scikit-learn 中提供的 PCA 类。PCA 也是 scikit-learn 中的一个数据转换类，在使用相同的模型参数对训练数据和测试数据进行转换之前，我们首先使用训练数据对模型进行拟合。下面，我们使用 scikit-learn 中的 PCA 对葡萄酒数据集做预处理，然后使用逻辑斯特回归对转换后的数据进行分类，最后用第2章定义的 plot_decision_region 函数对决策区域进行可视化展示。

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):
    # setup markder generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors=  ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
            np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
            alpha=0.8, c=cmap(idx), 
            marker=markers[idx], label=cl)

from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
lr = LogisticRegression()
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
lr.fit(X_train_pca, y_train)
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()

执行上述代码后，我们可以看到将训练数据转换到两个主成分轴后生成的决策区域。

如果比较 scikit-learn 中提供的 PCA 类与我们自己实现的 PCA 类的分析结果，可以发现: 上图可以看作是我们此前自己完成 PCA 方法所得到沿 y 轴进行旋转后的结果。出现此差异的原因，不是两种方法在实现中出现了什么错误，而在于特征分析方法: 特征向量可以为正或者为负。这不是重点，因为有需要时我们可以在数据上乘以 -1 来实现图像的镜像。注意，特征向量通常会缩放到单位长度1。为了保证整个分析过程的完整性，我们绘制一下逻辑斯特回归在转换后的测试数据上所得到的决策区域，看其是否能很好地将各类分开:

plot_decision_regions(X_test_pca, y_test, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()

执行上述代码，我们绘制出了测试集的决策区域，可以看到: 逻辑斯特回归在这个小的二维特征子空间上表现优异，只误判了一个样本。

如果我们队不同主成分的方差贡献率感兴趣，可以将 PCA 类中的 n_components 参数设置为 None。由此，可以保留所有的主成分，并且可以通过 explained_variance_ratio_ 属性得到相应的方差贡献率:

>>> pca = PCA(n_components=None)
>>> X_train_pca = pca.fit_transform(X_train_std)
>>> pca.explained_variance_ratio_
array([0.37329648, 0.18818926, 0.10896791, 0.07724389, 0.06478595,
       0.04592014, 0.03986936, 0.02521914, 0.02258181, 0.01830924,
       0.01635336, 0.01284271, 0.00642076])

请注意，在初始化 PCA 类时，如果我们将 n_components 设置为 None，那么它将按照方差贡献率递减顺序返回所有的主成分，而不是进行降维操作。

2. 通过线性判别分析压缩无监督数据

线性判别分析(Linear Discriminant Analysis，LDA) 是一种可作为特征抽取的技术，它可以提高数据分析过程中的计算效率，同时，对于不适用正则化的模型，它可以降低因维度灾难带来的过拟合。

LDA 的基本概念与 PCA 非常相似，PCA 视图在数据集中找到方差最大的正交的主成分分量的轴，而 LDA 的目标是发现可以最优化分类的特征子空间。LDA 与 PCA 都是可用于降低数据集维度的线性转换技巧。其中，PCA 是无监督算法，而 LDA 是监督算法。因此，我们可以这样直观地认为: 与 PCA 相比，LDA 是一种更优越的用于分类的特征提取技术。但是 A.M.Martinez 提出: 在图像识别任务中的某些情况下，如每个类别中只有少量样本，使用 PCA 作为预处理工具的分类结果更佳\(^{注1}\)。

「🔖 LDA 有时也称作 Fisher LDA，Ronald A.Fisher 于 1937年针对二类别问题对 Fisher 线性判别(Fisher Linear Discriminant) 做了最初的形式化\(^{注2}\)。1948年，基于类别方差相等和类内样本呈现标准正太分布的假设，Radhakrishna 将 Fisher LDA 泛化到了多类别分类问题上，即我们现在所说的 LDA\(^{注3}\)。」

注1: A.M.Martinez and A.C.Kak.PCA VerSUS LDA. Pattern Analysis and Machine Intelligence，IEEE Transactions on，23(2): 228-223，2001. 注2: R.A. Fisher. The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics，7(2): 179-188，1936. 注3: C.R. Rao. The Utilization of Multiple Measurements in Problems of Biological Classification. Journal of the Royal Statistical Society. Series B(Methodological)，10(2): 159-203，1948.

下图解释了二类别分类中 LDA 的概念。类别1、类别2 中的样本分别用叉号和原点来表示:

如上图所示，在 x 轴方向(LD1)，通过现象判定，可以很好地将正态分布的两个类分开。虽然沿 y 轴(LD2) 方向的线性判定保持了数据集的较大方差，但是沿此方向无法提供关于类别区分的任何信息，因此它不是一个好的线性判定。

一个关于 LDA 的假设就是数据呈正态分布。此外，我们还嘉定各类别中数据具有相同的协方差矩阵，且样本的特征从统计上来讲是相互独立的。不过，即使一个或多个假设没有满足，LDA 仍旧可以很好地完成降维工作\(^{注4}\)。

注4: RR.O.Duda，P.E.Hart，and D.G.Stork.Pattern Classification.2^nd.Edition.New York, 2001.

在进入下一节详细讨论 LDA 的原理之前，我们先来总结一个 LDA 方法的关键步骤:

1. 对 \(d\) 维数据集进行标准化处理(d 为特征的数量)
2. 对于每一类别，计算 \(d\) 维的均值向量。
3. 构造类间的散步矩阵 \(S_B\) 以及类内的散步矩阵 \(S_W\)。
4. 计算矩阵 \(S_W^{-1}S_B\) 的特征值及对应的特征向量。
5. 选取前 \(k\) 个特征值所对应的特征向量，构造一个 \(d \times k\) 维的转换矩阵 \(W\)，其中特征向量以列的形式排列。
6. 使用转换矩阵 \(W\) 将样本映射到新的特征子空间上。

🔖 特征呈正态分布且特征间相互独立是我们使用 LDA 时所作的假设。同时，LDA 算法假定各个类别的协方差矩阵是一致的。然而，即使我们违背了上述假设，LDA 算法仍旧能很好地完成数据降维及分类任务(R.O.Duda，P.E.Hart，and D.G.Stork.Pattern Classification.2^nd.Edition.New York, 2001)。

2.1 计算散步矩阵

本节伊始，在讲解 PCA 时就对普泡酒数据集做了标准化处理，因此我们将跳过第一步直接计算均值向量。计算中，我们将分别构建内散布矩阵和类间散布矩阵。均值向量 \(m_i\) 存储了类别 \(i\) 中样本的特征均值 \(\mu_m\):

$$ m_i = \frac {1}{n} \sum_{x \in D_i}^c x_m $$ 葡萄酒数据集的三个类别对应三个均值向量:

\[ m_i = \begin{bmatrix} \mu_{i, alcohol} \\ \mu_{i, malic acid} \\ \vdots \\ \mu_{i, proline} \end{bmatrix} \quad i \in \{1,2,3\} \]

np.set_printoptions(precision=4)
mean_vecs = []
for label in range(1, 4):
    mean_vecs.append(np.mean(X_train_std[y_train == label], axis=0))
    print('MV %s: %s\n' % (label, mean_vecs[label - 1]))

# 输出:
MV 1: [ 0.9259 -0.3091  0.2592 -0.7989  0.3039  0.9608  1.0515 -0.6306  0.5354
  0.2209  0.4855  0.798   1.2017]

MV 2: [-0.8727 -0.3854 -0.4437  0.2481 -0.2409 -0.1059  0.0187 -0.0164  0.1095
 -0.8796  0.4392  0.2776 -0.7016]

MV 3: [ 0.1637  0.8929  0.3249  0.5658 -0.01   -0.9499 -1.228   0.7436 -0.7652
  0.979  -1.1698 -1.3007 -0.3912]

通过均值向量，我们来计算一下类内散布矩阵 \(S_W\):

$$ S_W = \sum_{i=1}^c S_i $$ 这可以通过累加各类别 \(i\) 的散步矩阵 \(S_i\) 来计算:

\[ S_i = \sum_{x \in D_i}^c (x-m_i)(x-m_i)^T \]

d = 13 # number of features
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.zeros((d, d))
    for row in X[y == label]:
        row, mv = row.reshape(d, 1), mv.reshape(d, 1)
        class_scatter += (row - mv).dot((row - mv).T)
    S_W += class_scatter
print('Within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

# 输出:
Within-class scatter matrix: 13x13

此前我们队散布矩阵进行计算时，曾假设训练集的类标是均匀分布的。但是，通过打印类标的数量，可以看到在此并未遵循此假设:

>>> print('Class label distribution: %s' % np.bincount(y_train)[1:])
Class label distribution: [40 49 35]

因此，在我们通过累加方式计算散布矩阵 \(S_W\) 前，需要对各类别的散步矩阵 \(S_i\) 做缩放处理。当我们用各类别单独的散布矩阵除以此类别内样本的数量 \(N_i\) 时，可以发现计算散布矩阵的方式与计算协方差矩阵 \(\sum_i\) 的方式是一致的。协方差矩阵可以看作是归一化的散布矩阵:

\[ \sum_i = \frac {1}{N_i} S_w= \frac {1}{N_i} \sum_{x \in d_i}^c (x-m_i)(x-m_i)^T \]

d = 13 # number of features
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.cov(X_train_std[y_train == label].T)
    S_W += class_scatter

print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

# 输出:
Scaled within-class scatter matrix: 13x13

在完成类内散布矩阵(或协方差矩阵)的计算后，我们进入下一步骤，计算类间散布矩阵 \(S_B\):

\[ S_B = \sum_{i=1}^c N_i (m_i - m)(m_i-m)^T \]

其中，\(m\) 为全局均值，它在计算时用到了所有类别中的全部样本:

mean_overall = np.mean(X_train_std, axis=0)
d = 13 # number of features
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vecs):
    n = X[y == i+1, :].shape[0]
    mean_vec = mean_vec.reshape(d, 1)
    mean_overall = mean_overall.reshape(d, 1)
    S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
print('Between-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

# 输出:
Between-class scatter matrix: 13x13

2.2 在新特征子空间上选取线性判别算法

LDA 余下的步骤与 PCA 的步骤相似。不过，这里我们部队协方差矩阵做特征做特征分解，而是求解矩阵 \(S_W^{-1}S_B\) 广义特征值:

eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

在求得特征之后，我们按照降序对特征值进行排序:

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) 
                        for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
print('Eigenvalues in decreasing order:\n')
for eigen_val in eigen_pairs:
    print(eigen_val[0])

# 输出:
Eigenvalues in decreasing order:

643.015384346051
225.0869818541626
8.792307609971557e-14
5.250276792004119e-14
5.2491045182795644e-14
5.2491045182795644e-14
5.062616992290714e-14
4.281024234637359e-14
4.281024234637359e-14
1.1103061461361521e-14
5.3026791025565126e-15
5.3026791025565126e-15
4.953358662764158e-15

熟悉线性代数的读者应该知道: \(d \times d\) 维协方差矩阵的秩最大为 \(d-1\)，而且确实可以发现，我们只得到了两个非零特征值(实际得到的第 3~13 个特征值并非完全为零，而且趋近于0的实数，这个结果是由 NumPy 浮点运算导致的)。请注意，在极少的情况下可达到完美的共线性(所有样本的点落在一条直线上)，这时协方差矩阵的秩为1，将导致矩阵只有一个含非零特征值的特征向量。

为了度量线性判别(特征向量)可以获取多少区分类别的信息，与前面 PCA 小节中对累积方差的绘制类似，我们按照特征值降序绘制出特征对线性判别信息保持程序的图像。为了简便起见，我们在此使用了判定类别区分能力的相关信息(discriminability)。

tot = sum(eigen_vals.real)
discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]
cum_discr = np.cumsum(discr)
plt.bar(range(1, 14), discr, alpha=0.4, align='center',
        label='individual "discriminability" ')
plt.step(range(1, 14), cum_discr, where='mid',
        label='cumulative "discriminability" ')
plt.xlabel('Linear Discriminants')
plt.ylabel('"discriminability" ratio')
plt.ylim([-0.1, 1.1])
plt.legend(loc='best')
plt.tight_layout()
plt.show()

从结果图像中可以看到，前两个线性判别几乎获取到了葡萄酒训练数据集中全部有用信息:

下面我们叠加这两个判别能力最强的特征向量列来构建转换矩阵 \(W\):

w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
            eigen_pairs[1][1][:, np.newaxis].real))
print('Matrix W:\n', w)

# 输出:
Matrix W:
 [[-0.0707  0.3778]
 [ 0.0359  0.2223]
 [-0.0263  0.3813]
 [ 0.1875 -0.2955]
 [-0.0033 -0.0143]
 [ 0.2328 -0.0151]
 [-0.7719 -0.2149]
 [-0.0803 -0.0726]
 [ 0.0896 -0.1767]
 [ 0.1815  0.2909]
 [-0.0631 -0.2376]
 [-0.3794 -0.0867]
 [-0.3355  0.586 ]]

2.3 将样本映射到新的特征空间

通过上一小节中构建的矩阵 \(W\)，我们可以通过乘积的方式对训练集进行转换:

\[ X' = XW \]

X_train_lda = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_lda[y_train == l, 0],
        X_train_lda[y_train == l, 1],
        c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='best')
plt.show()

通过结果图像可见，三个葡萄酒类在新的特征子空间上是线性可分的:

2.4 使用 scikit-learn 进行 LDA 分析

自己写代码逐步实现 LDA，对于理解 LDA 内部的工作原理及其与 PCA 的差别是一种很好的体验。下面我们来看看 scikit-learn 中对 LDA 类的实现:

# from sklearn.lda import LDA # Ver. < 0.18
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
lda = LDA(n_components=2)
x_train_lda = lda.fit_transform(X_train_std, y_train)

接下来，在将训练数据通过 LDA 进行转换后，我们来看一下逻辑斯特回归在相对低维数据上的表现:

lr = LogisticRegression()
lr = lr.fit(X_train_lda, y_train)
plot_decision_regions(X_train_lda, y_train, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='best')
plt.show()

从结果图像中可以看到，逻辑斯特回归模型只错误地判断了类别2中的一个样本:

通过降低正则化强度，我们或许可以对决策树边界进行调整，以使得逻辑斯特回归模型可以正确地对训练数据进行分类。下面，我们再来看一下模型在测试数据集上的效果:

X_test_lda = lda.transform(X_test_std)
plot_decision_regions(X_test_lda, y_test, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='best')
plt.show()

正如从结果图像中看到的那样，当我们使用只有两维的特征子空间来替代原始数据集中的13个葡萄酒特征时，逻辑斯特回归在测试数据集上对样本的分类结果可谓完美:

3. 使用核主成分分析进行非线性映射

许多机器学习算法都嘉定输入数据是线性可分的。感知器为了保证其收敛性，甚至要求悬链数据是完美线性可分的。我们目前学习过的算法中，像 Adaline、逻辑斯特回归和(标准)支持向量机等，都无法实现完美的线性划分的原因归咎于噪声。然而，在现实世界中，大多数情况下我们面对的是非线性问题，针对此类问题，通过降维技术，如 PCA 和 LDA 等，将其转化为线性问题并不是最好的方法。在本小节中，我们将了解一下利用核技巧的 PCA，或者称为核 PCA，这与在第3章中我们介绍过的核支持向量机的概念有一定相关性。使用核 PCA，我们将学习如何将非线性可分的数据转换到一个适合对其进行线性分类的新的低维子空间上。

3.1 核函数与核技巧

会议一下我们再第三章中曾讨论的基于核的支持向量机，通过将非线性可分问题映射到维度更高的特征空间，使其在新的特征空间上线性可分。为了将样本 \(x \in R^d\) 转换到维度更高的 \(k\) 维子空间，我们定义如下非线性映射函数 \(\varphi\):

\[ \phi: R^d \rightarrow R^k (k >> d) \]

我们可以将 \(\varphi\)看作是一个函数，它能够对原始特征进行非线性组合，以将原始的 \(d\) 维数据集映射到更高维的 \(k\) 维特征空间。例如，对于二维(d=2) 特征向量\(x \in R^d\)(\(x\) 是包含 \(d\) 个特征的列向量)来说，可用如下映射将其转换到三维空间:

\[ x = \begin{bmatrix} x_1, x_2 \end{bmatrix}^T \\ \downarrow \phi \\ z = \begin{bmatrix} x_1^2, \sqrt {2x_1x_2}, x_2^2 \end{bmatrix}^T \]

话句话说，利用核 PCA，我们可以通过非线性映射将数据转换到一个高维空间，然后在此高维空间中使用标准 PCA 将其映射到另外一个低维空间中，并通过线性分类器对样本进行划分(前提条件是，样本可以根据输入空间的密度进行划分)。但是，这种方法的缺点是会带来高昂的计算成本，这也正是我们为什么要使用核技巧的原因。通过使用核技巧，我们可以在原始特征空间中计算两个高维特征空间中向量的相似度。

在更深入了解了使用核技巧解决计算成本高昂的问题之前，我们先回顾一下本章最初锁介绍的标准 PCA 方法。两个特征 \(k\) 和 \(j\) 之间协方差的计算公式如下:

\[ \sigma_{jk} = \frac {1}{n} \sum_{i=1}^n (x_j^{(i)} - \mu_j)(x_k^{(i)} - \mu_k) \]

由于在对特征做标准化处理后，其均值为0。例如: 我们可将上述公式简化为:

\[ \sigma_{jk} = \frac {1}{n} \sum_{i=1}^n x_j^{(i)} x_k^{(i)} \]

请注意，上述公式是两个特征值之间的协方差计算公式，下面给出计算协方差矩阵 \(\sum\) 的通用公式:

\[ \Sigma = \frac {1}{n} \sum_{i=1}^n x^{(i)} x^{(i)T} \]

Bernhard Schoellkopf 提出了一种方法\(^{注1}\)，可以使用 \(\varphi\) 通过在原始特征空间上的非线性特征组合来替代样本间点积的计算:

注1: B.Scholkopf, A. Smola, and K.-R. Muler. Kernel Principal Component Analysis. pages 583-588, 1997.

\[ \Sigma = \frac {1}{n} \sum_{i=1}^n \phi (x^{(i)}) \phi (x^{(i)})^T \]

为了求得此协方差矩阵的特征向量，也就是主成分，我们需要求解下述公式:

\[ \Sigma \nu = \lambda \nu \\ \Rightarrow \frac {1}{n} \sum_{i=1}^n \phi (x^{(i)}) \phi (x^{(i)})^T \nu = \lambda \nu \\ \Rightarrow \nu = \frac {1}{\lambda n} \sum_{i=1}^n \phi (x^{(i)}) \phi (x^{(i)})^T \nu = \lambda \nu = \frac {1}{n} \sum_{i=1}^n a^{(i)} \phi (x^{(i)}) \]

其中，\(\lambda\) 和 \(\nu\) 分别为协方差矩阵 \(\Sigma\) 的特征值和特征向量，这里的 \(a\) 可以通过提取核(相似) 矩阵 \(K\) 的特征向量来得到，具体内容在将后续段落中进行介绍。

首先，使用矩阵符号来表示协方差矩阵，其中 \(\varphi (X)\) 是一个 \(n \times k\)维的矩阵:

\[ \Sigma = \frac {1}{n} \sum_{i=1}^n \phi (x^{(i)}) \phi (x^{(i)})^T = \frac {1}{n} \phi (X)^T \phi (X) \]

现在，我们可以将特征向量的公式记为:

\[ \nu = \frac {1}{n} \sum_{i=1}^n a^{(i)} \phi (x^{(i)}) = \nu \phi(X)^T a \]

由于 \(\Sigma \nu = \lambda \nu\)，我们可以得到:

\[ \frac {1}{n} \phi(X)^T \phi(X) \phi(X)^T a = \lambda \phi(X) a \]

两边同乘以 \(\varphi(X)\)，可得:

\[ \frac {1}{n} \phi(X) \phi(X)^T \phi(X) \phi(X)^T a = \lambda \phi(X) \phi(X)^T a \\ \Rightarrow \phi(X) \phi(X)^T a = \lambda a \\ \Rightarrow \frac {1}{n} K a = \lambda a \]

其中，\(K\) 为相似(核) 矩阵:

\[ K = \phi(X) \phi(X)^T \]

回顾一下 3.4节的内容，通过核技巧，使用核函数 \(K\) 以避免使用 \(\varphi\) 来精确计算样本集合 \(x\) 中样本对之间的点积，这样我们就无需对特征向量进行精确的计算:

\[ K(x^{(i)}, x^{(j)}) = \phi(x^{(i)})^T \phi(x^{(j)}) \]

换句话说，通过核 PCA，我们能够得到已经映射到各成分的样本，而不像标准 PCA 方法那样去构建一个转换矩阵。简单地说，可以将核函数(或者简称核) 理解为: 通过两个向量点积来度量向量间相似度的函数。最常用的核函数有:

多项式核:

\[ K(x^{(i)}, x^{(j)}) = (x^{(i)T} x^{(j)} + \theta)^p \]

其中，阈值 \(\theta\) 和幂的值 \(p\) 需自行定义。

双曲正切(sigmoid)核:

\[ K(x^{(i)}, x^{(j)}) = thah(\eta x^{(i)T} x^{(j)} + \theta) \]

径向基核函数(Radial Basis Function, RBF) 或者成为高斯核函数，我们将在下一小节的示例中用到:

\[ K(x^{(i)}, x^{(j)}) = exp \large( -\frac {||x^{(i)} - x^{(j)}||^2}{2 \sigma^2} \large ) \]

也可以写作:

\[ K(x^{(i)}, x^{(j)}) = exp (-\gamma {||x^{(i)} - x^{(j)}||^2}) \]

综合上述讨论，我们可以通过如下三个步骤来实现一个基于 RBF 核的 PCA:

1) 为了计算核(相似)矩阵 \(k\)，我们需要做如下计算:

\[ K(x^{(i)}, x^{(j)}) = exp (-\gamma {||x^{(i)} - x^{(j)}||^2}) \]

我们需要计算任意两样本对之间的值:

\[ K = \begin{bmatrix} K(x^{(1)}, x^{(1)}) & K(x^{(1)}, x^{(2)}) & \cdots & K(x^{(1)}, x^{(n)}) \\ K(x^{(2)}, x^{(1)}) & K(x^{(2)}, x^{(2)}) & \cdots & K(x^{(2)}, x^{(n)}) \\ \vdots & \vdots & \ddots & \vdots \\ K(x^{(n)}, x^{(1)}) & K(x^{(n)}, x^{(2)}) & \cdots & K(x^{(n)}, x^{(n)}) \\ \end{bmatrix} \]

例如: 如果数据集包含 100 个训练样本，将得到一个 \(100 \times 100\) 维的对称矩阵。

2) 通过如下公式进行计算，使核矩阵 \(K\) 更为聚集:

\[ K' = K-1_nK - K1_n + 1_nK1_n \]

其中，\(l_n\) 是一个 \(n \times n\) 维的矩阵(与核矩阵维度相同)，其所有的值均为 \(\frac {1}{n}\)。

3) 将聚集后的核矩阵的特征值按照降序排列，选择前 \(k\) 个特征值所对应的特征向量。与标准 PCA 不同，这里的特征向量不是主成分轴，而是将样本映射到这些轴上。

至此，读者可能会感到困惑: 为什么要在第2步对核矩阵进行聚集处理? 我们曾经假定，数据需要经过标准化处理，当在生成协方差矩阵并通过 \(\varphi\) 以非线性特征的组合替代点积时，所有特征的矩阵为0。由此，在第2部轴聚集核矩阵就显得很有必要，因为我们并没有精确计算新的特征空间，而且也不能确定新特征空间的中心在零点。

下一小节中，我们将基于这三个步骤使用 Python 实现 PCA。

3.2 使用 Python 实现核主成分分析

在千米你的小节中，我们讨论了核 PCA 相关的核心概念。现在根据之前总结过的实现核 PCA 方法的三个步骤，使用 Python 实现基于 RBF 核的 PCA。借助于 SciPy 和 NumPy 的函数，我们将会看到，实现 PCA 实际上很简单。

from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigen
import numpy as np


def rbf_kernel_pca(X, gamma, n_components):
    """
    RBF kernel PCA implementation.

    Parameters
    -----------
    X: {NumPy ndarray}, shape = [n_samples, n_features]

    gamma: float
        Turning parameter of the RBF kernel

    n_components: float
        Number of principal components to return

    Returns
    -----------
    X_pc: {NumPy ndarray}, shape = [n_samples, k_features]
        Projected dataset
    """
    # Calculate pairwise squared Euclidean distances
    # in the MxN dimensional dataset.
    sq_dists = pdist(X, 'sqeuclidean')

    # Convert pairwise distance into a square matrix
    mat_sq_dists = squreform(sq_dists)

    # Compute the symmetric kernel matrix
    K = exp(-gamma * mat_sq_dists)

    # Center the kernel matrix
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K)  - K.dot(one_n) + one_n.dot(K).dot(one_n)

    # Obtaining eigenpairs from the centered kernel matrix
    # numpy.eigh returns them in sorted order
    eigvals, eigvecs = eigh(K)

    # Collect the top k eigenvectors (projected samples)
    X_pc = np.column_stack((eigvecs[:, i] for i in range(1, n_components + 1)))
    return X_pc

采用 RBF 核函数实现的 PCA 进行降维时存在一个问题，就是我们必须制定先验证参数 r，需要通过实验来找到一个合适的 r 值，最好是通过参数调优来确定，例如网格搜索算法，我们将在第6章中对其进行深入探讨。

示例1: 分离半月形数据

现在，我们将实现 rbf_kernel_pca 方法应用于非线性示例数据集。我们首先创建一个包含 100 个样本点的二维数据集，以两个半月形状表示:

from sklearn.datasets import make_moons
X, y = make_moons(n_samples=100, random_state=123)
plt.scatter(X[y==0, 0], X[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y==1, 0], X[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.show()

处于演示的需要，使用三角符号标识的标识一个类别中的样本，使用圆形符号标识的表示另一类别的样本:

显然，这两个半月形不是线性可分的，而我们的目标是通过核 PCA 将这两个半月形数据展开，使得数据成为适用于某一线性分类器的输入数据。首先，我们通过标准的 PCA 将数据映射到主成分上，并观察其形状:

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

scikit_pca = PCA(n_components=2)
X_spca = scikit_pca.fit_transform(X)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))
ax[0].scatter(X_spca[y==0, 0], X_spca[y==0, 1],
             color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y == 1, 0], X_spca[y == 1, 1],
              color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y == 0, 0], np.zeros((50, 1)) + 0.02,
              color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y == 1, 0], np.zeros((50, 1)) - 0.02,
              color='blue', marker='o', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.tight_layout()
plt.show()           

很明显，经过标准 PCA 的转换后，线性分类器未必能很好地发挥作用:

请注意，当我们仅回执第一主成分的图像时(见右子图)，我们分别将三角形和圆形代表的样本向上或向下做了轻微调整，以更好地展示类间重叠。

🔖 PCA 是无监督方法，与 LDA 相比，它在使得方差最大化的过程中未使用类标信息。出于增强可视化效果的考虑，为了显示分类的程度，我们才在此使用了三角形和圆形符号。」

现在我们将使用前一小节中实现的和核 PCA 函数 rbf_kernel_pca:

from matplotlib.ticker import FormatStrFormatter

X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)

fig, ax = plt.subplots(nrows=1,ncols=2, figsize=(7,3))
ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], 
            color='red', marker='^', alpha=0.5)
ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1],
            color='blue', marker='o', alpha=0.5)

ax[1].scatter(X_kpca[y==0, 0], np.zeros((50,1))+0.02, 
            color='red', marker='^', alpha=0.5)
ax[1].scatter(X_kpca[y==1, 0], np.zeros((50,1))-0.02,
            color='blue', marker='o', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
ax[0].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))

plt.tight_layout()
plt.show()

可以看到，两个类别(圆形和三角形) 此时是线性可分的，这使得转换后的数据适合作为线性分类器的训练数据集:

不过，对于可调整参数 \(\gamma\)，没有一个通用的值使其适用于不同的数据集。针对给定问题找到一个适宜的参数值需要通过实验来解决。在第6章中，我们将讨论可自动进行参数优化等任务的技术。这里，我使用了一个已有的能够生成良好结果的值。

示例2. 分离同心圆

在上一小节，我们演示了如何通过核 PCA 分离半月形数据。既然我们已经投入了如此多的经历去理解核 PCA 的概念，就再看一下另外一个关于非线性问题的有趣例子——同心圆。

代码如下:

from sklearn.datasets import make_circles
X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)

plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue', marker='o', alpha=0.5)
plt.tight_layout()
plt.show()

同样，我们假设了一个涉及两个类别的问题，三角形和圆形分别标识不同类别的样本:

首先使用标准 PCA 方法，以便将结果与基于 RBF 核的 PCA 生成的结果进行比较:

scikit_pca = PCA(n_components=2)
X_spca = scikit_pca.fit_transform(X)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))

ax[0].scatter(X_spca[y == 0, 0], X_spca[y == 0, 1],
              color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y == 1, 0], X_spca[y == 1, 1],
              color='blue', marker='o', alpha=0.5)

ax[1].scatter(X_spca[y == 0, 0], np.zeros((500, 1)) + 0.02,
              color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y == 1, 0], np.zeros((500, 1)) - 0.02,
              color='blue', marker='o', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')

plt.tight_layout()
plt.show()

再一次发现，通过标准 PCA 无法得到适合于线性分类器的训练数据:

给定一个合适的 \(\gamma\) 值，来看看基于 RBF 的核 PCA 实现能够得到令人满意的结果:

X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
ax[0].scatter(X_kpca[y == 0, 0], X_kpca[y == 0, 1],
              color='red', marker='^', alpha=0.5)
ax[0].scatter(X_kpca[y == 1, 0], X_kpca[y == 1, 1],
              color='blue', marker='o', alpha=0.5)

ax[1].scatter(X_kpca[y == 0, 0], np.zeros((500, 1)) + 0.02,
              color='red', marker='^', alpha=0.5)
ax[1].scatter(X_kpca[y == 1, 0], np.zeros((500, 1)) - 0.02,
              color='blue', marker='o', alpha=0.5)

ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')

plt.tight_layout()
plt.show()

基于 RBF 的核 PCA 再一次将数据映射到了一个新的子空间中，使两个类别变得线性可分:

3.3 映射新的数据点

在上述两个核 PCA 应用的例子(半月形和同心圆)中，我们都将单一数据集映射到一个新的特征上。但在实际应用中，我们可能需要转换多个数据集，例如，训练数据、测试数据，以及在完成模型构建和评估后所要收集的新样本。本节将介绍如何映射训练数据集以外的数据点。

在本章开始时介绍过的标准 PCA 方法中，我们通过转换矩阵和输入样本之间的点积来对数据进行映射；映射矩阵的列是协方差矩阵中 \(k\) 个最大特征值所对应的特征向量(v)。现在的问题是: 如此将此概念应用于核 PCA？回忆一下核 PCA 的原理可以记得，我们从聚集核矩阵(不是协方差矩阵) 中得到了特征向量(a)。这一位这样本已经映射到了主成分轴v。由此，如果我们希望将新的样本(\(x'\)) 映射到此主成分轴，需要进行如下计算:

\[\phi(x')^T \nu\]

幸运的是，我们可以使用核技巧，这样就无需精确计算映射 \(\varphi(x')^T \nu\)。然而值得注意的是: 与标准 PCA 相比，核 PCA 是一种基于内存的方法，这意味着每次映射新的样本前，必须再次使用训练数据。我们需要计算训练数据集中每一个训练样本和新样本 \(x'\) 之间的 RBF 核(相似度):

\[ \begin{align} \phi(x')^T \nu &= \sum_i a^{(i)} \phi(x')^T \phi(x^{(i)}) \\ & = \sum_i a^{(i)} k(x', x^{(i)})^T \end{align} \]

其中，核矩阵 K 的特征向量 a 及特征值 \(\lambda\) 需满足如下等式:

\[ Ka = \lambda a \]

在完成新样本与训练数据集内样本间相似度的计算后，我们还需通过特征向量对应的特征值来对其进行归一化处理。可以通过修改早前实现过的 rbf_kernel_pca 函数来让其返回矩阵的特征值:

from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):
    """
    RBF kernel PCA implementation.

    Parameters
    ------------
    X: {NumPy ndarray}, shape = [n_samples, n_features]

    gamma: float
      Tuning parameter of the RBF kernel

    n_components: int
      Number of principal components to return

    Returns
    ------------
     X_pc: {NumPy ndarray}, shape = [n_samples, k_features]
       Projected dataset   

     lambdas: list
       Eigenvalues

    """
    # Calculate pairwise squared Euclidean distances
    # in the MxN dimensional dataset.
    sq_dists = pdist(X, 'sqeuclidean')

    # Convert pairwise distances into a square matrix.
    mat_sq_dists = squareform(sq_dists)

    # Compute the symmetric kernel matrix.
    K = exp(-gamma * mat_sq_dists)

    # Center the kernel matrix.
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

    # Obtaining eigenpairs from the centered kernel matrix
    # numpy.eigh returns them in sorted order
    eigvals, eigvecs = eigh(K)

    # Collect the top k eigenvectors (projected samples)
    alphas = np.column_stack((eigvecs[:, -i]
                              for i in range(1, n_components + 1)))

    # Collect the corresponding eigenvalues
    lambdas = [eigvals[-i] for i in range(1, n_components + 1)]

    return alphas, lambdas

至此，我们可以创建一个新的半月形数据集，并使用更新过的 RBF 核PCA 实现来将其映射到一个一维的空间上:

X, y = make_moons(n_samples=100, random_state=123)
alphas, lambdas = rbf_kernel_pca(X, gamma=15, n_components=1)

为了确保我们已经完成了实现新样本映射的代码，嘉定半月形数据集中的第26个点事一个新的数据点 \(x'\)，现在要将其映射到新的子空间中:

>>> x_new = X[-1]
>>> x_new
array([ 0.4816, -0.3551])

>>> x_proj = alphas[-1] # original projection
>>> x_proj
array([ 0.1192])

通过执行下面的代码，我们可以重现原始映射。使用 project_x 函数，还可以映射新的数据样本。代码如下:

def project_x(x_new, X, gamma, alphas, lambdas):
    pair_dist = np.array([np.sum((x_new - row)**2) for row in X])
    k = np.exp(-gamma * pair_dist)
    return k.dot(alphas / lambdas)

# projection of the "new" datapoint
x_reproj = project_x(x_new, X, gamma=15, alphas=alphas, lambdas=lambdas)
x_reproj 

# 输出:
array([ 0.1192])

最后，将第一主成分上的映射进行可视化:

plt.scatter(alphas[y == 0, 0], np.zeros((50)),
            color='red', marker='^', alpha=0.5)
plt.scatter(alphas[y == 1, 0], np.zeros((50)),
            color='blue', marker='o', alpha=0.5)
plt.scatter(x_proj, 0, color='black',
            label='original projection of point X[25]', marker='^', s=100)
plt.scatter(x_reproj, 0, color='green',
            label='remapped point X[25]', marker='x', s=500)
plt.legend(scatterpoints=1)

plt.tight_layout()
plt.show()

从下图可见，我们将样本 \(x'\) 正确映射到了第一主成分上:

X, y = make_moons(n_samples=100, random_state=123)
alphas, lambdas = rbf_kernel_pca(X[:-1, :], gamma=15, n_components=1)

def project_x(x_new, X, gamma, alphas, lambdas):
    pair_dist = np.array([np.sum((x_new - row)**2) for row in X])
    k = np.exp(-gamma * pair_dist)
    return k.dot(alphas / lambdas)

# projection of the "new" datapoint
x_new = X[-1]

x_reproj = project_x(x_new, X[:-1], gamma=15, alphas=alphas, lambdas=lambdas)

plt.scatter(alphas[y[:-1] == 0, 0], np.zeros((50)),
            color='red', marker='^', alpha=0.5)
plt.scatter(alphas[y[:-1] == 1, 0], np.zeros((49)),
            color='blue', marker='o', alpha=0.5)
plt.scatter(x_reproj, 0, color='green',
            label='new point [ 100.0,  100.0]', marker='x', s=500)
plt.legend(scatterpoints=1)

plt.scatter(alphas[y[:-1] == 0, 0], np.zeros((50)),
            color='red', marker='^', alpha=0.5)
plt.scatter(alphas[y[:-1] == 1, 0], np.zeros((49)),
            color='blue', marker='o', alpha=0.5)
plt.scatter(x_proj, 0, color='black',
            label='some point [1.8713,  0.0093]', marker='^', s=100)
plt.scatter(x_reproj, 0, color='green',
            label='new point [ 100.0,  100.0]', marker='x', s=500)
plt.legend(scatterpoints=1)

plt.tight_layout()
plt.show()

5.3.4 scikit-learn 中的核主成分分析

scikit-learn 的 sklearn.decomposition 子模块中已经实现了一种核 PCA 类。其使用方法与标准 PCA 类类似，我们可以通过 kernel 参数来选择不同的核函数:

from sklearn.decomposition import KernelPCA

X, y = make_moons(n_samples=100, random_state=123)
scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_skernpca = scikit_kpca.fit_transform(X)

为了验证得到的结果与我们自己实现的核 PCA 是否一致，我们来回执半月形数据映射到前两个主成分的图像:

plt.scatter(X_skernpca[y == 0, 0], X_skernpca[y == 0, 1],
            color='red', marker='^', alpha=0.5)
plt.scatter(X_skernpca[y == 1, 0], X_skernpca[y == 1, 1],
            color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')
plt.ylabel('PC2')
plt.tight_layout()
plt.show()

从结果图像中可见，通过 scikit-learn 中 KernelPCA 得到的结果与我们自己实现的结果一致:

🔖 scikit-learn 实现了一些高级的非线性降维技术，这些内容已经超出了本书的范围。读者可以通过链接 http://scikit-learn.org/stable/modules.mainfold.html 来了解相关内容概述及其示例。

4. 本章小结

在本章中，读者学习了三种不同的基于特征提取的基本降维技术: 标准 PCA、LDA，以及核 PCA。使用 PCA，我们可以在忽略类标的情况下，将数据映射到一个低维的子空间上，并沿正交的特征坐标方向使方差最大化。与 PCA 不同，LDA 是一种监督降维技术，这意味着: 在线性特征中尝试类别最大可分时，需要使用训练数据集中的类别信息。最后，我们学习了核 PCA，通过核 PCA 可以将非线性数据集映射到一个低维的特征空间中，使得数据线性可分。

在掌握了这些数据预处理的基本技术之后，读者可以准备进入下一章，学习如何有效地组合不同的预处理技术，以及模型性能评估等方面的内容。