在第4章 构建优秀的训练数据集 – 数据预处理中，我们学习了使用特征选择技术对数据集降维的不同方法。特征选择以外的另一种降维方法是特征提取。本章中我们会学习两种基本技术，可帮助我们通过将其变换为比原来更低维度的特征子空间总结出数据集中的信息内容。数据压缩是机器学习中非常重要的课题，它有助于我们存储和分析现代技术时代生产和收集的与日俱增的数据。
本章中我们会讲解如下内容：
- 用于无监督数据压缩的主成分分析
- 最大化类别分割性监督降维的线性判别分析
- 非线性降维技术的概览及用于数据可视化的t-分布随机近邻嵌入
类似于特征选择，我们可以使用不同的特征提取技术来减少数据集中的特征数量。特征选择与特征提取的不同之处在于使用特征选择算法时我们保留原始特征，如序列后向选择，我们使用特征提取变换或投射数据到新的特征空间上。
在进行降维时，特征提取可以理解为一种数据压缩技术，目标是保留大部分相关信息。实操时，特征提取不仅用于改善存储空间或学习算法的计算效率，还通过降低维数灾难来改善预测性能，在处理非正则化模型时尤其如此。
本节中我们会讨论主成分分析(PCA)，一种广泛用于不同领域的无监督线性变换技术，对特征提取和降维尤为突出。其它知名的PCA应用有股票市场交易中的探索性数据分析及信号去噪，以及生物信息领域的基因数据分析和基因表达水平。
PCA帮助我们根据特征间关联识别数据中的模式。总之，PCA旨在找到高维数据中的最大方差方向并将数据投射到新的子空间上，维数等于或小于原空间。新子空间的正交轴（主成分）可解释为限定新特征轴彼此正交时最大方差的方向，如图5.1所示：
图5.1：使用PCA查找数据集中的最大方差方向
在图5.1中，x1和x2是原始特征轴，PC 1和PC 2是主成分。
如果使用PCA降维，我们会构建一个d×k维的变换矩阵W，可以将训练样本特征的向量x映射到一个新的k-维特征子空间，其维数少于原始的d-维特征空间。例如下面的流程。假设我们有一个特征向量x：
然后通过变换矩阵进行变换，$W\in \mathbb{R}^{d\times k}$：
xW = z
产生输出向量：
将原始d-维数据转换到新的k-维子空间(通常k << d)上的结果是，第一个主成分会拥有最大的方差。后续的主成分在这些成分与其它主成分不相关（正交）时会具有最大方差，即例输入特征相关联，得到的主成分也互为正交（不相关）。注意PCA的方向对数据缩放超级敏感，如果特征使用不同量级度量而我们又希望为所有特征赋相同的重要性，那么应需要在做PCA之前标准化特征。
在更进一步学习用PCA算法实现降维前，我们先将该方法总结成简单的步骤：
- 标准化d-维数据集。
- 构建协方差矩阵。
- 将协方差分解为特征向量（eigenvectors）和特征值（eigenvalues）。
- 通过特征值降序排列相应的特征向量。
- 选取k最大特征值的k特征向量，其中k是新特征字空间的维数（$k\le d$）。
- 通过“头部”k特征向量构建投影矩阵W。
- 使用投影矩阵W变换d-维输入数据集X，获取一新的 k-维特征子空间。
在下面的小节中，作为练习我们会使用Python逐步执行PCA。然后我们会学习如何使用scikit-learn更便利地执行PCA。
特征分解：将矩阵分解为特征向量和特征值
特征分解，将方阵分解为特征值和特征向量，是本节中所述的PCA处理的信心。
协方差矩阵是方阵的一种特例：它是一个对称矩阵，也就是矩阵与其转置相等，A = AT。
在分解这种对称矩阵时，特征值是实数（而非复数），特征向量彼此正交（垂直相交）。此外，特征值与特征向量成对出现。如果将协方差矩阵分解为特征向量和特征值，与最高特征值相关联的特征向量对数据集中最大方差的方向。这里的“方向”是数据集特征列的线性变换。
有关特征值和特征向量更深入的讨论不在本书范畴内，比较详尽的参考资料请见维基百科：https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors。
在本小节中，我们会处理PCA的前四个步骤：
- 标准化数据
- 构建协方差矩阵
- 获取协方差矩阵的特征值和特征向量
- 以特征值降序排列特征向量
首先我们会加载第4章 构建优秀的训练数据集 - 数据预处理中使用的葡萄酒数据集：
>>> import pandas as pd
>>> df_wine = pd.read_csv(
... 'https://archive.ics.uci.edu/ml/'
... 'machine-learning-databases/wine/wine.data',
... header=None
... )
获取葡萄酒数据集
可以在本书代码库中找到一份葡萄酒数据集（以及本书中使用的其它数据集）的拷贝，以防读者离线操作或是UCI服务器上https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data临时断网。比如我们从本地目录读者葡萄酒数据集，只需将如下行：
df = pd.read_csv( 'https://archive.ics.uci.edu/ml/' 'machine-learning-databases/wine/wine.data', header=None )替换为：
df = pd.read_csv( 'your/local/path/to/wine.data', header=None )
接下来我们会将葡萄酒数据分成训练集和测试集，分别占70%和30%，并将其标准化为单位方差：
>>> from sklearn.model_selection import train_test_split
>>> X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
>>> X_train, X_test, y_train, y_test = \
... train_test_split(X, y, test_size=0.3,
... stratify=y,
... random_state=0)
>>> # standardize the features
>>> from sklearn.preprocessing import StandardScaler
>>> sc = StandardScaler()
>>> X_train_std = sc.fit_transform(X_train)
>>> X_test_std = sc.transform(X_test)
通过以上代码完成自主预处理后，我们进入第2步：构建协方差矩阵。对称的d×d-维协方差矩阵，其中d是数据集中的维数，存储着不同特征间的成对协方差。例如，特征xj和xk之间对总体的协方差可通过如下等式计算：
这里的$\mu_j$和$\mu_k$分别是样本j和k的均值。注意如果我们标准化了数据集那么样本均值就是0。两个特征间的协方差为正值表示特征同时上升或下降，或负协方差则表示特征的方向相反。例如，三个特征的协方差矩阵可以写成如下这样（注意这里的$\Sigma$是大写希腊字母，不要与加和符号相混淆）：
协方差矩阵的特征向量表示主成分（最大方差的方向），而对应的特征值则定义它们的量级。就葡萄酒一例来说，我们获取了12个特征向量以及13×13-维协方差矩阵中的值。
下面到第3步，我们来获取协方差矩阵的特征对。如果读者学过线性代码，可能学过满足以下条件的特征向量v：
这里的$\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', eigen_vals)
Eigenvalues
[ 4.84274532 2.41602459 1.54845825 0.96120438 0.84166161
0.6620634 0.51828472 0.34650377 0.3131368 0.10754642
0.21357215 0.15362835 0.1808613 ]
使用numpy.cov函数，我产计算了标准化的训练集的协方差矩阵。使用linalg.eig函数，我们执行了特征分解，这会产生一个包含13个特征值和对应的存储为13×13-维矩阵的特征向量(eigen_vecs)的向量(eigen_vals)。
NumPy中的特征分解
numpy.linalg.eig设计用于操作对称及非对称方阵。但在有些情况下会发现它返回的是复数特征值。相关的函数
numpy.linalg.eigh，是为分解共轭（Hermetian）矩阵而实现的，对于处理协方差矩阵这样的对称矩阵在数值上更为稳定，numpy.linalg.eigh一定会返回实数特征值。
因为我们希望通过将数据集压缩到新的特征子空间上实现降维，所以只选择了包含大部分信息（方差）的特征向量子集（主成分）。特征值定义了特征向量的量级，所以我们按降序对特征值排序，我们感兴趣的是相应特征值得到的k个最高特征向量。但在采集这k个信息最丰富的特征向量前，我们先绘制特征值的方差解释率。特征值$\lambda_j$的方差解释率，只需使用特征值$\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, align='center',
... label='Individual explained variance')
>>> plt.step(range(1,14), cum_var_exp, where='mid',
... label='Cumulative explained variance')
>>> plt.ylabel('Explained variance ratio')
>>> plt.xlabel('Principal component index')
>>> plt.legend(loc='best')
>>> plt.tight_layout()
>>> plt.show()
生成的图表明仅第一个主成分就占了方差的大约40%。
同时，我们可以看出前两个主成分合起来用解释了数据集中差不多60%的方差：
图5.2：通过主成分捕获的总方差比例
虽然可解释方差图让我们想起了第4章 构建优秀的训练数据集 - 数据预处理中随机森林的特征重要性值，但读者要记住PCA是一种无监督方法，表示它会忽略类标签的信息。而随机森林使用类成员信息来计算节点杂度，方差度量沿特征轴的扩展程度。
我们已经成功将协方差矩阵分解为了特征对，下面继续将葡萄酒数据集变换为主成分轴的最后3个步骤。本节中将处理的剩余步骤为：
- 选取k最大特征值的k特征向量，其中k是新特征子空间的维数（$k\le d$）。
- 通过“头部”k个特征向量构建投影矩阵W。
- 使用投影矩阵W变换d-维输入数据集X，获取一个新的k-维特征子空间。
换句不那么技术的话来说，我们会按特征值降序排列特征对，通过选中的特征向量构建投影矩阵，并使用投影矩阵将数据变换到更低维的子空间上。
我们先通过特征值降序排列特征对：
>>> # Make a list of (eigenvalue, eigenvector) tuples
>>> eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
... for i in range(len(eigen_vals))]
>>> # Sort the (eigenvalue, eigenvector) tuples from high to low
>>> eigen_pairs.sort(key=lambda k: k[0], 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.13724218 0.50303478]
[ 0.24724326 0.16487119]
[-0.02545159 0.24456476]
[ 0.20694508 -0.11352904]
[-0.15436582 0.28974518]
[-0.39376952 0.05080104]
[-0.41735106 -0.02287338]
[ 0.30572896 0.09048885]
[-0.30668347 0.00835233]
[ 0.07554066 0.54977581]
[-0.32613263 -0.20716433]
[-0.36861022 -0.24902536]
[-0.29669651 0.38022942]]
执行以上代码，我们通过最高的两个特征向量创建了一个13×2维的投影矩阵W。
镜像投影
根据所使用的NumPy和LAPACK的版本，可能会获得符号相反的矩阵W。这并不是什么问题，如果v是矩阵$\Sigma$的特征向量：
这里的v是特征向量，–v也是一个特征向量，可在下面看到。使用基础代数，可以在等式两边都乘上一个标量$\alpha$:
因矩阵乘法对标量具有结合律，可重新排列如下：
现在可以看到在$\alpha=1$和$\alpha=-1$时$\alpha v$是相同特征值$\lambda$的特征向量。因此v和–v 都是特征向量
使用投影矩阵，我们可以将示例x（表现为13-维的行向量）变换到PCA子空间（主成分1和2）上获取 x′，这时它是一个包含两个新特征的二维示例向量：
x′ = xW
>>> X_train_std[0].dot(w)
array([ 2.38299011, 0.45458499])
类似地，我们可以通过计算矩阵点乘把整个124×13-维训练集变换成两个主成分：
X′ = XW
>>> X_train_pca = X_train_std.dot(w)
最后，我们来可视化变换后的葡萄酒训练集，现在以一个二维散点图存储着一个124×2-维矩阵：
>>> colors = ['r', 'b', 'g']
>>> markers = ['o', 's', '^']
>>> 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=f'Class {l}', marker=m)
>>> plt.xlabel('PC 1')
>>> plt.ylabel('PC 2')
>>> plt.legend(loc='lower left')
>>> plt.tight_layout()
>>> plt.show()
由图5.3可以看出，数据沿第一个主成分(x轴)要比第二个主成分(y轴)散得更开，这与前面小节中创建的可解释方差率图是一致的。但我们能看出一个线性分类器很可能可以很好地进行分类：
图5.3：通过PCA将葡萄酒数据记录投射到2D特征空间
虽然我们编码类标签信息的目的是绘制上面的散点图，但要谨记PCA是一种不使用任何类标签信息的无监督技术。
虽然上一小节中繁琐的介绍有助于我们了解PCA的内部原理，但我们还是要谈谈scikit-learn中实现的PCA类。
PCA类是scikit-learn中又一个变换器类，我们在使用相同模型参数谈的训练数据和测试集前先使用训练数据拟合模型。下面对葡萄酒训练数据集使用scikit-learn中的PCA类，通过逻辑回归分类变换后的样本，再使用第2章 为分类训练简单机器学习算法中定义的plot_decision_regions函数可视化决策区域：
from matplotlib.colors import ListedColormap
def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):
# setup marker generator and color map
markers = ('o', 's', '^', '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))
lab = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
lab = lab.reshape(xx1.shape)
plt.contourf(xx1, xx2, lab, alpha=0.3, cmap=cmap)
plt.xlim(xx1.min(), xx1.max())
plt.ylim(xx2.min(), xx2.max())
# plot class examples
for idx, cl in enumerate(np.unique(y)):
plt.scatter(x=X[y == cl, 0],
y=X[y == cl, 1],
alpha=0.8,
c=colors[idx],
marker=markers[idx],
label=f'Class {cl}',
edgecolor='black')
为方便起见，可以将前面的plot_decision_regions代码放到当前工作目录中单独的文件内，例如，plot_decision_regions_script.py，然后在当前Python会话中导入它：
>>> from sklearn.linear_model import LogisticRegression
>>> from sklearn.decomposition import PCA
>>> # initializing the PCA transformer and
>>> # logistic regression estimator:
>>> pca = PCA(n_components=2)
>>> lr = LogisticRegression(multi_class='ovr',
... random_state=1,
... solver='lbfgs')
>>> # dimensionality reduction:
>>> X_train_pca = pca.fit_transform(X_train_std)
>>> X_test_pca = pca.transform(X_test_std)
>>> # fitting the logistic regression model on the reduced dataset:
>>> lr.fit(X_train_pca, y_train)
>>> plot_decision_regions(X_train_pca, y_train, classifier=lr)
>>> plt.xlabel('PC 1')
>>> plt.ylabel('PC 2')
>>> plt.legend(loc='lower left')
>>> plt.tight_layout()
>>> plt.show()
通过执行这段代码，应该会看训练数据的决策区域降低到两个主成分轴：
图5.4：使用scikit-learn的PCA降维后的训练样本和逻辑回归决策区域
对比scikit-learn的PCA投射和我们自己实现的PCA，你可能会发现产生的图彼此互为镜像。请注意这不是说两种实现中有什么错误，出现不同的原因是根据特征值求解的不同，特征向量可能会为正值或负值。
到不是说这很重要，但我们可以轻松通过将数据乘以-1来翻转镜像图片，特征向量一般会缩放到单位长度1。为完整起见，我们对变换的测试数据集绘制逻辑回归决策区域，查看其分类状况：
>>> plot_decision_regions(X_test_pca, y_test, classifier=lr)
>>> plt.xlabel('PC 1')
>>> plt.ylabel('PC 2')
>>> plt.legend(loc='lower left')
>>> plt.tight_layout()
>>> plt.show()
执行如上代码绘制好测试集的决策区域后，可以看到逻辑回归对这个小型二维特征子空间的执行效果很好，只对测试集中很少的几个样本分类有误：
图5.5：通过基于PCA的特征子空间中的逻辑回归决策区域测试数据点
如果对不同主成分的可解释方差比感兴趣的话，可以通过将参数n_components设置为None来初始化PCA，这样会保留所有主成分，然后可通过explained_variance_ratio_属性访问可解释方差比：
>>> pca = PCA(n_components=None)
>>> X_train_pca = pca.fit_transform(X_train_std)
>>> pca.explained_variance_ratio_
array([ 0.36951469, 0.18434927, 0.11815159, 0.07334252,
0.06422108, 0.05051724, 0.03954654, 0.02643918,
0.02389319, 0.01629614, 0.01380021, 0.01172226,
0.00820609])
注意我们在初始化PCA类时设置了n_components=None，这样所有主成分按排序返回，不会执行降维操作。
本节中我们会简要地学习如何评估原始特征对主成分的贡献。我们已经学习过，通过PCA，我们创建了表示特征线性组合的主成分。有时，我们想知道每个原始特征对某一主成分贡献了几成。这些贡献通常称为载荷（loadings）。
因子载荷可通过以特征值平方根缩放特征向量进行计算。生成的值可解释为原始特征与主成分之间的关联。为进行讲解，我们来绘制第一个主成分的载荷。
首先，通过对特征向量乘以特征值平方根来计算13×13-维载荷矩阵。
>>> loadings = eigen_vecs * np.sqrt(eigen_vals)
然后我们绘制第一个主成分的载荷，loadings[:, 0]，它也是矩阵中的第一列：
>>> fig, ax = plt.subplots()
>>> ax.bar(range(13), loadings[:, 0], align='center')
>>> ax.set_ylabel('Loadings for PC 1')
>>> ax.set_xticks(range(13))
>>> ax.set_xticklabels(df_wine.columns[1:], rotation=90)
>>> plt.ylim([-1, 1])
>>> plt.tight_layout()
>>> plt.show()
在图5.6中，可以看到，比如酒精度（Alcohol）与第一个主成分负相关（约–0.3），而苹果酸（ Malic acid）与其正相关（约0.54）。注意值1表示完全正相关，而–1对应着完全负相关：
图5.6：与第一个主成分的特征相关性
在以下的代码示例中，我们对自己的PCA实现计算了因子载荷。可以用类似的方式获取拟合scikit-learn PCA对象的载荷，其中pca.components表示特征向量，pca.explained_variance表示特征值：
>>> sklearn_loadings = pca.components_.T * np.sqrt(pca.explained_variance_)
为比较scikit-learn PCA载荷与些前我们自己所创建的，我们绘制一个类似的柱状图：
>>> fig, ax = plt.subplots()
>>> ax.bar(range(13), sklearn_loadings[:, 0], align='center')
>>> ax.set_ylabel('Loadings for PC 1')
>>> ax.set_xticks(range(13))
>>> ax.set_xticklabels(df_wine.columns[1:], rotation=90)
>>> plt.ylim([-1, 1])
>>> plt.tight_layout()
>>> plt.show()
可以看到这个柱状图基本一样：
图5.7：使用scikit-learn得到的与第一个主成分的特征相关性
在以无监督提取技术探讨了PCA之后，下一节中我们会介绍线性判别分析（linear discriminant analysis (LDA)），这是一种考虑了类标签的线性变换技术。
LDA可用作特征提取技术来提升计算效率和降低因非正则化模型中维度灾难所导致的过拟合程序。LDA背后的基本概念与PCA非常类似，PCA尝试找到数据集中最大方差的正交分量轴，而LDA的目标是找到优化类分离的特征子空间。在下面的小节中，我们会详细讨论LDA和PCA的相似性并逐步实现LDA方法。
PCA和LDA都是可降低数据集中维数的线性变换技术，前者是无监督算法，后者是有监督的。因此，我们可以把LDA看做是比PCA更高级的分类特征提取技术。但A.M. Martinez指出通过PCA对某些图像识别任务做预处理会产生更好的分类结果，比如，每类仅包含很少的样本(PCA Versus LDA by A. M. Martinez and A. C. Kak, IEEE Transactions on Pattern Analysis and Machine Intelligence, 23(2): 228-233, 2001)。
Fisher LDA
LDA有时也称Fisher’s LDA。Ronald A. Fisher首先于1936年为二类分类问题制定了Fisher线性判别（The Use of Multiple Measurements in Taxonomic Problems, R. A. Fisher, Annals of Eugenics, 7(2): 179-188, 1936）。1948年，Fisher线性判别由C. Radhakrishna Rao泛化到多类问题上，前提是相等类协方差和正态分布分类，现在称之为LDA（The Utilization of Multiple Measurements in Problems of Biological Classification by C. R. Rao, Journal of the Royal Statistical Society. Series B (Methodological), 10(2): 159-203, 1948）。
图5.8总结了用于二类问题的LDA的概念。类1的样本使用圆圈表示，类2的样本使用十字符号表示：
图5.8：用于二类问题的LDA的概念
在x轴(LD 1)上所示的线性判别，会很好地分出两个正态分布的类。虽然y轴 (LD 2)上示例线性判别项捕获了数据集中的大量方差，它并不能作为一个很好的线性判别，因为没有捕获到任何类判别信息。
LDA中的一个假设是数据是正态分布的。同时我们假定类具有相同的协方差矩阵，并且训练样本在统计层面上彼此独立。但即使一个或多个假设有（此许的）不符合条件，用于降维的LDA仍可以很好地运行(Pattern Classification 2nd Edition by R. O. Duda, P. E. Hart, and D. G. Stork, New York, 2001)。
在进入代码实现前，我们简要地总结下执行LDA所需的主要步骤：
- 标准化d-维数据集（d为特征数）。
- 为每个类计算d-维平均向量。
- 构建类间散度矩阵SB和类内散度矩阵SW。
- 计算矩阵$S_W^{-1}S_B$的特征向量和对应的特征值。
- 通过特征值降序来排列对应的特征向量。
- 选择对应k最大特征值的k特征向量来构建d×k-维变换矩阵W，特征向量为矩阵的列。
- 使用变换矩阵W将样本投射到新的特征子空间上。
可以看到，LDA与PCA在将矩阵解构为特征值和特征向量层面非常类似，它会形成新的更低维特征空间。但像前面所说的，LDA会考虑类标签信息，这体现在第2步中计算的平均向量。在下一节中，我们会更详细讨论这7个步骤，以及其代码实现。
在本章开始的PCA小节中我们已经对葡萄酒数据集中的特征做过标准化，所以跳过第一步直接计算平均向量，我们会使用它们分别构建类内散度矩阵和类间散度矩阵。每个平均向量mi，存储着类i中样本的平均特征值$\mu_m$：
这会产生三个平均向量：
这些平均向量可通过如下代码计算，其中计算了三个标签各自的平均向量：
>>> 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(f'MV {label}: {mean_vecs[label - 1]}\n')
MV 1: [ 0.9066 -0.3497 0.3201 -0.7189 0.5056 0.8807 0.9589 -0.5516
0.5416 0.2338 0.5897 0.6563 1.2075]
MV 2: [-0.8749 -0.2848 -0.3735 0.3157 -0.3848 -0.0433 0.0635 -0.0946
0.0703 -0.8286 0.3144 0.3608 -0.7253]
MV 3: [ 0.1992 0.866 0.1682 0.4148 -0.0451 -1.0286 -1.2876 0.8287
-0.7795 0.9649 -1.209 -1.3622 -0.4013]
使用平均向量我们就可以计算出类内散度矩阵SW：
这是通过单个类i的单独散度矩阵Si加和而得：
>>> 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_train_std[y_train == 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: '
... f'{S_W.shape[0]}x{S_W.shape[1]}')
Within-class scatter matrix: 13x13
我们在计算散度矩阵时的假设是训练集中的类标签是均匀分布的。但如果打印出类标签数，会发现不符合这一假定：
>>> print('Class label distribution:',
... np.bincount(y_train)[1:])
Class label distribution: [41 50 33]
因此在加总散度矩阵SW前我们会希望缩放各散度矩阵Si。将散度矩阵除以类样本数ni时，可以看到计算散度矩阵实际上和计算协方差矩阵$\Sigma_i$一样-协方差矩阵是散度矩阵的归一化版本：
计算缩放后的类内散度矩阵的代码如下：
>>> 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: '
... f'{S_W.shape[0]}x{S_W.shape[1]}')
Scaled within-class scatter matrix: 13x13
在计算好缩放后类内散度矩阵（或协方差矩阵）后，我们可以进入下一步计算类间散度矩阵SB：
这里m是所计算的整体均值，包含所有c类中的样本：
>>> mean_overall = np.mean(X_train_std, axis=0)
>>> mean_overall = mean_overall.reshape(d, 1)
>>> d = 13 # number of features
>>> S_B = np.zeros((d, d))
>>> for i, mean_vec in enumerate(mean_vecs):
... n = X_train_std[y_train == i + 1, :].shape[0]
... mean_vec = mean_vec.reshape(d, 1) # make column vector
... S_B += n * (mean_vec - mean_overall).dot(
... (mean_vec - mean_overall).T)
>>> print('Between-class scatter matrix: '
... f'{S_B.shape[0]}x{S_B.shape[1]}')
Between-class scatter matrix: 13x13
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 descending order:\n')
>>> for eigen_val in eigen_pairs:
... print(eigen_val[0])
Eigenvalues in descending order:
349.617808906
172.76152219
3.78531345125e-14
2.11739844822e-14
1.51646188942e-14
1.51646188942e-14
1.35795671405e-14
1.35795671405e-14
7.58776037165e-15
5.90603998447e-15
5.90603998447e-15
2.25644197857e-15
0.0
在LDA中，线性判别的数量最多为c – 1，其中c是类标签的数量，因为类内类间散度矩阵SB，是秩为1或更低的c矩阵的和。可以看到我们只有两个非零特征值（特征值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, align='center',
... label='Individual discriminability')
>>> plt.step(range(1, 14), cum_discr, where='mid',
... label='Cumulative discriminability')
>>> plt.ylabel('"Discriminability" ratio')
>>> plt.xlabel('Linear Discriminants')
>>> plt.ylim([-0.1, 1.1])
>>> plt.legend(loc='best')
>>> plt.tight_layout()
>>> plt.show()
在图5.9中可以看出，前两个线性判别就捕获取了葡萄酒训练集100%的有用信息：
图5.9：头两个判别捕获了100%的有用信息
我们来堆叠最具判别力的两个特征向量列，创建变换矩阵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.1481 -0.4092]
[ 0.0908 -0.1577]
[-0.0168 -0.3537]
[ 0.1484 0.3223]
[-0.0163 -0.0817]
[ 0.1913 0.0842]
[-0.7338 0.2823]
[-0.075 -0.0102]
[ 0.0018 0.0907]
[ 0.294 -0.2152]
[-0.0328 0.2747]
[-0.3547 -0.0124]
[-0.3915 -0.5958]]
使用我们在前一小节中创建的变换矩阵W，现在可以通过矩阵相乘来变换训练集：
X′ = XW
>>> X_train_lda = X_train_std.dot(w)
>>> colors = ['r', 'b', 'g']
>>> markers = ['o', 's', '^']
>>> 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] * (-1),
... c=c, label= f'Class {l}', marker=m)
>>> plt.xlabel('LD 1')
>>> plt.ylabel('LD 2')
>>> plt.legend(loc='lower right')
>>> plt.tight_layout()
>>> plt.show()
从图5.10中可以看出，三个葡萄酒类现在可以完美地在新特征子空间上进行线性分割：
图5.10：在数据投射到前两个判别上后葡萄酒类可完美分割
前面的分步实现对于理解LDA内部原理并掌握LDA和PCA之间的差别是一个很好的练习。下面我们来学习scikit-learn中LDA类的实现：
>>> # the following import statement is one line
>>> 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(multi_class='ovr', random_state=1,
... solver='lbfgs')
>>> 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='lower left')
>>> plt.tight_layout()
>>> plt.show()
在图5.11中可以看出逻辑回归模型对类2中的一个样本分错了类：
图5.11：逻辑回归模型分错了一个类
通过降低正则化强度，我们可以移动决策边界，以使用逻辑回归模型正确分类训练集中的所有样本。但更重要的是，我们来看看它对测试集的结果：
>>> 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='lower left')
>>> plt.tight_layout()
>>> plt.show()
在图5.12中可以看到，逻辑回归分类器只需使用一个二维特征子空间就可以很精准地对测试集样本分类，根本不需要原来的12个葡萄酒特征：
图5.12：逻辑回归模型可很好地处理测试数据
在前面一节中，我们讲解了用线性变换技术，比如PCA和LDA进行特征提取。本节中，我们会讨论为什么考虑使用非线性降维技术是有价值的。
尤其值得重点介绍的一种非线性降维技术是t-分布随机近邻嵌套(t-SNE) ，因为文献中经常使用它来以二维或三维可视化高维数据集。我们会学习如何应用t-SNE来以二维特征空间对手写图片绘图。
很多机器学习算法都假定输入数据是线性可分割的。我们学习过感知机甚至要求完全线性可分割的训练数据收敛。我们至此学习过的其它算法假定无法完全线性分割的原因是噪声：Adaline、逻辑回归和（标准）SVM等皆是如此。
但如果要处理非线性问题，这在现实世界的应用中经常会碰到，进行降维的线性变换技术，比如PCA和LDA可能就不是最佳选择了：
图5.13：线性和非线性问题之间的不同
scikit-learn库实现了很多用于非线性降维的高级技术，这些不在本书的讲解范围内。感举的读者可以在scikit-learn当前实现中进行总览，辅以一些图示：http://scikit-learn.org/stable/modules/manifold.html。
非线性降维技术的开发和应用常被称作流形学习（manifold learning），流形是指低维拓扑空间嵌套在高维空间中。流形学习的算法需要捕获数据的复杂结构，以投射到低维空间上并且保留数据点间的关系。
经典的流形学习示例是图5.14中所示的3维瑞士卷：
图5.14：三维瑞士卷投射到更低的二维空间上
虽然非线性降维和流形学习算法很强大，但需要注意这些技术出了名的难用，并且如果选择了不理想的超参数，它们可能会弊大于利。这种难度背后的原因是我们常常处理尚未能可视化的高维数据集，其中的结果尚不清晰（不同于图5.14中的瑞士卷示例）。此外，除非我们将数据集投射到二维或三维上（通常不足以捕获更复杂的关系），否则会很难甚至是无法评估结果的质量。因此，很多人仍然依靠PCA和LDA这种更简单的技术来实现降维。
在介绍非线性降维并讨论了它的一些挑战后，我们来学习一个涉及t-SNE的实操案例，它通常用于将复杂数据可视化至二维或三维。
简单地说，t-SNE是基于高维（原始）特征空间中成对距离的一些模型数据点。然后，它在新的低维空间中找到接近原始空间成对距概率分布的成对距离概率分布。或者换句话说，t-SNE学习将数据点嵌套至低维空间，使得原始空间的成对距离得以保留。可以原研究论文中找到该方法是多详细的说明：Visualizing data using t-SNE by Maaten and Hinton, Journal of Machine Learning Research, 2018 (https://www.jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf)。但就像研究论文标题所说的，t-SNE是一种意在用于可视化的技术，因为它需要整个数据集来进行投射。因其直接投射各点（不同于PCA，它并不涉及到投射矩阵），我们无法对新数据点应用t-SNE。
以下代码快速演示了如何将t-SNE应用于一个64-维数据集。首先，我们通过scikit-learn加载Digits数据集，其中包含一些低分辨率的手写数字（数字0-9）：
>>> from sklearn.datasets import load_digits
>>> digits = load_digits()
数字是8×8的灰度图片。以下代码绘制数据集中的前4组图片，其中总共包含1,797个图像：
>>> fig, ax = plt.subplots(1, 4)
>>> for i in range(4):
>>> ax[i].imshow(digits.images[i], cmap='Greys')
>>> plt.show()
从图5.15中可以看出，图像的像素较低，8×8像素（也就是每张图64像素）：
图5.15：低分辨率手写数据图像
注意digits.data属性让我们可以访问该数据集的平坦版本，其中样本表现为行，列对应像素：
>>> digits.data.shape
(1797, 64)
接下来我们对新变量X_digits赋特征（像素）并将标签赋值新的变量y_digits：
>>> y_digits = digits.target
>>> X_digits = digits.data
然后，我们从scikit-learn导入t-SNE类并拟合一个新的tsne对象。使用fit_transform，我们用一步执行了t-SNE拟合及数据变换：
>>> from sklearn.manifold import TSNE
>>> tsne = TSNE(n_components=2, init='pca',
... random_state=123)
>>> X_digits_tsne = tsne.fit_transform(X_digits)
通过这段代码，我们将64维数据集投射到了一个二维空间上。我们指定了init='pca'，它使用PCA初始化了t-SNE嵌套，因为这是研究文章中所推荐的：nitialization is critical for preserving global data structure in both t-SNE and UMAP by Kobak and Linderman, Nature Biotechnology Volume 39, pages 156–157, 2021 (https://www.nature.com/articles/s41587-020-00809-z)。
注意t-SNE还包含其它超参数，比如困惑度（perplexity）和学习率（通常称为epsilon），在我们的示例中省略了（使用了scikit-learn的默认值）。在实操中，我们推荐读者也研究下这些参数。更多有关这些参数的信息及其对结果的影响请见这篇高质量的文章：How to Use t-SNE Effectively by Wattenberg, Viegas, and Johnson, Distill, 2016 (https://distill.pub/2016/misread-tsne/)。
最后，我们使用如下代码可视化二维t-SNE嵌套：
>>> import matplotlib.patheffects as PathEffects
>>> def plot_projection(x, colors):
... f = plt.figure(figsize=(8, 8))
... ax = plt.subplot(aspect='equal')
... for i in range(10):
... plt.scatter(x[colors == i, 0],
... x[colors == i, 1])
... for i in range(10):
... xtext, ytext = np.median(x[colors == i, :], axis=0)
... txt = ax.text(xtext, ytext, str(i), fontsize=24)
... txt.set_path_effects([
... PathEffects.Stroke(linewidth=5, foreground="w"),
... PathEffects.Normal()])
>>> plot_projection(X_digits_tsne, y_digits)
>>> plt.show()
同PCA一样，t-SNE是一种无监督方法，在上面的代码中，我们只使用了类标签y_digits (0-9)是为了通过函数颜色参数进行可视化。Matplotlib的PathEffects用于可视化，（通过np.median）使得类标签显示在属于各自数字数据点的中央 。结果图如下：
图5.16：如何用t-SNE将手写数字嵌套到二维特征空间的可视化
可以看到，t-SNE可以很好地分割不同的数字（类），虽然仍不完美。通过调优超参数很可能达到更好的分离。但模糊的手写体可能会使用一定程度的类交叉不可避免。例如，通过查看单独的图片，我们可能发现有些数字3的图片看起来很像9，诸如此类。
统一流形逼近和投影
另一种著名的可视化技术是统一流形逼近和投影（uniform manifold approximation and projection (UMAP)）。UMAP可以产生类似t-SNE一样的好结果（比如，参见早前引用的Kobak和Linderman的论文），且通常更快，它也可用于投射新数据，这使用其在机器学习中用作降维技术时更具吸引力，类似PCA。感兴趣的读者可以在原论文中找到UMAP更多的信息：UMAP: Uniform manifold approximation and projection for dimension reduction by McInnes, Healy, and Melville, **2018 **(https://arxiv.org/abs/1802.03426) 。 兼容scikit-learn的UMAP实现请见https://umap-learn.readthedocs.io。
本章中，我们学习了用于特征提取的两个基本降维技术：PCA和LDA。使用PCA，我们将数据投射到更低维的子空间上，以沿正交的特征轴最大化方差，同时忽略类标签。LDA则与PCA不同，是一种监督降维技术，也就是说它会考虑训练集中的类信息，尝试在线性特征空间中最大化类的分割。最后，我们还学习了t-SNE，这是一种非线性特征提取技术，可用于在二维或三维中可视化数据。
配备了基本预处理技术PCA和LDA，我们为在下一章中学习有效结合各种预处理技术最佳实践及评估各种模型表现做好了准备。