《机器学习》从零开始学(24) 距离模型专题之降维¶

参考教材：

机器学习-周志华-清华大学出版社-2016  
图解机器学习-[日]杉山将-许永伟(译)-人民邮电出版社-2015

所谓“从零开始”面向的对象其实仅指我本人。也就是我的学习笔记。

注意：此笔记是在我理解基础上的梳理和表述。另外本笔记代码由我额外提供。 为了理解原理，我提供的源码尽可能不用机器学习算法包。

本节重点在 降维，涉及：主成分分析、核化主成分分析。

5）低维嵌入¶

如果我们获得的样本集是高维的，即有众多属性，并且数据样本是稀疏的，那么会导致距离计算的困难，或者距离计算偏离我的预期。比如，两个离得较远的样本之间的直线距离连线附近没有任何其它样本，如图箭头所示的三维距离：

高维样本集的样本稀疏、距离计算困难等问题，被称为"维数灾难"。而降维是解决这个问题重要途径，即通过某种数学变换将原始的高维属性空间变成一个低维子空间，样本之间的真实距离（测地线距离，不一定是直线距离）在降维前后保持不变。这种低维子空间称作是原高维空间的低维嵌入。  如图：

我们的目标是：将原始数据集$\boldsymbol{X} = (\boldsymbol{x}_1^T;\dots;\boldsymbol{x}_m^T)_{m \times n}$，降维成$\boldsymbol{Z} = (\boldsymbol{z}_1^T;\dots;\boldsymbol{z}_m^T)_{m \times d}$ 。

下面这个算法，基于要求原始空间中的样本之间的距离在低维空间中得以保持，被称作多维缩放（Multiple Dimensional Scaling, MDS），是一种经典的降维方法。

书中是从保距离角度开始论述， 但实际上保内积肯定就保距离(反之不正确)，所以我们只需对原始数据集$\boldsymbol{X}$中心化后(后面假设$\boldsymbol{X}$已经中心化)，直接算内积矩阵，期望也是$\boldsymbol{Z}$的内积矩阵，进而降维后保距离:

$$\boldsymbol{B} = \boldsymbol{X} \boldsymbol{X}^T$$

然后对$\boldsymbol{B}$做特征值*， 内积矩阵是对称矩阵*后的左右特征矢量矩阵式相等的，所以*后可表示成：

$$\boldsymbol{B} = \boldsymbol{V}^* \boldsymbol{\Lambda}^* \boldsymbol{V}^{*T}$$

其中，$\boldsymbol{\Lambda}^*=\mathrm{diag}(\lambda_1,\dots, \lambda_m)$是特征值构成的对角矩阵，$\boldsymbol{V}_{m \times m}^*=(\boldsymbol{v}_1,\dots,\boldsymbol{v}_n)$是特征向量矩阵。如果我们希望维降到$d$维，我们只需要取最大的$d$个特征值，及其对应的特征向量(假设特征值已经从大到小排好序)：

$$\boldsymbol{\lambda} = (\lambda_1,\dots, \lambda_d),\quad \boldsymbol{\Lambda} = \mathrm{diag}(\boldsymbol{\lambda}) , \quad   \boldsymbol{V} =  (\boldsymbol{v}_1,\dots,\boldsymbol{v}_d)$$

进而获得降维后的样本集：

$$\boldsymbol{Z} = \boldsymbol{V}  \boldsymbol{\Lambda}^{1/2} \in \mathbb{R}^{m \times d}$$

注意:降维后不一定有效果，其效果是需要评估的。 如果性能有所提高，则认为是有效果的。

【算法流程】：

首先，根据给定样本集$\boldsymbol{X}_{m \times n}$，然后计算中心化后的内积矩阵$\boldsymbol{B} = \boldsymbol{X} \boldsymbol{X}^T$；

