唐宇迪的《跟着迪哥学：Python数据分析与机器学习实战》，2019年9月出版，本系列为读书笔记。主要是为了系统整理，加深记忆。 第15章降维算法。如果拿到的数据特征过于庞大，一方面会使得计算任务变得繁重；另一方面，如果数据特征还有问题，可能会对结果造成不利的影响。降维是机器学习领域中经常使用的数据处理方法，一般通过某种映射方法，将原始高维空间中的数据点映射到低维度的空间中，本章将从原理和实践的角度介绍两种经典的降维算法——线性判别分析和主成分分析。
第15章降维算法
如果拿到的数据特征过于庞大，一方面会使得计算任务变得繁重；另一方面，如果数据特征还有问题，可能会对结果造成不利的影响。降维是机器学习领域中经常使用的数据处理方法，一般通过某种映射方法，将原始高维空间中的数据点映射到低维度的空间中，本章将从原理和实践的角度介绍两种经典的降维算法——线性判别分析和主成分分析。
15.1线性判别分析
线性判别式分析（Linear Discriminant Analysis,LDA），也叫作Fisher线性判别（Fisher Linear Discriminant,FLD），最开始用于处理机器学习中的分类任务，但是由于其对数据特征进行了降维投影，使其成为一种经典的降维方法。
15.1.1降维原理概述
线性判别分析属于有监督学习算法，也就是数据中必须要有明确的类别标签，它不仅能用来降维，还可以处理分类任务，不过，更多用于降维。下面通过一个小例子来感受下降维的作用。这个游戏需要通过不断地寻找最合适的投影面，来观察原始物体的形状，如图15-1所示。降维任务与之非常相似，也是通过找到最合适的投影方向，使得原始数据更容易被计算机理解并利用。
图15-1 投影的意义
接下来，通过实例说明降维的过程。假设有两类数据点，如图15-2所示。由于数据点都是二维特征，需要将其降到一维，也就是需要找到一个最合适的投影方向把这些数据点全部映射过去。图15-2（a）、（b）分别是选择两个投影方向后的结果，那么，究竟哪一个更好呢？
图15-2 降维的目的
从投影结果上观察，图15-2（a）中的数据点经过投影后依旧有一部分混在一起，区别效果有待提高。图15-2（b）中的数据点经过投影后，没有混合，区别效果比图15-2（a）更好一些。因此，我们当然会选择图15-2（b）所示的降维方法，由此可见，降维不仅要压缩数据的特征，还需要寻找最合适的方向，使得压缩后的数据更有利用价值。
由图15-2可知，线性判别分析的原理可以这样理解：任务目标就是要找到最合适的投影方向，这个方向可以是多维的。
为了把降维任务做得更圆满，提出了两个目标。
1.对于不同类别的数据点，希望其经过投影后能离得越远越好，也就是两类数据点区别得越明显越好，不要混在一起。
2.对于同类别的数据点，希望它们能更集中，离组织的中心越近越好。
接下来的任务就是完成这两个目标，这也是线性判别分析的核心优化目标，降维任务就是找到能同时满足这两个目标的投影方向。
15.1.2优化的目标
投影就是通过矩阵变换的方式把数据映射到最适合做分类的方向上：
其中，x表示当前数据所在空间，也就是原始数据；y表示降维后的数据。最终的目标也很明显，就是找到最合适的变换方向，即求解出参数W。除了降维任务中提出的两个目标，还需要定义一下距离这个概念，例如“扎堆”该怎么体现呢？这里用数据点的均值来表示其中心位置，如果每一个数据点都离中心很近，它们就“扎堆”在一起了。中心点位置计算方法如下：
对于多分类问题，也可以得到各自类别的中心点，不同类别需要各自计算，这就是为什么要强调线性判别分析是一个有监督问题，因为需要各个类别分别进行计算，所以每一个数据点是什么类别，必须在标签中给出。
在降维算法中，其实我们更关心的并不是原始数据集中数据点的扎堆情况，而是降维后的结果，因此，可知投影后中心点位置为：
由式（15.3）可以得到投影后的中心点计算方法，按照之前制定的目标，对于一个二分类任务来说，应当使得这两类数据点的中心离得越远越好，这样才能更好地区分它们：
现在可以把当作目标函数，目标是希望其值能够越大越好，但是只让不同类别投影后的中心点越远可以达到我们期望的结果吗？
对于图15-3所示的样本数据，假设只能在x1和x2两方向进行投影，如果按照之前定义的J(w)，显然x1方向更合适，但是投影后两类数据点依旧有很多重合在一起，而x2方向上的投影结果是两类数据点重合较少；因此，x2方向更好。
这个问题就涉及要优化的另一个目标，不仅要考虑不同类别之间要区分开，还要考虑同类样本点应当尽可能聚集在一起。显然在图15-3中，x1方向不满足这个条件，因为在x1方向上，同类样本变得更分散，不够集中。
图15-3 投影方向选择
我们还可以使用另一个度量指标——散列值（scatter），表示同类数据样本点的离散程度，定义如下：
其中，y表示经过投影后的数据点，从式（15.5）中可以看出，散列值表示样本点的密集程度，其值越大，表示越分散；反之，则越集中。定义好要优化的两个目标后，接下来就是求解了。
15.1.3线性判别分析求解
上一小节已经介绍了降维后想要得到的目标，现在把它们综合在一起，但是优化的目标有两个，那么如何才能整合它们呢？
既然要最大化不同类别之间的距离，那就把它当作分子；最小化同类样本之间的离散程度，那就把它当作分母，最终整体的J(W)依旧求其极大值即可。
在公式推导过程中，牢记最终的要求依旧是寻找最合适的投影方向，先把散列值公式展开：
为了化简方便，则令：
式（15.8）称为散布矩阵（scatter matrices）。
由此可以定义类内散布矩阵为：
Sw=S1+S2 (15.9)
将式（15.9）代入式（15.7），可得：
对式（15.6）分子进行展开可得：
其中，SB为类间散布矩阵。
目标函数J(w)最终可以表示为：
对于散列矩阵SB和SW，只要有数据和标签即可求解。但是如何求解最终的结果呢？如果对分子和分母同时求解，就会有无穷多解，通用的解决方案是先固定分母，经过放缩变换后，将其值限定为1，则令：
在此条件下求WTSBW的极大值点，利用拉格朗日乘子法可得：
既然目标是找投影方向，也就是w，将式（15.14）左右两边同时乘以Sw-1可得：
观察一下式（15.15），它与线性代数中的特征向量有点像，如果把Sw-1当作一个整体，那么w就是其特征向量，问题到此迎刃而解。在线性判别分析中，其实只需要得到类内和类间散布矩阵，然后求其特征向量，就可以得到投影方向，然后，只需要对数据执行相应的矩阵变换，就完成全部降维操作。
15.1.4Python实现线性判别分析降维
下面要在非常经典的“鸢尾花”数据集上使用线性判别分析完成降维任务。鸢尾花数据集可从该链接下载：https://archive.ics.uci.edu/ml/datasets/Iris，也可以从sklearn库内置的数据集中获取。数据集中含有3类、共150条鸢尾花基本数据，其中山鸢尾、变色鸢尾、维吉尼亚鸢尾各有50条数据，每条数据包括萼片长度（单位：厘米）、萼片宽度、花瓣长度、花瓣宽度4种特征。
首先读取数据集，代码如下：
1 # 自己来定义列名 2 feature_dict = {i:label for i,label in zip( 3 range(4), 4 (\'sepal length in cm\', 5 \'sepal width in cm\', 6 \'petal length in cm\', 7 \'petal width in cm\', ))} 8 9 label_dict = {i:label for i,label in zip( 10 range(1,4), 11 (\'Setosa\', 12 \'Versicolor\', 13 \'Virginica\' 14 ))} 15 import pandas as pd 16 # 数据读取，大家也可以先下载下来直接读取 17 df = pd.io.parsers.read_csv( 18 filepath_or_buffer=\'https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data\', 19 header=None, 20 sep=\',\', 21 ) 22 # 指定列名 23 df.columns = [l for i,l in sorted(feature_dict.items())] + [\'class label\'] 24 25 df.head()
数据集共有150条数据，每条数据有4个特征，现在需要将四维特征降成二维。观察输出结果可以发现，其特征已经是数值数据，不需要做额外处理，但是需要转换一下标签：
1 from sklearn.preprocessing import LabelEncoder 2 3 X = df[[\'sepal length in cm\',\'sepal width in cm\',\'petal length in cm\',\'petal width in cm\']].values 4 y = df[\'class label\'].values 5 6 # 制作标签{1: \'Setosa\', 2: \'Versicolor\', 3:\'Virginica\'} 7 enc = LabelEncoder() 8 label_encoder = enc.fit(y) 9 y = label_encoder.transform(y) + 1
上述代码使用了sklearn工具包中的LabelEncoder用于快速完成标签转换，可以发现基本上所有sklearn中的数据处理操作都是分两步走，先fit再transform，变换结果如图15-4所示。
在计算过程中需要基于均值来判断距离，因此先要对数据中各个特征求均值，但是只求4个特征的均值能满足要求吗？不要忘记任务中还有3种花，相当于3个类别，所以也要对每种花分别求其各个特征的均值：
图15-4 标签转换
1 import numpy as np
2 #设置小数点的位数
3 np.set_printoptions(precision=4)
4 #这里会保存所有的均值
5 mean_vectors = []
6 # 要计算3个类别
7 for cl in range(1,4):
8 # 求当前类别各个特征均值
9 mean_vectors.append(np.mean(X[y==cl], axis=0))
10 print(\'均值类别 %s: %s\n\' %(cl, mean_vectors[cl-1]))
均值类别 1: [5.006 3.418 1.464 0.244]
均值类别 2: [5.936 2.77 4.26 1.326]
均值类别 3: [6.588 2.974 5.552 2.026]
接下来计算类内散布矩阵：
1 # 原始数据中有4个特征
2 S_W = np.zeros((4,4))
3 # 要考虑不同类别，自己算自己的
4 for cl,mv in zip(range(1,4), mean_vectors):
5 class_sc_mat = np.zeros((4,4))
6 # 选中属于当前类别的数据
7 for row in X[y == cl]:
8 # 这里相当于对各个特征分别进行计算，用矩阵的形式
9 row, mv = row.reshape(4,1), mv.reshape(4,1)
10 # 跟公式一样
11 class_sc_mat += (row-mv).dot((row-mv).T)
12 S_W += class_sc_mat
13 print(\'类内散布矩阵:\n\', S_W)
类内散布矩阵:
[[38.9562 13.683 24.614 5.6556]
[13.683 17.035 8.12 4.9132]
[24.614 8.12 27.22 6.2536]
[ 5.6556 4.9132 6.2536 6.1756]]
继续计算类间散布矩阵：
式中，m为全局均值；mi为各个类别的均值；Ni为样本个数。
1 # 全局均值 2 overall_mean = np.mean(X, axis=0) 3 # 构建类间散布矩阵 4 S_B = np.zeros((4,4)) 5 # 对各个类别进行计算 6 for i,mean_vec in enumerate(mean_vectors): 7 #当前类别的样本数 8 n = X[y==i+1,:].shape[0] 9 mean_vec = mean_vec.reshape(4,1) 10 overall_mean = overall_mean.reshape(4,1) 11 # 如上述公式进行计算 12 S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T) 13 14 print(\'类间散布矩阵:\n\', S_B)
类间散布矩阵: [[ 63.2121 -19.534 165.1647 71.3631] [-19.534 10.9776 -56.0552 -22.4924] [165.1647 -56.0552 436.6437 186.9081] [ 71.3631 -22.4924 186.9081 80.6041]]
得到类内和类间散布矩阵后，还需将它们组合在一起，然后求解矩阵的特征向量：
1 #求解矩阵特征值，特征向量 2 eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B)) 3 # 拿到每一个特征值和其所对应的特征向量 4 for i in range(len(eig_vals)): 5 eigvec_sc = eig_vecs[:,i].reshape(4,1) 6 print(\'\n特征向量 {}: \n{}\'.format(i+1, eigvec_sc.real)) 7 print(\'特征值 {:}: {:.2e}\'.format(i+1, eig_vals[i].real))
特征向量 1: [[ 0.2049] [ 0.3871] [-0.5465] [-0.7138]] 特征值 1: 3.23e+01 特征向量 2: [[-0.009 ] [-0.589 ] [ 0.2543] [-0.767 ]] 特征值 2: 2.78e-01 特征向量 3: [[-0.8379] [ 0.1696] [ 0.1229] [ 0.5041]] 特征值 3: -4.13e-15 特征向量 4: [[ 0.2 ] [-0.3949] [-0.4567] [ 0.7717]] 特征值 4: 1.20e-14
输出结果得到4个特征值和其所对应的特征向量。特征向量直接观察起来比较麻烦，因为投影方向在高维上很难理解；特征值还是比较直观的，这里可以认为特征值代表的是其所对应特征向量的重要程度，也就是特征值越大，其所对应的特征向量就越重要，所以接下来可以对特征值按大小进行排序，排在前面的越重要，排在后面的就没那么重要了。
1 #特征值和特征向量配对 2 eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))] 3 4 # 按特征值大小进行排序 5 eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True) 6 7 print(\'特征值排序结果:\n\') 8 for i in eig_pairs: 9 print(i[0])
特征值排序结果:
32.27195779972981
0.27756686384004264
1.1953730364935478e-14
4.1311796919088535e-15
1 print(\'特征值占总体百分比:\n\') 2 eigv_sum = sum(eig_vals) 3 for i,j in enumerate(eig_pairs): 4 print(\'特征值 {0:}: {1:.2%}\'.format(i+1, (j[0]/eigv_sum).real))
特征值占总体百分比: 特征值 1: 99.15% 特征值 2: 0.85% 特征值 3: 0.00% 特征值 4: 0.00%
可以看出，打印出来的结果差异很大，第一个特征值占据总体的99.15%，第二个特征值只占0.85%，第三和第四个特征值，看起来微不足道。这表示对鸢尾花数据进行降维时，可以把特征数据降到二维甚至一维，但没必要降到三维。
既然已经有结论，选择把数据降到二维，只需选择特征值1、特征值2所对应的特征向量即可：
1 W = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1))) 2 print(\'矩阵W:\n\', W.real)
矩阵W: [[ 0.2049 -0.009 ] [ 0.3871 -0.589 ] [-0.5465 0.2543] [-0.7138 -0.767 ]]
这也是最终所需的投影方向，只需和原始数据组合，就可以得到降维结果：
1 # 执行降维操作 2 X_lda = X.dot(W) 3 X_lda.shape
(150, 2)
现在可以看到数据维度从原始的（150,4）降到（150,2），到此就完成全部的降维工作。接下来对比分析一下降维后结果，为了方便可视化展示，在原始四维数据集中随机选择两维进行绘图展示：
1 from matplotlib import pyplot as plt 2 # 可视化展示 3 def plot_step_lda(): 4 5 ax = plt.subplot(111) 6 for label,marker,color in zip( 7 range(1,4),(\'^\', \'s\', \'o\'),(\'blue\', \'red\', \'green\')): 8 9 plt.scatter(x=X[:,0].real[y == label], 10 y=X[:,1].real[y == label], 11 marker=marker, 12 color=color, 13 alpha=0.5, 14 label=label_dict[label] 15 ) 16 17 plt.xlabel(\'X[0]\') 18 plt.ylabel(\'X[1]\') 19 20 leg = plt.legend(loc=\'upper right\', fancybox=True) 21 leg.get_frame().set_alpha(0.5) 22 plt.title(\'Original data\') 23 24 # 把边边角角隐藏起来 25 plt.tick_params(axis="both", which="both", bottom="off", top="off", 26 labelbottom="on", left="off", right="off", labelleft="on") 27 28 # 为了看的清晰些，尽量简洁 29 ax.spines["top"].set_visible(False) 30 ax.spines["right"].set_visible(False) 31 ax.spines["bottom"].set_visible(False) 32 ax.spines["left"].set_visible(False) 33 34 plt.grid() 35 plt.tight_layout 36 plt.show() 37 38 plot_step_lda()
从上述输出结果可以发现，如果对原始数据集随机取两维数据，数据集并不能按类别划分开，很多数据都堆叠在一起（尤其是图中方块和圆形数据点）。再来看看降维后的数据点分布，绘图代码保持不变，只需要传入降维后的两维数据即可，可以生成图15-5的输出。
1 from matplotlib import pyplot as plt 2 # 可视化展示 3 def plot_step_lda(): 4 5 ax = plt.subplot(111) 6 for label,marker,color in zip( 7 range(1,4),(\'^\', \'s\', \'o\'),(\'blue\', \'red\', \'green\')): 8 9 plt.scatter(x=X_lda[:,0].real[y == label], 10 y=X_lda[:,1].real[y == label], 11 marker=marker, 12 color=color, 13 alpha=0.5, 14 label=label_dict[label] 15 ) 16 17 plt.xlabel(\'LD1\') 18 plt.ylabel(\'LD2\') 19 20 leg = plt.legend(loc=\'upper right\', fancybox=True) 21 leg.get_frame().set_alpha(0.5) 22 plt.title(\'LDA on iris\') 23 24 # 把边边角角隐藏起来 25 plt.tick_params(axis="both", which="both", bottom="off", top="off", 26 labelbottom="on", left="off", right="off", labelleft="on") 27 28 # 为了看的清晰些，尽量简洁 29 ax.spines["top"].set_visible(False) 30 ax.spines["right"].set_visible(False) 31 ax.spines["bottom"].set_visible(False) 32 ax.spines["left"].set_visible(False) 33 34 plt.grid() 35 plt.tight_layout 36 plt.show() 37 38 plot_step_lda()
图15-5 线性判别分析降维后数据点分布
可以明显看到，坐标轴变成LD1与LD2，这就是降维后的结果，从数据点的分布来看，混杂在一起的数据不多，划分起来就更容易。这就是经过一步步计算得到的最终降维结果。
当拿到一份规模较大的数据集时，如何选定降维的维数呢？一方面可以通过观察特征值排序结果来决定，另一方面还是要通过实验来进行交叉验证。
我们肯定希望用更高效、稳定的方法来完成一个实际任务，再来看看sklearn工具包中如何调用线性判别分析进行降维：
1 from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA 2 3 # LDA 4 sklearn_lda = LDA(n_components=2) 5 X_lda_sklearn = sklearn_lda.fit_transform(X, y)
1 def plot_scikit_lda(X, title): 2 3 ax = plt.subplot(111) 4 for label,marker,color in zip( 5 range(1,4),(\'^\', \'s\', \'o\'),(\'blue\', \'red\', \'green\')): 6 7 plt.scatter(x=X[:,0][y == label], 8 y=X[:,1][y == label] * -1, # flip the figure 9 marker=marker, 10 color=color, 11 alpha=0.5, 12 label=label_dict[label]) 13 14 plt.xlabel(\'LD1\') 15 plt.ylabel(\'LD2\') 16 17 leg = plt.legend(loc=\'upper right\', fancybox=True) 18 leg.get_frame().set_alpha(0.5) 19 plt.title(title) 20 21 # hide axis ticks 22 plt.tick_params(axis="both", which="both", bottom="off", top="off", 23 labelbottom="on", left="off", right="off", labelleft="on") 24 25 # remove axis spines 26 ax.spines["top"].set_visible(False) 27 ax.spines["right"].set_visible(False) 28 ax.spines["bottom"].set_visible(False) 29 ax.spines["left"].set_visible(False) 30 31 plt.grid() 32 plt.tight_layout 33 plt.show()
1 plot_scikit_lda(X_lda_sklearn, title=\'Default LDA via scikit-learn\')
很简单，仅仅两步外加一个指定降维到两维的参数即可，上述代码可以生成图15-6的输出。
可以发现，使用sklearn工具包降维后的结果与自己一步步计算的结果完全一致。
图15-6 sklearn工具包降维结果
15.2主成分分析
主成分分析（Principal Component Analysis，PCA）是在降维中使用特别广泛的算法。在使用主成分分析降维的时候没有束缚，不像线性判别分析，必须要有数据标签，只要拿到数据，没有标签也可以用主成分分析进行降维。所以应该先有一个直观的认识，主成分分析本质上属于无监督算法，这也是它流行的主要原因。
15.2.1PCA降维基本知识点
既然不需要标签，就很难去分析类间与类内距离因素，那么该怎么办呢？PCA的基本思想就是方差，可以想象一下哪些特征更有价值？应当是那些区别能力更强的特征。例如我们想比较两个游戏玩家的战斗力水平。第一个特征是其所在帮派等级：A玩家，5级帮派；B玩家，4级帮派，A、B玩家的帮派等级看起来差别不大。第二个特征是其充值金额：A玩家，10000；B玩家，100。A、B玩家的充值金额的差距好像有些大。通过这两个特征就可以预估一下哪个玩家战斗力更强，答案肯定是A玩家。
现在再来观察一下这两个特征，帮派等级似乎相差不大，不能拉开差距，但是充值金额的差异却很大。我们希望得到充值金额这种能把不同玩家区分开的特征。在数学上可以用方差来描述这种数据的离散程度，所以在主成分析中主要依靠方差。
为了让大家更好地理解主成分分析，下面介绍一些基本概念。
（1）向量的表示。假设有向量（3，2），如图15-7所示。为什么向量可以表示为（3，2），这是在直角坐标系中的表示，如果坐标系变了，向量的表示形式也要发生变换。
实际上该向量可以表示成线性组合3·(1，0)T+2·(0，1)T，其中（1，0）和（0，1）就称为二维空间中的一组基。
（3）基变换。大家常见的坐标系都是正交的，即内积为0，两两相互垂直，并且线性无关。为什么基都是这样的呢？如果不垂直，那么肯定线性相关，能用一个表示另一个，此时基就会失去意义，所以基的出发点就是要正交。
图15-7 向量的组成
基也可以进行变换，将一个向量从一组基变换到另一组基中。例如新的坐标系的两个基分别是(
)，因此向量（3，2）映射到这个新的坐标系中，可以通过下面变换实现：
（3）方差和协方差。方差（variance）相当于特征辨识度，其值越大越好。协方差（covariance）就是不同特征之间的相关程度，协方差的计算式式为：
如果两个变量的变化趋势相同，例如随着身高的增长，体重也增长，此时它们的协方差值就会比较大，表示正相关。而方差又描述了各自的辨识能力，接下来就要把这些知识点穿插在一起。
15.2.2PCA优化目标求解
对于降维任务，无非就是将原始数据特征投影到一个更合适的空间，结合基的概念，这就相当于由一组基变换到另一组基，变换的过程要求特征变得更有价值，也就是方差能够更大。所以现在已经明确基本目标了：找到一组基，使得变换后的特征方差越大越好。
假设找到了第一个合适的投影方向，这个方向能够使得方差最大，对于降维任务来说，一般情况下并不是降到一维，接下来肯定要找方差第二大的方向。方差第二大的方向理论上应该与第一方向非常接近，甚至重合，这样才能保证方差最大，如图15-8所示。
在这种情况下，看似可以得到无数多个方差非常大的方向，但是想一想它们能组成想要的基吗？不能，因为没有满足基的最基本要求——线性无关，也就是相互垂直正交。所以在寻找方差最大的方向的同时，还要使得各个投影方向能够正交，即协方差应当等于0，表示完全独立无关。所以在选择基的时候，一方面要尽可能地找方差的最大方向，另一方面要在其正交方向上继续寻找方差第二大的方向，以此类推。
图15-8 方差方向选择
解释PCA中要求解的目标后，接下来就是在数学上将它表达出来。先来看一下协方差矩阵，为了简便，可以把数据中各个特征的均值默认为0，也可以认为数据已经进行过标准化处理。其计算式如下：
其中，X为实际的数据。包含2个特征a和b，一共有m个样本。
此时协方差矩阵为：
先观察一下协方差矩阵结果，其主对角线上的元素就是各个特征的方差（均值为0时），而非主对角线的上元素恰好是特征之间的协方差。按照目标函数的要求，首先应当使得方差越大越好，并且确保协方差为0，这就需要对协方差矩阵做对角化。
从一个n行n列的实对称矩阵中一定可以找到n个单位正交特征向量E=(e1,e2 ,… ,en)，以完成对角化的操作：
式（15.21）中的协方差矩阵恰好满足上述要求。假设需要将一组N维向量降为K维（K大于0，小于N），目标是选择K个单位正交基，使原始数据变换到这组基上后，各字段两两间协方差为0，各字段本身的方差尽可能大。当得到其协方差矩阵后，对其进行对角化操作，即可使得除主对角线上元素之外都为0。
其中对角线上的元素就是矩阵的特征值，这与线性判别分析很像，还是先把特征值按从大到小的顺序进行排列，找到前K个最大特征值对应的特征向量，接下来就是进行投影变换。
按照给定PCA优化目标，基本流程如下。
- 第①步：数据预处理，只有数值数据才可以进行PCA降维。
- 第②步：计算样本数据的协方差矩阵。
- 第③步：求解协方差矩阵的特征值和特征向量。
- 第④步：将特征值按照从大到小的顺序排列，选择其中较大的K个，然后将其对应的K个特征向量组成投影矩阵。
- 第⑤步：将样本点投影计算，完成PCA降维任务。
15.2.3Python实现PCA降维
接下来通过一个实例介绍如何使用PCA处理实际问题，同样使用鸢尾花数据集，目的依旧是完成降维任务，下面就来看一下PCA是怎么实现的。
第①步：导入数据。
1 import numpy as np 2 import pandas as pd 3 4 # 读取数据集 5 df = pd.read_csv(\'iris.data\') 6 # 原始数据没有给定列名的时候需要我们自己加上 7 df.columns=[\'sepal_len\', \'sepal_wid\', \'petal_len\', \'petal_wid\', \'class\'] 8 df.head()
sepal_len sepal_wid petal_len petal_wid class 0 4.9 3.0 1.4 0.2 Iris-setosa 1 4.7 3.2 1.3 0.2 Iris-setosa 2 4.6 3.1 1.5 0.2 Iris-setosa 3 5.0 3.6 1.4 0.2 Iris-setosa 4 5.4 3.9 1.7 0.4 Iris-setosa
第②步：展示数据特征。
1 # 把数据分成特征和标签 2 3 X = df.iloc[:,0:4].values 4 y = df.iloc[:,4].values 5 6 from matplotlib import pyplot as plt 7 8 # 展示我们标签用的 9 label_dict = {1: \'Iris-Setosa\', 10 2: \'Iris-Versicolor\', 11 3: \'Iris-Virgnica\'} 12 13 # 展示特征用的 14 feature_dict = {0: \'sepal length [cm]\', 15 1: \'sepal width [cm]\', 16 2: \'petal length [cm]\', 17 3: \'petal width [cm]\'} 18 19 # 指定绘图区域大小 20 plt.figure(figsize=(8, 6)) 21 for cnt in range(4): 22 # 这里用子图来呈现4个特征 23 plt.subplot(2, 2, cnt+1) 24 for lab in (\'Iris-setosa\', \'Iris-versicolor\', \'Iris-virginica\'): 25 plt.hist(X[y==lab, cnt], 26 label=lab, 27 bins=10, 28 alpha=0.3,) 29 plt.xlabel(feature_dict[cnt]) 30 plt.legend(loc=\'upper right\', fancybox=True, fontsize=8) 31 32 plt.tight_layout() 33 plt.show()
上述代码可以生成图15-9的输出。可以看出，有些特征区别能力较强，能把3种花各自呈现出来；有些特征区别能力较弱，部分特征数据样本混杂在一起。
第③步：数据的标准化。一般情况下，在进行训练前，数据经常需要进行标准化处理，可以使用sklearn库中的StandardScaler方法进行标准化处理，代码如下：
图15-9 鸢尾花数据集特征
1 from sklearn.preprocessing import StandardScaler
2 X_std = StandardScaler().fit_transform(X)
第④步：计算协方差矩阵。按照式（15.19）定义的协方差矩阵公式计算：
1 mean_vec = np.mean(X_std, axis=0) 2 cov_mat = (X_std - mean_vec).T.dot((X_std - mean_vec)) / (X_std.shape[0]-1) 3 print(\'协方差矩阵 \n%s\' %cov_mat)
协方差矩阵 [[ 1.00675676 -0.10448539 0.87716999 0.82249094] [-0.10448539 1.00675676 -0.41802325 -0.35310295] [ 0.87716999 -0.41802325 1.00675676 0.96881642] [ 0.82249094 -0.35310295 0.96881642 1.00675676]]
或者可以直接使用Numpy工具包来计算协方差，结果是一样的：
1 print(\'NumPy 计算协方差矩阵: \n%s\' %np.cov(X_std.T))
NumPy 计算协方差矩阵: [[ 1.00675676 -0.10448539 0.87716999 0.82249094] [-0.10448539 1.00675676 -0.41802325 -0.35310295] [ 0.87716999 -0.41802325 1.00675676 0.96881642] [ 0.82249094 -0.35310295 0.96881642 1.00675676]]
第⑤步：求特征值与特征向量。
1 cov_mat = np.cov(X_std.T) 2 3 eig_vals, eig_vecs = np.linalg.eig(cov_mat) 4 5 print(\'特征向量 \n%s\' %eig_vecs) 6 print(\'\n特征值 \n%s\' %eig_vals)
特征向量 [[ 0.52308496 -0.36956962 -0.72154279 0.26301409] [-0.25956935 -0.92681168 0.2411952 -0.12437342] [ 0.58184289 -0.01912775 0.13962963 -0.80099722] [ 0.56609604 -0.06381646 0.63380158 0.52321917]] 特征值 [2.92442837 0.93215233 0.14946373 0.02098259]
第⑥步：按照特征值大小进行排序。
1 # 把特征值和特征向量对应起来 2 eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))] 3 print (eig_pairs) 4 print (\'----------\') 5 # 把它们按照特征值大小进行排序 6 eig_pairs.sort(key=lambda x: x[0], reverse=True) 7 8 # 打印排序结果 9 print(\'特征值又大到小排序结果:\') 10 for i in eig_pairs: 11 print(i[0])
[(2.9244283691111135, array([ 0.52308496, -0.25956935, 0.58184289, 0.56609604])),
(0.9321523302535066, array([-0.36956962, -0.92681168, -0.01912775, -0.06381646])),
(0.1494637348981336, array([-0.72154279, 0.2411952 , 0.13962963, 0.63380158])),
(0.02098259276427038, array([ 0.26301409, -0.12437342, -0.80099722, 0.52321917]))] ---------- 特征值又大到小排序结果: 2.9244283691111135 0.9321523302535066 0.1494637348981336 0.02098259276427038
第⑦步：计算累加贡献率。同样可以用累加的方法，将特征向量累加起来，当其超过一定百分比时，就选择其为降维后的维度大小：
1 # 计算累加结果 2 tot = sum(eig_vals) 3 var_exp = [(i / tot)*100 for i in sorted(eig_vals, reverse=True)] 4 print (var_exp) 5 cum_var_exp = np.cumsum(var_exp) 6 cum_var_exp
[72.6200333269203, 23.14740685864414, 3.7115155645845284, 0.5210442498510098]
array([ 72.62003333, 95.76744019, 99.47895575, 100. ])
可以发现，使用前两个特征值时，其对应的累计贡献率已经超过95%，所以选择降到二维。也可以通过画图的形式，这样更直接：
1 a = np.array([1,2,3,4]) 2 print (a) 3 print (\'-----------\') 4 print (np.cumsum(a)) 5 6 plt.figure(figsize=(6, 4)) 7 8 plt.bar(range(4), var_exp, alpha=0.5, align=\'center\', 9 label=\'individual explained variance\') 10 plt.step(range(4), cum_var_exp, where=\'mid\', 11 label=\'cumulative explained variance\') 12 plt.ylabel(\'Explained variance ratio\') 13 plt.xlabel(\'Principal components\') 14 plt.legend(loc=\'best\') 15 plt.tight_layout() 16 plt.show()
上述代码可以生成图15-10的输出。
图15-10 累加特征值
第⑧步：完成PCA降维。接下来把特征向量组合起来完成降维工作：
1 matrix_w = np.hstack((eig_pairs[0][1].reshape(4,1), 2 eig_pairs[1][1].reshape(4,1))) 3 4 print(\'Matrix W:\n\', matrix_w)
Matrix W: [[ 0.52308496 -0.36956962] [-0.25956935 -0.92681168] [ 0.58184289 -0.01912775] [ 0.56609604 -0.06381646]]
1 Y = X_std.dot(matrix_w) 2 Y
array([[-2.10795032, 0.64427554], [-2.38797131, 0.30583307], [-2.32487909, 0.56292316], [-2.40508635, -0.687591 ], [-2.08320351, -1.53025171], [-2.4636848 , -0.08795413], [-2.25174963, -0.25964365], [-2.3645813 , 1.08255676], [-2.20946338, 0.43707676], [-2.17862017, -1.08221046], [-2.34525657, -0.17122946], [-2.24590315, 0.6974389 ], [-2.66214582, 0.92447316], [-2.2050227 , -1.90150522], [-2.25993023, -2.73492274], [-2.21591283, -1.52588897], [-2.20705382, -0.52623535], [-1.9077081 , -1.4415791 ], [-2.35411558, -1.17088308], [-1.93202643, -0.44083479], [-2.21942518, -0.96477499], [-2.79116421, -0.50421849], [-1.83814105, -0.11729122], [-2.24572458, -0.17450151], [-1.97825353, 0.59734172], [-2.06935091, -0.27755619], [-2.18514506, -0.56366755], [-2.15824269, -0.34805785], [-2.28843932, 0.30256102], [-2.16501749, 0.47232759], [-1.8491597 , -0.45547527], [-2.62023392, -1.84237072], [-2.44885384, -2.1984673 ], [-2.20946338, 0.43707676], [-2.23112223, 0.17266644], [-2.06147331, -0.6957435 ], [-2.20946338, 0.43707676], [-2.45783833, 0.86912843], [-2.1884075 , -0.30439609], [-2.30357329, -0.48039222], [-1.89932763, 2.31759817], [-2.57799771, 0.4400904 ], [-1.98020921, -0.50889705], [-2.14679556, -1.18365675], [-2.09668176, 0.68061705], [-2.39554894, -1.16356284], [-2.41813611, 0.34949483], [-2.24196231, -1.03745802], [-2.22484727, -0.04403395], [ 1.09225538, -0.86148748], [ 0.72045861, -0.59920238], [ 1.2299583 , -0.61280832], [ 0.37598859, 1.756516 ], [ 1.05729685, 0.21303055], [ 0.36816104, 0.58896262], [ 0.73800214, -0.77956125], [-0.52021731, 1.84337921], [ 0.9113379 , -0.02941906], [-0.01292322, 1.02537703], [-0.15020174, 2.65452146], [ 0.42437533, 0.05686991], [ 0.52894687, 1.77250558], [ 0.70241525, 0.18484154], [-0.05385675, 0.42901221], [ 0.86277668, -0.50943908], [ 0.33388091, 0.18785518], [ 0.13504146, 0.7883247 ], [ 1.19457128, 1.63549265], [ 0.13677262, 1.30063807], [ 0.72711201, -0.40394501], [ 0.45564294, 0.41540628], [ 1.21038365, 0.94282042], [ 0.61327355, 0.4161824 ], [ 0.68512164, 0.06335788], [ 0.85951424, -0.25016762], [ 1.23906722, 0.08500278], [ 1.34575245, -0.32669695], [ 0.64732915, 0.22336443], [-0.06728496, 1.05414028], [ 0.10033285, 1.56100021], [-0.00745518, 1.57050182], [ 0.2179082 , 0.77368423], [ 1.04116321, 0.63744742], [ 0.20719664, 0.27736006], [ 0.42154138, -0.85764157], [ 1.03691937, -0.52112206], [ 1.015435 , 1.39413373], [ 0.0519502 , 0.20903977], [ 0.25582921, 1.32747797], [ 0.25384813, 1.11700714], [ 0.60915822, -0.02858679], [ 0.31116522, 0.98711256], [-0.39679548, 2.01314578], [ 0.26536661, 0.85150613], [ 0.07385897, 0.17160757], [ 0.20854936, 0.37771566], [ 0.55843737, 0.15286277], [-0.47853403, 1.53421644], [ 0.23545172, 0.59332536], [ 1.8408037 , -0.86943848], [ 1.13831104, 0.70171953], [ 2.19615974, -0.54916658], [ 1.42613827, 0.05187679], [ 1.8575403 , -0.28797217], [ 2.74511173, -0.78056359], [ 0.34010583, 1.5568955 ], [ 2.29180093, -0.40328242], [ 1.98618025, 0.72876171], [ 2.26382116, -1.91685818], [ 1.35591821, -0.69255356], [ 1.58471851, 0.43102351], [ 1.87342402, -0.41054652], [ 1.23656166, 1.16818977], [ 1.45128483, 0.4451459 ], [ 1.58276283, -0.67521526], [ 1.45956552, -0.25105642], [ 2.43560434, -2.55096977], [ 3.29752602, 0.01266612], [ 1.23377366, 1.71954411], [ 2.03218282, -0.90334021], [ 0.95980311, 0.57047585], [ 2.88717988, -0.38895776], [ 1.31405636, 0.48854962], [ 1.69619746, -1.01153249], [ 1.94868773, -0.99881497], [ 1.1574572 , 0.31987373], [ 1.007133 , -0.06550254], [ 1.7733922 , 0.19641059], [ 1.85327106, -0.55077372], [ 2.4234788 , -0.2397454 ], [ 2.31353522, -2.62038074], [ 1.84800289, 0.18799967], [ 1.09649923, 0.29708201], [ 1.1812503 , 0.81858241], [ 2.79178861, -0.83668445], [ 1.57340399, -1.07118383], [ 1.33614369, -0.420823 ], [ 0.91061354, -0.01965942], [ 1.84350913, -0.66872729], [ 2.00701161, -0.60663655], [ 1.89319854, -0.68227708], [ 1.13831104, 0.70171953], [ 2.03519535, -0.86076914], [ 1.99464025, -1.04517619], [ 1.85977129, -0.37934387], [ 1.54200377, 0.90808604], [ 1.50925493, -0.26460621], [ 1.3690965 , -1.01583909], [ 0.94680339, 0.02182097]])
输出结果显示，使用PCA降维算法把原数据矩阵从150×4降到150×2。
第⑨步：可视化对比降维前后数据的分布。由于数据具有4个特征，无法在平面图中显示，因此只使用两维特征显示数据，代码如下：
1 plt.figure(figsize=(6, 4)) 2 for lab, col in zip((\'Iris-setosa\', \'Iris-versicolor\', \'Iris-virginica\'), 3 (\'blue\', \'red\', \'green\')): 4 plt.scatter(X[y==lab, 0], 5 X[y==lab, 1], 6 label=lab, 7 c=col) 8 plt.xlabel(\'sepal_len\') 9 plt.ylabel(\'sepal_wid\') 10 plt.legend(loc=\'best\') 11 plt.tight_layout() 12 plt.show()
上面代码只使用前两个特征显示3类数据，如图15-11所示，看起来3类鸢尾花相互交叠在一起，不容易区分开。
图15-11 原始数据集数据样本分布
下面看看使用PCA降维后的情况，代码如下：
1 plt.figure(figsize=(6, 4)) 2 for lab, col in zip((\'Iris-setosa\', \'Iris-versicolor\', \'Iris-virginica\'), 3 (\'blue\', \'red\', \'green\')): 4 plt.scatter(Y[y==lab, 0], 5 Y[y==lab, 1], 6 label=lab, 7 c=col) 8 plt.xlabel(\'Principal Component 1\') 9 plt.ylabel(\'Principal Component 2\') 10 plt.legend(loc=\'lower center\') 11 plt.tight_layout() 12 plt.show()
上述代码使用降维以后的二维特征作为x,y轴，显示如图15-12所示，对比这两个结果，可以看出经过PCA降维后的结果更容易区别。
图15-12 PCA降维结果
本章小结：
本章介绍了两种非常实用的降维方法：线性判别分析和主成分分析。其中线性判别分析是有监督算法，需要有标签才能计算；而主成分析是无监督算法，无须标签，直接就可以对数据进行分析。那么，在实际建模任务中，到底用哪种降维方法呢？没有固定的模式，需要大家通过实验对比确定，例如，取得一份数据后可以分多步走，以对比不同策略的结果。
降维可以大大减少算法的计算量，加快程序执行的速度，遇到特征非常多的数据时就可以大显身手。但是降维算法本身最大的问题就是降维后得到结果的物理含义很难解释，只能说这就是计算机看来最好的降维特征，而不能具体化其含义。此时如果想进一步对结果进行分析就有些麻烦，因为其中每一个特征指标的含义都只是数值，而没有具体指代。
Python案例中使用非常简单的鸢尾花数据集，按照原理推导的流程一步步完成了整个任务，大家在练习的时候也可以选用稍微复杂一点的数据集来复现算法。
第15章完。
该书资源下载，请至异步社区：https://www.epubit.com