2.2 高维数据降维

高维数据降维是指采用某种映射方法，降低随机变量的数量，例如将数据点从高维空间映射到低维空间中，从而实现维度减少。降维分为特征选择和特征提取两类，前者是从含有冗余信息以及噪声信息的数据中找出主要变量，后者是去掉原来数据，生成新的变量，可以寻找数据内部的本质结构特征。

降维的过程是通过对输入的原始数据特征进行学习，得到一个映射函数，实现将输入样本映射后到低维空间中之后，原始数据的特征并没有明显损失，通常情况下新空间的维度要小于原空间的维度。目前大部分降维算法是处理向量形式的数据。

2.2.1 主成分分析

主成分分析（Principal Component Analysis，PCA）是最常用的线性降维方法，它的目标是通过某种线性投影，将高维的数据映射到低维的空间中，并期望在所投影的维度上数据的方差最大，以此使用较少的维度，同时保留较多原数据的维度。

主成分分析的降维是指经过正交变换后，形成新的特征集合，然后从中选择比较重要的一部分子特征集合，从而实现降维。这种方式并非是在原始特征中选择，所以PCA这种线性降维方式最大程度保留了原有的样本特征。

设有m条n维数据，PCA的一般步骤如下。

① 将原始数据按列组成n行m列矩阵X；

② 计算矩阵X中每个特征属性（n维）的平均向量M（平均值）；

③ 将X的每一行（代表一个属性字段）进行零均值化，即减去M；

④ 按照公式求出协方差矩阵；

⑤ 求出协方差矩阵的特征值及对应的特征向量；

⑥ 将特征向量按对应特征值从大到小按行排列成矩阵，取前k(k<n)行组成基向量P；

⑦ 通过Y=PX计算降维到k维后的样本特征。

PCA算法目标是求出样本数据的协方差矩阵的特征值和特征向量，而协方差矩阵的特征向量的方向就是PCA需要投影的方向。使样本数据向低维投影后，能尽可能表征原始的数据。协方差矩阵可以用散布矩阵代替，协方差矩阵乘以(n-1)就是散布矩阵，n为样本的数量。协方差矩阵和散布矩阵都是对称矩阵，主对角线是各个随机变量（各个维度）的方差。

【例2.3】基于sklearn（Python语言下的机器学习库）和numpy随机生成2个类别共40个三维空间的样本点，生成的代码如下：

　　mu_vec1 = np.array([0,0,0])
　　cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]])
　　class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T
　　mu_vec2 = np.array([1,1,1])
　　cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
　　class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T

其中multivariate_normal()方法生成多元正态分布样本，参数mu_vec1是设定的样本均值向量，cov_mat1是指定的协方差矩阵，每个类别数量为20个。

生成的两个类别class1_sample和class2_sample的样本数据维度为三维，即样本数据的特征数量为3个，将其置于三维空间中展示，可视化结果如图2-7所示。

其中每20个点作为一个类别，平均分成class1和class2两个类别，可以看到三角形和圆点在空间中的分布并没有明显分离。用PCA技术将其投射到二维的空间中，查看其分布情况。

计算40个点在3个维度上的平均向量，首先将两个类别的数据合并到all_samples变量中，然后计算平均向量，代码如下：

　　all_samples = np.concatenate((class1_sample, class2_sample), axis=1)
　　mean_x = np.mean(all_samples[0,:])
　　mean_y = np.mean(all_samples[1,:])
　　mean_z = np.mean(all_samples[2,:])

计算平均向量mean_x、mean_y、mean_z的结果分别为：0.41667492、0.69848315、0.49242335，基于平均向量计算散布矩阵，计算方法如下所示，其中m就是之前计算的平均向量。

所有向量与m的差值经过点积并求和后即可获得散布矩阵的值，代码如下：

　　scatter_matrix = np.zeros((3,3))
　　for i in range(all_samples.shape[1]):
　　     scatter_matrix += (all_samples[:,i].reshape(3,1) - mean_vector).dot((all_samples[:,i].
　　     reshape(3,1) - mean_vector).T)

计算后scatter_matrix的结果为：

　　[[ 38.4878051   10.50787213  11.13746016]
　　 [ 10.50787213  36.23651274  11.96598642]
　　[ 11.13746016  11.96598642  49.73596619]]

