使用不同降维算法对MNIST数字图片进行降维，并在二维空间中展示出不同类图片在低维空间的布局，可以观察降维后的数据是否利于分类。
本案例主要包括以下内容：
1. 手写数字图片的获取
2. 展示其中的一部分数字图片
3. 降维及可视化
3.1 Random projection
3.2 PCA
3.3 truncated SVD
3.4 LDA
3.5 MDS
3.6 Isomap
3.7 LLE
3.8 t-SNE
3.9 RandomTrees
4. 总结
手写数字识别是数据分析师入门计算机视觉（computer vision）的第一个“大作业”。
绝大多数的竞赛以及相关文章中的手写数字图片数据的来源为MNIST database（Modified National Institute of Standards and Technology database），之所以叫modified，是因为该数据库为超大数据库NIST的一个子集，经过调整之后这些数字图片的分辨率相同，并且将数字置于图片的中心。MNIST database中的训练集有60000个样本，测试集有10000个样本，完整的数据可从杨立昆教授（Yann Lecun）的MNIST主页获得，并且可以从该主页看到许许多多识别算法的尝试。该图片数据集自1999年公开以来，已经成为众多突破性的分类算法的试炼场，从线性分类器（linear classifier）的7.6%到2016年Convolutional neural network（CNN，卷积神经网络）的0.21%，众多机器学习领域的高手孜孜不倦地为减小识别错误率（error rate）而努力。
本案例提供了1083张手写数字图片的降维及可视化过程，具体的手写数字识别模型的训练（从传统的kNN、SVM到深度学习模型）可参考kaggle上的Digit Recognizer竞赛。
0 包版本信息¶
通过在notebook中执行shell命令查看包的版本。
!pip freeze | grep matplotlib;
!pip freeze | grep scikit-learn;
!pip freeze | grep numpy
1 手写数字图片获取¶
在本案例中，我们选择直接从sklearn.datasets模块中通过load_digits导入手写数字图片数据集，该数据集是UCI datasets的Optical Recognition of Handwritten Digits Data Set中的测试集，并且只是MNIST的很小的子集，一共有1797张分辨率为8$\times$8的手写数字图片。同时，该图片有从0到9共十类数字。
我们先导入load_digits模块及本案例所需的相关的包。
from time import time # 用于计算运行时间
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import offsetbox # 定义图形box的格式
from sklearn import (manifold, datasets, decomposition, ensemble,
discriminant_analysis, random_projection)
%matplotlib inline
%config InlineBackend.figure_format = 'retina' # 调整notebook中输出图片的清晰度，使得在retina屏幕上显示得更清楚
load_digits中有n_class参数，可以指定选择提取多少类的图片（从数字0开始），缺省值为10；还有一个return_X_y参数（sklearn 0.18版本的新参数），若该参数值为True，则返回图片数据data和标签target，默认为False。return_X_y为False的情况下，将会返回一个Bunch对象，该对象是一个类似字典的对象，其中包括了数据data、images和数据集的完整描述信息DESCR。下面我们就这两种读取方式分别展示：
- 返回Bunch对象
digits = datasets.load_digits(n_class=6)
print digits
# 获取bunch中的data,target
print digits.data
print digits.target
- 只返回data和target
data = datasets.load_digits(n_class=6)
print data
本案例只提取从数字0到数字5的图片进行降维，我们挑选图片数据集中的前六个照片进行展示：
# plt.gray()
fig, axes = plt.subplots(nrows=1, ncols=6, figsize=(8, 8))
for i,ax in zip(range(6),axes.flatten()):
ax.imshow(digits.images[i], cmap=plt.cm.gray_r)
plt.show()
为了能够方便的展示手写数字图片，我们使用返回Bunch对象的导入方式。
digits = datasets.load_digits(n_class=6)
X = digits.data
y = digits.target
n_samples, n_features = X.shape
n_neighbors = 30
2 展示其中的一部分数字图片¶
n_img_per_row = 30 # 每行显示30个图片
# 整个图形占 300*300，由于一张图片为8*8，所以每张图片周围包了一层白框，防止图片之间互相影响
img = np.zeros((10 * n_img_per_row, 10 * n_img_per_row))
for i in range(n_img_per_row):
ix = 10 * i + 1
for j in range(n_img_per_row):
iy = 10 * j + 1
img[ix:ix + 8, iy:iy + 8] = X[i * n_img_per_row + j].reshape((8, 8))
plt.figure(figsize=(6,6))
plt.imshow(img, cmap=plt.cm.binary)
plt.xticks([])
plt.yticks([])
plt.title('A selection from the 64-dimensional digits dataset')
3 降维及可视化¶
图片数据是一种高维数据（从几十到上百万的维度），如果把每个图片看成是高维空间的点，要把这些点在高维空间展示出来是极其困难的，所以我们需要将这些数据进行降维，在二维或者三维空间中看出整个数据集的内嵌结构。
本案例要展示的降维方法及所在sklearn的模块如下表所示：
| sklearn 模块 | 调用的方法 | 降维方法 |
|---|---|---|
| random_projecetion | SparseRandomProjection | Random projection |
| decomposition | PCA TruncatedSVD |
PCA truncated SVD |
| discriminant_analysis | LinearDiscriminantAnalysis | LDA |
| manifold | MDS Isomap LocallyLinearEmbedding TSNE SpectralEmbedding |
MDS Isomap LLE t-SNE Spectral Embedding |
| ensemble | RandomTreesEmbedding | Random Trees Embedding |
调用以上方法进行降维的流程都是类似的：
- 首先根据具体方法创建实例：
实例名 = sklearn模块.调用的方法(一些参数的设置) - 然后对数据进行转换：
转换后的数据变量名 = 实例名.fit_transform(X)，在某些方法如LDA降维中还需要提供标签y - 最后将转换后的数据进行可视化：输入转换后的数据以及标题，画出二维空间的图。
为了绘图方便并统一画图的风格，我们首先定义plot_embedding函数用于画出低维嵌入的图形。
配色方案：
- 绿色 #5dbe80
- 蓝色 #2d9ed8
- 紫色 #a290c4
- 橙色 #efab40
- 红色 #eb4e4f
- 灰色 #929591
# 首先定义函数画出二维空间中的样本点，输入参数：1.降维后的数据；2.图片标题
def plot_embedding(X, title=None):
x_min, x_max = np.min(X, 0), np.max(X, 0)
X = (X - x_min) / (x_max - x_min) # 对每一个维度进行0-1归一化，注意此时X只有两个维度
plt.figure(figsize= (6,6)) # 设置整个图形大小
ax = plt.subplot(111)
colors = ['#5dbe80','#2d9ed8','#a290c4','#efab40','#eb4e4f','#929591']
# 画出样本点
for i in range(X.shape[0]): # 每一行代表一个样本
plt.text(X[i, 0], X[i, 1], str(digits.target[i]),
#color=plt.cm.Set1(y[i] / 10.),
color=colors[y[i]],
fontdict={'weight': 'bold', 'size': 9}) # 在样本点所在位置画出样本点的数字标签
# 在样本点上画出缩略图，并保证缩略图够稀疏不至于相互覆盖
# 只有matplotlib 1.0版本以上，offsetbox才有'AnnotationBbox'，所以需要先判断是否有这个功能
if hasattr(offsetbox, 'AnnotationBbox'):
shown_images = np.array([[1., 1.]]) # 假设最开始出现的缩略图在(1,1)位置上
for i in range(digits.data.shape[0]):
dist = np.sum((X[i] - shown_images) ** 2, 1) # 算出样本点与所有展示过的图片（shown_images）的距离
if np.min(dist) < 4e-3: # 若最小的距离小于4e-3，即存在有两个样本点靠的很近的情况，则通过continue跳过展示该数字图片缩略图
continue
shown_images = np.r_[shown_images, [X[i]]] # 展示缩略图的样本点通过纵向拼接加入到shown_images矩阵中
imagebox = offsetbox.AnnotationBbox(
offsetbox.OffsetImage(digits.images[i], cmap=plt.cm.gray_r),
X[i])
ax.add_artist(imagebox)
plt.xticks([]), plt.yticks([]) # 不显示横纵坐标刻度
if title is not None:
plt.title(title)
3.1 Random projection¶
Random projection（随机投影）是最简单的一种降维方法。这种方法只能在很小的程度上展示出整个数据的空间结构，会丢失大部分的结构信息，所以这种降维方法很少会用到。
t0 = time()
rp = random_projection.SparseRandomProjection(n_components=2, random_state=66)
X_projected = rp.fit_transform(X)
plot_embedding(X_projected,
"Random Projection of the digits (time %.2fs)" %
(time() - t0))
3.2 PCA¶
PCA降维是最常用的一种线性的无监督的降维方法。PCA降维实际是对协方差矩阵进行SVD分解来进行降维的线性降维方法。
t0 = time()
pca = decomposition.PCA(n_components=2)
X_pca = pca.fit_transform(X)
plot_embedding(X_pca,
"Principal Components projection of the digits (time %.2fs)" %
(time() - t0))
print pca.explained_variance_ratio_ # 每一个成分对原数据的方差解释了百分之多少
3.3 truncated SVD¶
truncated SVD方法即利用截断SVD分解方法对数据进行线性降维。与PCA不同，该方法在进行SVD分解之前不会对数据进行中心化，这意味着该方法可以有效地处理稀疏矩阵如scipy.sparse定义的稀疏矩阵，而PCA方法不支持scipy.sparse稀疏矩阵的输入。在文本分析领域，该方法可以对稀疏的词频/tf-idf矩阵进行SVD分解，即LSA（隐语义分析）。
t0 = time()
svd = decomposition.TruncatedSVD(n_components=2)
X_svd = svd.fit_transform(X)
plot_embedding(X_svd,
"Principal Components projection of the digits (time %.2fs)" %
(time() - t0))
print svd.explained_variance_ratio_
3.4 LDA¶
LDA降维方法利用了原始数据的标签信息，所以在降维之后，低维空间的相同类的样本点聚在一起，不同类的点分隔比较开。
X2 = X.copy()
X2.flat[::X.shape[1] + 1] += 0.01 # 使得X可逆
t0 = time()
lda = discriminant_analysis.LinearDiscriminantAnalysis(n_components=2)
X_lda = lda.fit_transform(X2, y)
plot_embedding(X_lda,
"Linear Discriminant projection of the digits (time %.2fs)" %
(time() - t0))
print lda.explained_variance_ratio_
3.5 MDS¶
如果样本点之间本身存在距离上的某种类似城市之间地理位置的关系，那么可以使用MDS降维在低维空间保持这种距离关系。
clf = manifold.MDS(n_components=2, n_init=1, max_iter=100)
t0 = time()
X_mds = clf.fit_transform(X)
plot_embedding(X_mds,
"MDS embedding of the digits (time %.2fs)" %
(time() - t0))
3.6 Isomap¶
Isomap是流形学习方法的一种，是等度量映射（Isometric mapping）的简写。该方法可以看做是MDS方法的延伸。不同于MDS的是，该方法在降维前后保持的是所有数据点之间的测地距离（geodesic distances）的关系。
Isomap需要指定领域的最近邻点的个数，我们在提取图片数据时就已经指定为30了。由于需要计算样本点的测地距离，该方法在时间消耗上比较大。
t0 = time()
iso = manifold.Isomap(n_neighbors, n_components=2)
X_iso = iso.fit_transform(X)
plot_embedding(X_iso,
"Isomap projection of the digits (time %.2fs)" %
(time() - t0))
3.7 LLE¶
LLE降维同样需要指定领域样本点个数n_neighbors，LLE降维保持了邻域内的样本点之间的距离关系，它可以理解为一系列的局域PCA操作，但是它在全局上很好的保持了数据的非结构信息。LLE降维主要包括四种方法standard,modified,hessian和ltsa，下面进行一一展示，并且输出它们的重构误差（从低维空间的数据重构原始空间的数据时的误差）。
standard LLE¶
clf = manifold.LocallyLinearEmbedding(n_neighbors, n_components=2, method='standard')
t0 = time()
X_lle = clf.fit_transform(X)
plot_embedding(X_lle,
"Locally Linear Embedding of the digits (time %.2fs)" %
(time() - t0))
print("Reconstruction error: %g" % clf.reconstruction_error_)
modified LLE¶
standard LLE存在正则化的问题：当n_neighbors大于输入数据的维度时，局部邻域矩阵将出现rank-deficient（秩亏）的问题。为了解决这个问题，在standard LLE基础上引入了一个正则化参数$r$。通过设置参数methond='modified'，可以实现modified LLE降维。
clf = manifold.LocallyLinearEmbedding(n_neighbors, n_components=2, method='modified')
t0 = time()
X_mlle = clf.fit_transform(X)
plot_embedding(X_mlle,
"Modified Locally Linear Embedding of the digits (time %.2fs)" %
(time() - t0))
print("Reconstruction error: %g" % clf.reconstruction_error_)
hessian LLE¶
hessian LLE，也叫Hessian Eigenmapping，是另外一种解决LLE正则化问题的方法。该方法使用的前提是满足n_neighbors > n_components * (n_components + 3) / 2。
clf = manifold.LocallyLinearEmbedding(n_neighbors, n_components=2, method='hessian')
t0 = time()
X_hlle = clf.fit_transform(X)
plot_embedding(X_hlle,
"Hessian Locally Linear Embedding of the digits (time %.2fs)" %
(time() - t0))
print("Reconstruction error: %g" % clf.reconstruction_error_)
LTSA¶
LTSA（Local tangent space alignment）方法实际上并不是LLE的变体，只是由于它与LLE算法相似所以被分在了LocallyLinearEmbedding中。与LLE降维前后保持邻近点之间的距离关系不同，LTSA通过将数据放在tangent空间以刻画邻域内样本点之间的地理性质。
clf = manifold.LocallyLinearEmbedding(n_neighbors, n_components=2, method='ltsa')
t0 = time()
X_ltsa = clf.fit_transform(X)
plot_embedding(X_ltsa,
"Local Tangent Space Alignment of the digits (time %.2fs)" %
(time() - t0))
print("Reconstruction error: %g" % clf.reconstruction_error_)
可以看到，standard方法的LLE重构误差最小。
3.8 t-SNE¶
本示例使用PCA生成的嵌入来初始化t-SNE，而不是t-SNE的默认设置（即init='pca'）。 它确保嵌入的全局稳定性，即嵌入不依赖于随机初始化。
t-SNE方法对于数据的局部结构信息很敏感，而且有许多的优点：
- 揭示了属于不同的流形或者簇中的样本
- 减少了样本聚集在
当然，它也有许多缺点：
- 计算代价高，在百万级别的图片数据上需要花费好几小时，而对于同样的任务，PCA只需要花费几分钟或者几秒；
- 该算法具有随机性，不同的随机种子会产生不同的降维结果。当然通过选择不同的随机种子，选取重构误差最小的那个随机种子作为最终的执行降维的参数是可行的；
- 全局结构保持较差，不过这个问题可以通过使用PCA初始样本点来缓解（
init='pca'）。
tsne = manifold.TSNE(n_components=2, init='pca', random_state=10) # 生成tsne实例
t0 = time() # 执行降维之前的时刻
X_tsne = tsne.fit_transform(X) # 降维得到二维空间中的数据
plot_embedding(X_tsne, "t-SNE embedding of the digits (time %.2fs)" % (time() - t0)) # 画出降维后的嵌入图形
plt.show()
3.9 RandomTrees¶
来自sklearn.ensemble模块的RandomTreesEmbedding在技术层面看并不是一种多维嵌入方法，但是它学习了一种数据的高维表示可以用于数据降维方法中。可以先使用RandomTreesEmbedding对数据进行高维表示，然后再使用PCA或者truncated SVD进行降维。
hasher = ensemble.RandomTreesEmbedding(n_estimators=200, random_state=0, max_depth=5)
t0 = time()
X_transformed = hasher.fit_transform(X)
pca = decomposition.TruncatedSVD(n_components=2)
X_reduced = pca.fit_transform(X_transformed)
plot_embedding(X_reduced,
"Random forest embedding of the digits (time %.2fs)" %
(time() - t0))
3.10 Spectral embedding¶
Spectral embedding，也叫Laplacian Eigenmaps（拉普拉斯特征映射）。它利用拉普拉斯图的谱分解找到数据在低维空间中的表示。
embedder = manifold.SpectralEmbedding(n_components=2, random_state=0, eigen_solver="arpack")
t0 = time()
X_se = embedder.fit_transform(X)
plot_embedding(X_se,
"Spectral embedding of the digits (time %.2fs)" %
(time() - t0))
4 总结¶
本案例使用多种降维方法对手写数字图片数据进行降维及可视化展示，包括PCA、LDA和基于流形学习的降维方法等。
线性降维方法包括PCA、LDA时间消耗较少，但是这种线性降维方法会丢失高维空间中的非线性结构信息。相比较而言，非线性降维方法（这里没有提到KPCA和KLDA，有兴趣的可以试一试这两类非线性降维方法）中的流形学习方法可以很好的保留高维空间中的非线性结构信息。虽然典型的流形学习是非监督的方法，但是也存在一些有监督方法的变体。
在进行数据降维时，我们一定要弄清楚我们降维的目的，是为了进行特征提取，使得之后的模型解释性更强或效果提升，还是仅仅为了可视化高维数据。在降维的方法的选择上，我们也要尽量平衡时间成本和降维效果。
另外，在降维时需要注意以下几点：
- 降维之前，所有特征的尺度是一致的；
- 重构误差可以用于寻找最优的输出维度$d$（此时降维不只是为了可视化），随着维度$d$的增加，重构误差将会减小，直到达到实现设定的阈值；
- 噪音点可能会导致流形出现“短路”，即原本在流形之中很容易分开的两部分被噪音点作为一个“桥梁”连接在了一起；
- 某种类型的输入数据可能导致权重矩阵是奇异的，比如在数据集有超过两个样本点是相同的，或者样本点被分在了不相交的组内。在这种情况下，特征值分解的实现
solver='arpack'将找不到零空间。解决此问题的最简单的办法是使用solver='dense'实现特征值分解，虽然dense可能会比较慢，但是它可以用在奇异矩阵上。除此之外，我们也可以通过理解出现奇异的原因来思考解决的办法：如果是由于不相交集合导致，则可以尝试n_neighbors增大；如果是由于数据集中存在相同的样本点，那么可以尝试去掉这些重复的样本点，只保留其中一个。