然后，对$\boldsymbol{B}$进行特征值*获得特征值$\boldsymbol{\lambda}_{m \times 1}^*$和特征向量矩阵$\boldsymbol{V}_{m \times m}^*$， 取最大的$d$个特征值$\boldsymbol{\lambda}_{d \times 1}$及其对应特征向量$\boldsymbol{V}_{m \times d}$;

于是，计算出降维后的数据集$\boldsymbol{Z}_{m \times d} = \boldsymbol{V}  \boldsymbol{\Lambda}^{1/2}$：

最后，在降维后的数据集$\boldsymbol{Z}$进行机器学习算法。

下面是这种降维方法的Python实现：

import numpy as npfrom scipy import linalgimport matplotlib.pyplot as pltimport mpl_toolkits.mplot3d.axes3d as p3from sklearn.cluster import KMeansfrom sklearn.datasets.samples_generator import make_swiss_roll# 用sklearn生成瑞士蛋卷数据集n_samples = 1500noise = 0.05X, _ = make_swiss_roll(n_samples, noise)X = X[:,:3] # 截取我们所需#X[:, 1] *= .5 # 使其更瘦# 绘制数据fig = plt.figure()ax = p3.Axes3D(fig)ax.view_init(7, -80)ax.plot(X[:, 0], X[:, 1], X[:, 2], 'o', markerfacecolor='b', markeredgecolor='k', markersize=6)plt.title(u"原始未分类的样本数据")plt.show()# 数据集中心化(不影响距离矩阵，但会影响内积矩阵，但没关系)_X = X - X.mean(axis = 0)# 生成中心化后的内积矩阵(保内积矩阵 => 保了距离矩阵)B = np.dot(_X, _X.T)# 对矩阵B做特征值*lam_star, V_star =linalg.eig(B)# 在给定精度下近似处理特征值矢量lam_star[np.abs(lam_star) < 0.001] = 0if (np.abs(lam_star.imag)< 0.001).prod():lam_star = lam_star.realV_star[np.abs(V_star) < 0.001] = 0# 对lam从大到小排序(索引)idxLam = np.argsort(-lam_star)print u"特征值矢量(最大10个): ",list(lam_star[idxLam])[:10]# 降成2维，取最大的两个特征值lam = lam_star[idxLam[:2]]V = V_star[:,idxLam[:2]]if (np.abs(V.imag)< 0.001).prod():V = V.real# 最后获得降维后的数据集Z = np.dot(V, np.diag(np.sqrt(lam)))m,n = Z.shape# 这里我们的目标是降维。至于降维后的聚类，不妨直接调用机器学习算法包（采用k均值聚类）y_pred = KMeans(n_clusters=9, random_state=170).fit_predict(Z)# 最后绘制降维后的数据集上的分类plt.scatter(Z[:, 0], Z[:, 1], c=y_pred)plt.title(u"在降维空间中的样本数据的分类")plt.show()# 绘制原始数据集上的分类fig = plt.figure()ax = p3.Axes3D(fig)ax.view_init(7, -80)ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y_pred)plt.title(u"原始样本数据的分类")plt.show()