应用Python的numpy库内置np.linalg.eig(scatter_matrix)方法计算特征向量和特征值。此外，也可以使用协方差矩阵求解（可用numpy.cov()方法计算协方差矩阵）。代码如下：

　　eig_val_sc, eig_vec_sc = np.linalg.eig(scatter_matrix)

计算出的3个维度的特征值（eig_val_sc）结果分别为：65.16936779、32.69471296、26.59620328，3个维度的特征向量（eig_vec_sc）结果分别为：

　　[-0.49210223  -0.64670286  0.58276136]
　　[-0.47927902  -0.35756937  -0.8015209 ]
　　[-0.72672348  0.67373552  0.13399043]

以平均向量为起点，在图2-8中绘出特征向量，可以看到特征向量的方向，这个方向确定了要进行转换的新特征空间的坐标系。

按照特征值和特征向量进行配对，并按照特征值的大小从高到低进行排序，由于需要将三维空间投射到二维空间中，选择前两个特征值-特征向量对作为坐标，并构建2×3的特征向量矩阵W。原来空间的样本通过与此矩阵相乘，使用公式：y=WTx的方法将所有样本转换到新的空间中。并将结果可视化，如图2-9所示。

可以看到两种类别的样本比三维空间区分度更大，从PCA的实现原理来看，这种变换并没有改变各样本之间的关系，只是应用了新的坐标系。在本例中是将三维空间转降到二维空间，如果有一个n维的数据，想要降低到k维。那么就取前k个特征值对应的特征向量即可。PCA的主要缺点是当数据量和数据维数非常大的时候，用协方差矩阵的方法解PCA会变得非常低效，解决办法是采用奇异值分解（Singular Value Decomposition，SVD）技术。

2.2.2 奇异值分解

对于任意m×n的输入矩阵A，SVD分解结果为：

分解结果中U为左奇异矩阵，S为奇异值矩阵，除主对角线上的元素外全为0，主对角线上的每个元素都称为奇异值，V为右奇异矩阵。矩阵U、V中的列向量均为正交单位向量，而矩阵S为对角阵，并且从左上到右下以递减的顺序排序，可以直接借用SVD的结果来获取协方差矩阵的特征向量和特征值。

【例2.4】基于奇异值分解对Iris数据集降维。

Python语言的numpy库中已实现SVD方法，位于linalg模块（包含核心线性代数工具，如计算逆矩阵、求特征值、解线性方程组以及求解行列式等）中。首先，引入相应模块包，并使用pandas加载Iris数据集（文本格式），代码如下：

　　     from numpy import *
　　import matplotlib.pyplot as plt
　　import pandas as pd
　　from numpy.linalg import *
　　df = pd.read_csv('iris.data.csv', header=None)
　　df[4] = df[4].map({'Iris-setosa':0, 'Iris-versicolor':1, 'Iris-virginica':2})
　　print(df.head())

输出Iris数据集中前5条样本，如图2-10所示，0～3列表示鸢尾花的花萼长度、花萼宽度、花瓣长度和花瓣宽度4个属性，第4列是将鸢尾花的文本分类转换数字格式后的结果。

然后，将Iris数据集的分类标签列排除，使用前4列数据作为linalg.svd()方法的输入进行计算，得到左奇异矩阵U、奇异值矩阵S、右奇异值矩阵V。选择U中前2个特征分别作为二维平面的x、y坐标进行可视化，代码如下：

　　data = df.iloc[:, :-1].values
　　samples,features = df.shape
　　U, s, V = linalg.svd( data )
　　newdata = U[:,:2]
　　fig = plt.figure()
　　ax = fig.add_subplot(1,1,1)
　　colors = ['o','^','+']
　　for i in range(samples):
　　     ax.scatter(newdata[i,0],newdata[i,1],c='black',marker=marks[int(data[i,-1])])
　　plt.xlabel('SVD1')
　　plt.ylabel('SVD2')
　　plt.show()

鸢尾花的3个类别分别用不同的形状表示，结果如图2-11所示。

2.2.3 线性判别分析

线性判别分析（Linear Discriminant Analysis，LDA）也称为Fisher Linear Discriminant，是一种有监督的线性降维算法。与PCA不同，LDA是为了使降维后的数据点尽可能容易地被区分。

