目錄
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.7.1、standard LLE
3.7.2、modified LLE
3.7.3、hessian LLE
3.7.4、LTSA
3.8、t-SNE降維
3.9、RandomTrees降維
3.10、Spectral embedding降維
4、總結
降維是通過單幅圖像資料的高維化，對單幅圖像轉化為高維空間中的資料集合進行的一種操作。機器學習領域中所謂的降維就是指采用某種映射方法，将原高維空間中的資料點映射到低次元的空間中。降維的本質是學習一個映射函數 f : x->y，其中x是原始資料點的表達，目前最多使用向量表達形式。 y是資料點映射後的低維向量表達，通常y的次元小于x的次元（當然提高次元也是可以的）。f可能是顯式的或隐式的、線性的或非線性的。
本項目将依托于MNIST資料集，手把手實作圖像資料集降維。
MNIST資料集來自美國國家标準與技術研究所，是入門級的計算機視覺資料集。它是由6萬張訓練圖檔和1萬張測試圖檔構成的，這些圖檔是手寫的從0到9的數字，50%采集美國中學生，50%來自人口普查局(the Census Bureau)的從業人員。這些數字圖檔進行過預處理和格式化，均為黑白色構成，做了大小調整（28×28像素）并居中處理。MNIST資料集效果如下圖所示：
在本案例中，選擇直接從sklearn.datasets子產品中通過load_digits導入手寫數字圖檔資料集，該資料集是UCI datasets的Optical Recognition of Handwritten Digits Data Set中的測試集，并且隻是MNIST的很小的子集，一共有1797張分辨率為8××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) 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) 輸出結果如下所示：
[[ 0. 0. 5. ..., 0. 0. 0.]
[ 0. 0. 0. ..., 10. 0. 0.]
[ 0. 0. 0. ..., 16. 9. 0.]
...,
[ 0. 0. 0. ..., 9. 0. 0.]
[ 0. 0. 0. ..., 4. 0. 0.]
[ 0. 0. 6. ..., 6. 0. 0.]]
[0 1 2 ..., 4 4 0] 方法二：隻傳回data和target，實作代碼如下所示：
data = datasets.load_digits(n_class=6)
print(data) {'images': array([[[ 0., 0., 5., ..., 1., 0., 0.],
[ 0., 0., 13., ..., 15., 5., 0.],
[ 0., 3., 15., ..., 11., 8., 0.],
...,
[ 0., 4., 11., ..., 12., 7., 0.],
[ 0., 2., 14., ..., 12., 0., 0.],
[ 0., 0., 6., ..., 0., 0., 0.]],
[[ 0., 0., 0., ..., 5., 0., 0.],
[ 0., 0., 0., ..., 9., 0., 0.],
[ 0., 0., 3., ..., 6., 0., 0.],
...,
[ 0., 0., 1., ..., 6., 0., 0.],
[ 0., 0., 1., ..., 6., 0., 0.],
[ 0., 0., 0., ..., 10., 0., 0.]],
[[ 0., 0., 0., ..., 12., 0., 0.],
[ 0., 0., 3., ..., 14., 0., 0.],
[ 0., 0., 8., ..., 16., 0., 0.],
...,
[ 0., 9., 16., ..., 0., 0., 0.],
[ 0., 3., 13., ..., 11., 5., 0.],
[ 0., 0., 0., ..., 16., 9., 0.]],
...,
[[ 0., 0., 0., ..., 6., 0., 0.],
[ 0., 0., 0., ..., 2., 0., 0.],
[ 0., 0., 8., ..., 1., 2., 0.],
...,
[ 0., 12., 16., ..., 16., 1., 0.],
[ 0., 1., 7., ..., 13., 0., 0.],
[ 0., 0., 0., ..., 9., 0., 0.]],
[[ 0., 0., 0., ..., 4., 0., 0.],
[ 0., 0., 4., ..., 0., 0., 0.],
[ 0., 0., 12., ..., 4., 3., 0.],
...,
[ 0., 12., 16., ..., 13., 0., 0.],
[ 0., 0., 4., ..., 8., 0., 0.],
[ 0., 0., 0., ..., 4., 0., 0.]],
[[ 0., 0., 6., ..., 11., 1., 0.],
[ 0., 0., 16., ..., 16., 1., 0.],
[ 0., 3., 16., ..., 13., 6., 0.],
...,
[ 0., 5., 16., ..., 16., 5., 0.],
[ 0., 1., 15., ..., 16., 1., 0.],
[ 0., 0., 6., ..., 6., 0., 0.]]]), 'data': array([[ 0., 0., 5., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 10., 0., 0.],
[ 0., 0., 0., ..., 16., 9., 0.],
...,
[ 0., 0., 0., ..., 9., 0., 0.],
[ 0., 0., 0., ..., 4., 0., 0.],
[ 0., 0., 6., ..., 6., 0., 0.]]), 'target_names': array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]), 'DESCR': "Optical Recognition of Handwritten Digits Data Set\n===================================================\n\nNotes\n-----\nData Set Characteristics:\n :Number of Instances: 5620\n :Number of Attributes: 64\n :Attribute Information: 8x8 image of integer pixels in the range 0..16.\n :Missing Attribute Values: None\n :Creator: E. Alpaydin (alpaydin '@' boun.edu.tr)\n :Date: July; 1998\n\nThis is a copy of the test set of the UCI ML hand-written digits datasets\nhttp://archive.ics.uci.edu/ml/datasets/Optical+Recognition+of+Handwritten+Digits\n\nThe data set contains images of hand-written digits: 10 classes where\neach class refers to a digit.\n\nPreprocessing programs made available by NIST were used to extract\nnormalized bitmaps of handwritten digits from a preprinted form. From a\ntotal of 43 people, 30 contributed to the training set and different 13\nto the test set. 32x32 bitmaps are divided into nonoverlapping blocks of\n4x4 and the number of on pixels are counted in each block. This generates\nan input matrix of 8x8 where each element is an integer in the range\n0..16. This reduces dimensionality and gives invariance to small\ndistortions.\n\nFor info on NIST preprocessing routines, see M. D. Garris, J. L. Blue, G.\nT. Candela, D. L. Dimmick, J. Geist, P. J. Grother, S. A. Janet, and C.\nL. Wilson, NIST Form-Based Handprint Recognition System, NISTIR 5469,\n1994.\n\nReferences\n----------\n - C. Kaynak (1995) Methods of Combining Multiple Classifiers and Their\n Applications to Handwritten Digit Recognition, MSc Thesis, Institute of\n Graduate Studies in Science and Engineering, Bogazici University.\n - E. Alpaydin, C. Kaynak (1998) Cascading Classifiers, Kybernetika.\n - Ken Tang and Ponnuthurai N. Suganthan and Xi Yao and A. Kai Qin.\n Linear dimensionalityreduction using relevance weighted LDA. School of\n Electrical and Electronic Engineering Nanyang Technological University.\n 2005.\n - Claudio Gentile. A New Approximate Maximal Margin Classification\n Algorithm. NIPS. 2000.\n", 'target': array([0, 1, 2, ..., 4, 4, 0])} 本案例隻提取從數字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 将資料集部分資料圖檔可視化顯示，實作代碼如下所示：
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') 圖檔資料是一種高維資料（從幾十到上百萬的次元），如果把每個圖檔看成是高維空間的點，要把這些點在高維空間展示出來是極其困難的，是以我們需要将這些資料進行降維，在二維或者三維空間中看出整個資料集的内嵌結構。
本案例要展示的降維方法及所在sklearn的子產品如下表所示：
調用以上方法進行降維的流程都是類似的：
首先根據具體方法建立執行個體：執行個體名 = 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) 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)) 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_ # 每一個成分對原資料的方差解釋了百分之多少 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)) 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)) 如果樣本點之間本身存在距離上的某種類似城市之間地理位置的關系，那麼可以使用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)) 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)) LLE降維同樣需要指定領域樣本點個數n_neighbors，LLE降維保持了鄰域内的樣本點之間的距離關系，它可以了解為一系列的局域PCA操作，但是它在全局上很好的保持了資料的非結構資訊。LLE降維主要包括四種方法standard,modified,hessian和ltsa，下面進行一一展示，并且輸出它們的重構誤差（從低維空間的資料重構原始空間的資料時的誤差）。
3.7.1、standard LLE
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)) 效果如下圖所示：
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)) 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)) 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)) 本示例使用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() 來自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)) 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)) 本案例使用多種降維方法對手寫數字圖檔資料進行降維及可視化展示，包括PCA、LDA和基于流形學習的降維方法等。
線性降維方法包括PCA、LDA時間消耗較少，但是這種線性降維方法會丢失高維空間中的非線性結構資訊。相比較而言，非線性降維方法（這裡沒有提到KPCA和KLDA，有興趣的可以試一試這兩類非線性降維方法）中的流形學習方法可以很好的保留高維空間中的非線性結構資訊。雖然典型的流形學習是非監督的方法，但是也存在一些有監督方法的變體。
在進行資料降維時，我們一定要弄清楚我們降維的目的，是為了進行特征提取，使得之後的模型解釋性更強或效果提升，還是僅僅為了可視化高維資料。在降維的方法的選擇上，我們也要盡量平衡時間成本和降維效果。
另外，在降維時需要注意以下幾點：
降維之前，所有特征的尺度是一緻的；
重構誤差可以用于尋找最優的輸出次元𝑑d（此時降維不隻是為了可視化），随着次元𝑑d的增加，重構誤差将會減小，直到達到實作設定的門檻值；
噪音點可能會導緻流形出現“短路”，即原本在流形之中很容易分開的兩部分被噪音點作為一個“橋梁”連接配接在了一起；
某種類型的輸入資料可能導緻權重矩陣是奇異的，比如在資料集有超過兩個樣本點是相同的，或者樣本點被分在了不相交的組内。在這種情況下，特征值分解的實作solver='arpack'将找不到零空間。解決此問題的最簡單的辦法是使用solver='dense'實作特征值分解，雖然dense可能會比較慢，但是它可以用在奇異矩陣上。除此之外，我們也可以通過了解出現奇異的原因來思考解決的辦法：如果是由于不相交集合導緻，則可以嘗試n_neighbors增大；如果是由于資料集中存在相同的樣本點，那麼可以嘗試去掉這些重複的樣本點，隻保留其中一個。