特征值矢量(最大10个):  [78740.961127822869, 62358.668593818271, 54344.633865828197, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

6）主成分分析¶

注意：书中的样本矩阵：行是属性，列是样本。 但本笔记始终坚持： 行是样本，列是属性。 所以下面和书中有点差异。

最简单的降维是线性降维， 即使用投影矩阵$\boldsymbol{W}_{d \times n}$：

$$\boldsymbol{z}_i = \boldsymbol{W} \boldsymbol{x}_i, \quad i = 1,\dots,m$$

进而得到样本集变换$\boldsymbol{W}_{d \times n}: \boldsymbol{X}_{m \times n} \rightarrow \boldsymbol{Z}_{m \times d}$：

$$\boldsymbol{Z} = \boldsymbol{X} \boldsymbol{W}^T$$

这里的变换矩阵$\boldsymbol{W}=(\boldsymbol{w}_1^T;\dots;\boldsymbol{w}_d^T)$可以看成$d$个$n$维基向量，$\boldsymbol{z}_i$就是第$i$个样本$\boldsymbol{x}_i$与这个组基向量的内积构成$d$维向量。 也就是说新向量$\boldsymbol{z}_i$的分量是老向量$\boldsymbol{x}_i$在新坐标系$\{\boldsymbol{w}_1,\dots,\boldsymbol{w}_d\}$下的坐标。 如果$\boldsymbol{w}_i$与$\boldsymbol{w}_j \quad (i \ne j)$正交，则新坐标系是一个正交坐标系，此时$\boldsymbol{W}$为正交变换。

主成分分析（Principal Component Analysis, PCA）是最常用的降维方法。此方法的基本目标：对一个正交属性空间中的样本点，如何用一个超平面对所有样本进行恰当表达?

下面两种方式可接近这个目标：

• 最近重构性： 样本点到这个超平面都足够近。

• 最大可分性:  样本点在这个超平面上的投影都尽可能分开。

有趣的是，这两种方法得到的主成分分析是相互等价的。

基于最近重构性¶

不妨假设原始样本数据集$\boldsymbol{X}$已经中心化：$\sum_{i=1}^m \boldsymbol{x}_i = 0$；

假设属性空间$\boldsymbol{X}$存在一个标准正交变换矩阵$\boldsymbol{\mathbb{W}}_{n \times n} = (\boldsymbol{W}_{d \times n}; \boldsymbol{W'}_{(n-d) \times n})$，其行子集矩阵$\boldsymbol{W}_{d \times n}$可以把数据集$\boldsymbol{X}_{m \times n}$降维成$\boldsymbol{Z}_{m \times d}$，即： $\|\boldsymbol{w}_i\| = 1 , \quad \boldsymbol{w}_i^T \boldsymbol{w}_j=0, \quad i \ne j , \quad i,j = 1,\dots,d$；

因为$\boldsymbol{\mathbb{W}}_{n \times n}$是正交变换，尽管无法从$\boldsymbol{Z}$复原到$\boldsymbol{X}$，但复原到$\boldsymbol{X}$在$\boldsymbol{Z}$所属$d$维子空间上的是没有问题的(下面第二个式是投影部分)：

$$\boldsymbol{X} = \boldsymbol{Z} \boldsymbol{W} + \boldsymbol{X} \boldsymbol{W'}^T \boldsymbol{W'}$$$$\hat{\boldsymbol{X}} = \boldsymbol{Z} \boldsymbol{W}$$

我们希望原始样本集$\boldsymbol{X}$到投影样本集$\hat{\boldsymbol{X}}$的距离越小越好，即： $$ \begin{array} &\boldsymbol{W}^* &=& \underset{\boldsymbol{W}}{\mathrm{argmin}} \ \mathrm{Tr}\left( (\boldsymbol{Z} \boldsymbol{W} - \boldsymbol{X}) (\boldsymbol{Z} \boldsymbol{W} - \boldsymbol{X})^T \right)\\ \quad &s.t.&  \boldsymbol{W} \boldsymbol{W}^T = \boldsymbol{I} \end{array} $$

其中：

$ \begin{array} & \quad  \mathrm{Tr}\left( (\boldsymbol{Z} \boldsymbol{W} - \boldsymbol{X}) (\boldsymbol{Z} \boldsymbol{W} - \boldsymbol{X})^T \right) \quad {展开剔除常数项并利用迹性质}\\ = \mathrm{Tr}\left( \boldsymbol{Z} \boldsymbol{W} \boldsymbol{W}^T \boldsymbol{Z}^T - 2 \boldsymbol{X}  \boldsymbol{W}^T \boldsymbol{Z}^T \right) \quad {带入 \boldsymbol{X} *式} \\ = - \mathrm{Tr}\left(\boldsymbol{Z} \boldsymbol{W} \boldsymbol{W}^T \boldsymbol{Z}^T + 2 \boldsymbol{X} \boldsymbol{W'}^T \boldsymbol{W'} \boldsymbol{W}^T \boldsymbol{Z}^T \right)  \quad {利用 \boldsymbol{\mathbb{W}} 正交性: \boldsymbol{W} \boldsymbol{W}^T = \boldsymbol{I}, \boldsymbol{W'} \boldsymbol{W}^T = \boldsymbol{0} }\\ = - \mathrm{Tr}\left(\boldsymbol{Z} \boldsymbol{Z}^T \right)  \quad {利用迹性质} \\ = - \mathrm{Tr}\left(\boldsymbol{Z}^T \boldsymbol{Z} \right)  \quad {带入 \boldsymbol{Z} 的投影式} \\ = - \mathrm{Tr}\left(\boldsymbol{W} \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{W}^T \right) \end{array} $

于是获得主成分分析的优化目标：

$$ \begin{array} &\boldsymbol{W}^* &=& \underset{\boldsymbol{W}}{\mathrm{argmin}} \ - \mathrm{Tr}\left(\boldsymbol{W} \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{W}^T \right)\\ \quad &s.t.&  \boldsymbol{W} \boldsymbol{W}^T = \boldsymbol{I} \end{array} $$

基于最大可分性¶

最大可分性要求： 样本点在这个超平面上的投影都尽可能分开, 即投影后的方差$\mathrm{Tr}\left(\boldsymbol{Z}^T \boldsymbol{Z} \right)$最大。 我们注意到这个表达式在前一小段的推导中出现过，所以必然也推导出相同的主成分分析的优化目标。

主成分分析优化问题的求解¶

利用拉格朗日乘子法（并利用迹的求导性质），容易得到：

$$ \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{W}^T = \lambda \boldsymbol{W}^T $$

注意区别：$\boldsymbol{X}^T \boldsymbol{X}$是$n \times n$矩阵，是$n$个属性间协方差矩阵；而$\boldsymbol{X} \boldsymbol{X}^T$是$m \times m$矩阵，是$m$个样本间内积矩阵。

到这里我们发现， 只需要对中心化后的协方差矩阵$\boldsymbol{X}^T \boldsymbol{X}$进行特征值*获得$\boldsymbol{\lambda}^*=(\lambda_1,\dots, \lambda_n)$和特征向量矩阵$\boldsymbol{V}^*=(\boldsymbol{v}_1,\dots,\boldsymbol{v}_n)$。然后，从中取得最大的$d$个特征值，及其对应的特征向量(假设特征向量已经从大到小排好序)：

$$\boldsymbol{\lambda} = (\lambda_1,\dots, \lambda_d),\quad \boldsymbol{\Lambda} = \mathrm{diag}(\boldsymbol{\lambda}) , \quad   \boldsymbol{W} = \boldsymbol{V}^T = (\boldsymbol{v}_1,\dots,\boldsymbol{v}_d)^T$$

进而获得降维后的样本集：

$$\boldsymbol{Z} = \boldsymbol{X} \boldsymbol{V}$$

降维后的维数$d$通常是由用户指定，或者通过交叉验证选取合适的维数$d$。

【算法流程】：

首先，根据给定样本集$\boldsymbol{X}_{m \times n}$，然后计算中心化后的协方差矩阵$\boldsymbol{X}^T \boldsymbol{X}$；

然后，对$\boldsymbol{X}^T \boldsymbol{X}$进行特征值*获得特征值$\boldsymbol{\lambda}_{n \times 1}^*$和特征向量矩阵$\boldsymbol{V}_{n \times n}^*$， 取最大的$d$个特征值$\boldsymbol{\lambda}_{d \times 1}$及其对应特征向量$\boldsymbol{V}_{n \times d}$;

于是，计算出降维后的数据集$\boldsymbol{Z}_{m \times d} = \boldsymbol{X} \boldsymbol{V}$：

最后，在降维后的数据集$\boldsymbol{Z}$进行机器学习算法。

下面是主成分分析的Python实现（明显要比上一节快）：

import numpy as npfrom scipy import linalgimport matplotlib.pyplot as pltimport mpl_toolkits.mplot3d.axes3d as p3from sklearn.cluster import KMeansfrom sklearn.datasets.samples_generator import make_swiss_roll# 用sklearn生成瑞士蛋卷数据集n_samples = 1500noise = 0.05X, _ = make_swiss_roll(n_samples, noise)X = X[:,:3] # 截取我们所需#X[:, 1] *= .5 # 使其更瘦# 绘制数据fig = plt.figure()ax = p3.Axes3D(fig)ax.view_init(7, -80)ax.plot(X[:, 0], X[:, 1], X[:, 2], 'o', markerfacecolor='b', markeredgecolor='k', markersize=6)plt.title(u"原始未分类的样本数据")plt.show()# 数据集中心化_X = X - X.mean(axis = 0)# 生成中心化后的协方差矩阵C = np.dot(_X.T, _X)# 对协方差矩阵C做特征值*lam_star, V_star =linalg.eig(C)# 在给定精度下近似处理特征值矢量lam_star[np.abs(lam_star) < 0.001] = 0if (np.abs(lam_star.imag)< 0.001).prod():lam_star = lam_star.realV_star[np.abs(V_star) < 0.001] = 0# 对lam从大到小排序(索引)idxLam = np.argsort(-lam_star)print u"特征值矢量(最大10个): ",list(lam_star[idxLam])[:10]# 降成2维，取最大的两个特征值lam = lam_star[idxLam[:2]]V = V_star[:,idxLam[:2]]if (np.abs(V.imag)< 0.001).prod():V = V.real# 最后获得降维后的数据集Z = np.dot(X, V)m,n = Z.shape# 这里我们的目标是降维。至于降维后的聚类，不妨直接调用机器学习算法包（采用k均值聚类）y_pred = KMeans(n_clusters=9, random_state=170).fit_predict(Z)# 最后绘制降维后的数据集上的分类plt.scatter(Z[:, 0], Z[:, 1], c=y_pred)plt.title(u"在降维空间中的样本数据的分类")plt.show()# 绘制原始数据集上的分类fig = plt.figure()ax = p3.Axes3D(fig)ax.view_init(7, -80)ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y_pred)plt.title(u"原始样本数据的分类")plt.show()

特征值矢量(最大10个):  [76858.131250895298, 61536.005552176728, 55150.64580363025]

7）核化线性降维¶

线性降维方法是假设从高维空间到低维空间的映射是线性的，但在现实任务中，可能需要非线性映射才能找到恰当的低维嵌入。

可以采用支持向量机的核技巧类似的方法， 将原样本属性空间非线性映射到高维特征空间，然后再进行线性降维，这种方法称作核化线性降维。

假设样本属性空间到高维特征空间的一个非线性映射$\phi:\boldsymbol{x} \rightarrow \boldsymbol{u}$：

$$\phi(\boldsymbol{x}) = \boldsymbol{u}$$$$\phi(\boldsymbol{x}_i) = \boldsymbol{u}_i, \quad i=1,\dots,m$$

于是也把原样本集$\boldsymbol{X}=(\boldsymbol{x}_1^T;\dots;\boldsymbol{x}_m^T)$映射高维特征空间的样本集为$\boldsymbol{U}=(\boldsymbol{u}_1^T;\dots;\boldsymbol{u}_m^T)$：

$$\phi(\boldsymbol{X}) = \boldsymbol{U}$$

于是我们试图在高维特征空间对$\boldsymbol{U}$降维，套用前面主成分分析降维方法有：

$$ \boldsymbol{W} \sum_{i=1}^m \phi(\boldsymbol{x}_i) \phi(\boldsymbol{x}_i)^T = \boldsymbol{W} \sum_{i=1}^m \boldsymbol{u}_i \boldsymbol{u}_i^T = \boldsymbol{W} \boldsymbol{U}^T \boldsymbol{U} = \lambda \boldsymbol{W} $$

于是有：

$$ \boldsymbol{W} = \sum_{i=1}^m \boldsymbol{\alpha}_i \phi(\boldsymbol{x}_i)^T, \quad \boldsymbol{\alpha}_i  = \dfrac{1}{\lambda} \boldsymbol{W} \phi(\boldsymbol{x}_i) $$

一般并不清楚$\phi(.)$的具体形式，于是引入核函数：

$$ \kappa(\boldsymbol{x}_i,\boldsymbol{x}_j) = \phi(\boldsymbol{x}_i)^T \phi(\boldsymbol{x}_j) $$

于是可推导出：

$$ \sum_{j=1}^m \kappa(\boldsymbol{x}_i,\boldsymbol{x}_j) \boldsymbol{\alpha}_j = \lambda \boldsymbol{\alpha}_i, \quad i=1,\dots,m $$

简化为：

$$ \boldsymbol{K} \boldsymbol{A} = \lambda \boldsymbol{A}, \quad \boldsymbol{K} = \left\{\kappa(\boldsymbol{x}_i,\boldsymbol{x}_j)\right\}_{m \times m}, \quad \boldsymbol{A} = (\boldsymbol{\alpha}_1^T;\dots;\boldsymbol{\alpha}_m^T) $$

然后对$\boldsymbol{K}$做特征值*，获得特征值$\boldsymbol{\lambda}^*=(\lambda_1,\dots, \lambda_m)$及其对应的特征向量矩阵$\boldsymbol{V}_{m \times m}^*=(\boldsymbol{v}_1,\dots,\boldsymbol{v}_n)$。如果我们希望维降到$d$维，我们只需要取最大的$d$个特征值，及其对应的特征向量(假设特征值已经从大到小排好序)：

$$\boldsymbol{\lambda} = (\lambda_1,\dots, \lambda_d),\quad \boldsymbol{\Lambda} = \mathrm{diag}(\boldsymbol{\lambda}) , \quad   \boldsymbol{A} = \boldsymbol{V} =  (\boldsymbol{v}_1,\dots,\boldsymbol{v}_d)$$

于是可以将样本集降维成：

$$ \boldsymbol{Z} = \boldsymbol{U} \boldsymbol{W}^T, \quad \boldsymbol{z}_i = \boldsymbol{W} \boldsymbol{u}_i = \boldsymbol{W} \phi(\boldsymbol{x}_i) = \sum_{j=1}^m \boldsymbol{\alpha}_j \kappa(\boldsymbol{x}_j, \boldsymbol{x}_i), \quad i=1,\dots,d $$

进而：

$$ \boldsymbol{Z} = \boldsymbol{K} \boldsymbol{A} $$

特别地， 对新样本$\boldsymbol{x}$则有：

$$ \boldsymbol{z} = \sum_{i=1}^m \boldsymbol{\alpha}_i \kappa(\boldsymbol{x}_i, \boldsymbol{x}), \quad z_j = \sum_{i=1}^m \alpha_{i,j} \kappa(\boldsymbol{x}_i, \boldsymbol{x}), \quad j=1,\dots,d $$

至此，核化线性降维的算法就完整了。

【算法流程】：

首先，根据给定样本集，计算根据指定核函数计算Gram矩阵$\boldsymbol{K}_{m \times m}$；

然后，对Gram矩阵$\boldsymbol{K}$进行特征值*获得特征值$(\lambda_1,\dots,\lambda_m)$和特征矢量$\boldsymbol{A}_{m \times m}^*$， 取最大的$d$个特征值$(\lambda_1,\dots,\lambda_d)$及其对应特征向量$\boldsymbol{A}_{m \times d}$;

于是，计算出降维后的数据集$\boldsymbol{Z}_{m \times d} = \boldsymbol{K} \boldsymbol{A}$ （这步比较耗时）

最后，在降维后的数据集$\boldsymbol{Z}$进行机器学习算法。

下面用Python代码实现这套算法（效果不好，可能核选的不好，可能样本集不适合，以后我会找一个更好的例子）：

import numpy as npfrom scipy import linalgimport matplotlib.pyplot as pltimport mpl_toolkits.mplot3d.axes3d as p3from sklearn.cluster import KMeansfrom sklearn.datasets.samples_generator import make_swiss_rolldef guassKappa(X,Z=None):'''        选择高斯核函数计算Gram矩阵            矩阵X，Z必须有相同列（对应属性），行可以不一样        X是必须训练数据集或其子集        Z可是训练数据集或其子集，也可是测试数据集或其子集    '''    if Z is None:Z = Xassert X.shape[1] == Z.shape[1] 
    # norm_{i,j} = (x_i-z_j)^T (x_i-z_j) = x_i^T x_i - 2 x_i^T z_j + z_j^T z_j# i = 1,...,m1# j = 1,...,m2# m1是X的行数， m2是Z的行数norm = ((X**2).sum(axis=1) - 2*np.dot(Z,X.T)).T + (Z**2).sum(axis=1)# 取 2 * \sigma^2  = X.shape[1]return np.exp(-norm/X.shape[1])# 用sklearn生成瑞士蛋卷数据集n_samples = 1500noise = 0.05X, _ = make_swiss_roll(n_samples, noise)X = X[:,:3] # 截取我们所需#X[:, 1] *= .5 # 使其更瘦# 绘制数据fig = plt.figure()ax = p3.Axes3D(fig)ax.view_init(7, -80)ax.plot(X[:, 0], X[:, 1], X[:, 2], 'o', markerfacecolor='b', markeredgecolor='k', markersize=6)plt.title(u"原始未分类的样本数据")plt.show()# 数据集中心化_X = X - X.mean(axis = 0)# 用高斯核函数计算Gram矩阵kappa = guassKappa(_X)# 对Gram矩阵做特征值*lam_star, V_star =linalg.eig(kappa)# 在给定精度下近似处理特征值矢量lam_star[np.abs(lam_star) < 0.001] = 0if (np.abs(lam_star.imag)< 0.001).prod():lam_star = lam_star.realV_star[np.abs(V_star) < 0.001] = 0# 对lam从大到小排序(索引)idxLam = np.argsort(-lam_star)print u"特征值矢量(最大10个): ",list(lam_star[idxLam])[:10]# 降成2维，取最大的两个特征值lam = lam_star[idxLam[:2]]V = V_star[:,idxLam[:2]]if (np.abs(V.imag)< 0.001).prod():V = V.real# 最后获得降维后的数据集Z = kappa.dot(V).dot(np.diag(np.sqrt(lam)))m,n = Z.shape# 这里我们的目标是降维。至于降维后的聚类，不妨直接调用机器学习算法包（采用k均值聚类）y_pred = KMeans(n_clusters=9, random_state=170).fit_predict(Z)# 最后绘制降维后的数据集上的分类plt.scatter(Z[:, 0], Z[:, 1], c=y_pred)plt.title(u"在降维空间中的样本数据的分类")plt.show()# 绘制原始数据集上的分类fig = plt.figure()ax = p3.Axes3D(fig)ax.view_init(7, -80)ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y_pred)plt.title(u"原始样本数据的分类")plt.show()

特征值矢量(最大10个):  [12.977825088028101, 12.284602029244516, 12.122891787045363, 11.517778685059247, 11.320614914028665, 11.20866092493381, 11.020577139718808, 10.138779442394148, 10.105915184070023, 10.085170928986706]