线性判别分析在训练过程中，通过将训练样本投影到低维度上，使得同类别的投影点尽可能接近，异类别样本的投影点尽可能远离，即同类点的方差尽可能小，而类之间的方差尽可能大；对新样本，将其投影到低维空间，根据投影点的位置来确定其类别。PCA主要是从特征的协方差角度，去找到比较好的投影方式。LDA更多地考虑了标注，即希望投影后不同类别之间数据点的距离更大，同一类别的数据点更紧凑，如图2-12所示。

计算每一项观测结果的判别分值，对其所处的目标变量所属类别进行判断。这些分值是通过寻找自变量的线性组合得到的。假设每类中的观测结果来自于一个多变量高斯分布，而预测变量的协方差在响应变量y的所有k级别都是通用的。

LDA的降维过程如下：

① 计算数据集中每个类别下所有样本的均值向量；

② 通过均值向量，计算类间散布矩阵SB和类内散布矩阵SW；

③ 依据公式进行特征值求解，计算的特征向量和特征值；

④ 按照特征值排序，选择前k个特征向量构成投影矩阵U；

⑤ 通过Y=X×U的特征值矩阵将所有样本转换到新的子空间中；

LDA在求解过程中需要类内散度矩阵SW和类间散度矩阵SB，其中SW由两类扩展得到，而SB的定义则与两类有所不同，是由每类的均值和总体均值的乘积矩阵求和得到的。目标是求得一个矩阵U使得投影后类内散度尽量小，而类间散度尽量大。在多类情况下，散度表示为一个矩阵。一般情况下，LDA之前会做一次PCA，保证SW矩阵的正定性。

PCA降维是直接与数据维度相关的，例如原始数据是n维，那么使用PCA后，可以任意选最佳的k（k<n）维。LDA降维是与类别个数相关的，与数据本身的维度没关系，例如原始数据是n维的，一共有C个类别，那么LDA降维之后，可选的一般不超过C-1维。例如，假设图像分类，有两个类别为正例和反例，每个图像有1024维特征，那么LDA降维之后，就只有1维特征，而PCA可以选择降到100维。

【例2.5】应用LDA对鸢尾花（Iris）的样本数据进行分析，鸢尾花数据集是20世纪30年代的经典数据集，它由Fisher收集整理，数据集包含150个数据集，分为3类，每类50个数据，每个数据包含4个属性。可通过花萼长度、花萼宽度、花瓣长度和花瓣宽度4个属性预测鸢尾花卉属于山鸢尾（Iris Setosa）、杂色鸢尾（Iris Versicolour）、维吉尼亚鸢尾（Iris Virginica）中的哪种类别，将类别文字转化为数字类别，如表2-2所示。

数据集中有4个特征：萼片长、萼片宽、花瓣长和花瓣宽，总共150行，每一行是一个样本，这就构成了一个4×150的输入矩阵，输出是1列，即花的类别，构成了1×150的矩阵。分析的目标就是通过LDA算法将输入矩阵映射到低维空间中进行分类。

首先计算数据集的平均向量，即计算每种类别下各输入特征的平均值，代码如下：

　　class_labels = np.unique(y)
　　n_classes = class_labels.shape[0]
　　mean_vectors = []
　　for cl in class_labels:
　　     mean_vectors.append(np.mean(X[y==cl], axis=0))

3个类别的平均向量计算结果存于数组变量mean_vectors中，如下所列：

　　[ 5.006,  3.418,  1.464,  0.244],
　　[ 5.936,  2.77 ,  4.26 ,  1.326],
　　 [ 6.588,  2.974,  5.552,  2.026]

计算数据集的类内散布矩阵，其计算方法如下，其中m是上一步中得到的平均向量，xi是每一类别下的样本特征数值，c为类别数量。

计算类内散布矩阵SW的代码如下，其中sc_matrix_class是每个类别下的散布向量值，使用row.reshape()将其转化为列向量格式，按照类别计算后，将结果合并到S_W中。

　　S_W = np.zeros((4,4))
　　for cl,mv in zip(range(1,4), mean_vectors):
　　     sc_matrix_class = np.zeros((4,4))
　　     for row in X[y == cl]:
　　          row, mv = row.reshape(4,1), mv.reshape(4,1)
　　          sc_matrix_class += (row-mv).dot((row-mv).T)
　　     S_W += sc_matrix_class

当然，也可以计算不同类别的协方差矩阵，方法与散布矩阵相似，只是再将结果除以n-1即可。计算后得到类内散布矩阵结果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是上一步中得到的平均向量，xi是每一类别下的样本特征数值，c为类别数量，Ni是类别i的样本数量。

计算类间散布矩阵SB的代码如下，其中mean_vectors是在上一步中计算的均值向量，使用row.reshape()将其转化为列向量格式，n_features是输入变量特征数，本例鸢尾花特征数为4，按照类别计算后，将结果合并到S_B（4×4矩阵）中。

　　overall_mean = np.mean(X, axis=0)
　　n_features = X.shape[1]
　　S_B = np.zeros((n_features, n_features))
　　for i, mean_vec in enumerate(mean_vectors):
　　    n = X[y==i+1,:].shape[0]
　　    mean_vec = mean_vec.reshape(n_features, 1)            #转为列向量形式
　　    overall_mean = overall_mean.reshape(n_features, 1)    #转为列向量
　　    S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)
　　return 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]]

计算特征向量和特征值，应用Python的numpy库内置np.linalg.eig(scatter_matrix)方法计算：

　　eig_vals, eig_vecs = np。linalg.eig(np.linalg.inv(S_W).dot(S_B))

特征值和特征向量组成键值对，并按照由大到小进行排序，结果如下：

　　[(20.904622926374312, array([-0.20673448,-0.415927,0.56155039,0.68478226])),
　　(0.14283325667561703, array([-0.00176467,0.56263241,-0.22318422,0.79600908])),
　　(2.4119371059245178e-15, array([ 0.57400084,0.06633374,0.13588246,-0.80477253])),
　　(5.263703535421987e-16, array([-0.50588695,0.44445486,0.48663891,-0.55652568]))]

可以看到其中只有第一个特征值较大，其他几个接近0，说明只需要一个维度，向第一个特征向量上投射也完全可以。应用前两个特征向量和前一个特征向量的效果如图2-13所示。

其中，如图2-13（a）所示的是使用前2个特征向量构建了一个4×2矩阵转换到子空间中的效果，三种类别的鸢尾花基本上可以完整分离，而如图2-13（b）所示的是只使用第1个特征向量转换到一维空间中的效果，三种类别的鸢尾花也都分隔明显。

2.2.4 局部线性嵌入

流形学习（Manifold Learning）是机器学习中的一种维数约简方法，将高维的数据映射到低维，并依然能够反映原高维数据的本质结构特征。流形学习的前提是假设某些高维数据实际是一种低维的流形结构嵌入在高维空间中。流形学习分为线性流形算法和非线性流形算法，线性流形算法包括主成分分析和线性判别分析，非线性流形算法包括局部线性嵌入（Locally Linear Embedding，LLE）、拉普拉斯特征映射（LE）等。

局部线性嵌入是一种典型的非线性降维算法，这一算法要求每一个数据点都可以由其近邻点的线性加权组合构造得到，从而使降维后的数据也能基本保持原有流形结构。它是流形学习方法最经典的工作之一，后续的很多流形学习、降维方法都与其有密切联系。

局部线性嵌入寻求数据的低维投影，保留本地邻域内的距离。它可以被认为是一系列局部主成分分析，被全局比较以找到最佳的非线性嵌入。

算法的主要步骤分为三步：首先寻找每个样本点的k个近邻点；然后，由每个样本点的近邻点计算出该样本点的局部重建权值矩阵；最后，由该样本点的局部重建权值矩阵和近邻点计算出该样本点的输出值。

LLE在有些情况下也并不适用，例如数据分布在整个封闭的球面上，LLE则不能将它映射到二维空间，且不能保持原有的数据流形。因此在处理数据时，需要确保数据不是分布在闭合的球面或者椭球面上。

【例2.6】用LLE对“瑞士卷”数据集进行降维。下面是基于sklearn实现的LLE示例，数据集生成和LLE降维的核心代码如下：

　　from sklearn import manifold, datasets
　　X, color = datasets.samples_generator.make_swiss_roll(n_samples=1500)
　　X_r, err = manifold.locally_linear_embedding(X, n_neighbors=10, n_components=2)

其中dataset中内置了生成此类数据的方法make_swiss_roll，生成样本数为1500个，X为生成的结果，作为LLE的输入。locally_linear_embedding方法执行LLE运算，n_neighbors=10表示近邻点的数量为10，n_components=2表示降维到二维。

将生成的X_r结果和原始X数据进行可视化显示，核心代码如下，生成的对比效果如图2-14所示，其中如图2-14（a）所示的为原始数据，从侧面看像“瑞士卷”一样，如图2-14（b）所示为投射后样本的可视化效果。

　　ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=color, cmap=plt.cm.Spectral) #原始数据可视化
　　ax.scatter(X_r[:, 0], X_r[:, 1], c=color, cmap=plt.cm.Spectral) #投射数据可视化

可以看到经过LLE变换后，样本数据在低维空间上已经明显区分出来。

2.2.5 拉普拉斯特征映射

拉普拉斯特征映射（Laplacian Eigenmaps，LE）解决问题的思路和LLE相似，是一种基于图的降维算法，使相互关联的点在降维后的空间中尽可能地靠近。

通过构建邻接矩阵为W的图来重构数据流形的局部结构特征，如果两个数据实例i和j很相似，那么i和j在降维后目标子空间中也应该接近。设数据实例的数目为n，目标子空间（即降维后的维度）为m，定义n×m大小的矩阵Y，其中每一个行向量yi是数据实例i在目标子空间中的向量表示。为了让样本i和j在降维后的子空间里尽量接近，优化的目标函数如下：

其中，‖yi−yj‖2为两个样本在目标子空间中的距离，wij是两个样本的权重值，权重值可以用图中样本间的连接数来度量。经过推导，将目标函数转化为以下形式：

Ly=λDy

其中L和D均为对称矩阵，由于目标函数是求最小值，所以通过求得m个最小非零特征值所对应的特征向量，即可达到降维的目的。

拉普拉斯特征映射的具体步骤如下。

① 构建无向图，将所有的样本以点连接成一个图，例如使用KNN算法，将每个点最近的k个点进行连接，其中k是一个预先设定的值。

② 构建图的权值矩阵，通过点之间的关联程度来确定点与点之间的权重大小，例如，两个点之间如果相连接，则权重值为1，否则为0。

③ 特征映射，通过公式Ly=λDy计算拉普拉斯矩阵L的特征向量与特征值，用最小的m个非零特征值对应的特征向量作为降维的结果。

【例2.7】使用拉普拉斯特征映射的方法对“瑞士卷”数据集进行降维。基于sklearn实现“瑞士卷”样本生成和降维的核心代码如下：

　　X, color = datasets.samples_generator.make_swiss_roll(n_samples=1500)
　　#原始数据可视化
　　ax = fig.add_subplot(121, projection='3d')
　　ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=color, cmap=plt.cm.Spectral)
　　ax.view_init(4, -72)

首先，通过sklearn库中的make_swiss_roll()方法生成“瑞士卷”形状的测试样本，并用matplotlib库中的scatter方法以三维形式展示，同时用view_init方法调整三维视角，可以明显看到数据的形态，如图2-15所示。

对原始样本进行LE降维，代码如下：

　　se = manifold.SpectralEmbedding(n_components=2, n_neighbors=10)
　　Y = se.fit_transform(X)
　　#LE降维后可视化
　　ax = fig.add_subplot(122)
　　plt.scatter(Y[:, 0], Y[:, 1], c=color, cmap=plt.cm.Spectral)
　　ax.xaxis.set_major_formatter(NullFormatter())
　　ax.yaxis.set_major_formatter(NullFormatter())

其中，sklearn库中拉普拉斯特征映射是在SpectralEmbedding类中实现的，其参数n_components和n_neighbors分别表示目标子空间的特征数量（维度）和构建无向图时的最近邻数量。初始化se对象之后调用fit_transform()方法将原始样本投射到新的子空间中，效果如图2-16所示。

对比原始样本在三维空间和降维之后在二维空间的分布情况，可以看到，在高维和低维空间中的样本分布形状发生了变化，但降维后样本之间的关联关系并没有改变。