宣告: 本文轉譯自Data Camp上Manish Pathak的文章《Introduction to t-SNE》原文地址
譯者注: 本文言簡意賅的闡述了資料降維（ Dimensionality Reduction technique）技術中PCA以及t-Distributed Stochastic Neighbor Embedding(t-SNE)演算法的相關實現原理以及利弊，並且使用Python基於Fashion-MNIST資料集描述了對PCA以及t-SNE演算法的基本應用。本人覺得相關概念闡述的比較清晰因此特別轉譯在此部落格，但如果讀完本文後想深究t-SNE背後的數學原理還是推薦看原論文，論文地址會在附錄給出。本人會在原文基礎上新增一些相關注釋（本人學習和相關工作中的一些理解）以[注]為標記。原文在實現的過程中使用的是Python2，為了貼合當下Python生態體系，因此本人使用Python3重新復現。由於本人英文水平以及技術水平有限，我會盡最大的能力確保翻譯得當，倘若文中出現翻譯不恰當的地方還請諸位海涵且幫助我糾正錯誤，可通過文章底部的郵箱地址聯絡我。感謝！
以下為譯文內容 。
通過本教程，你將會對當下較為流行的t-SNE降維技術有一定的認知並能掌握其基礎應用。
本文涉及內容如下：
- 瞭解資料降維以及常用手段型別。
- 瞭解主成分分析（Principal Component Analysis (PCA)）技術原理以及如何在
Python上應用。 - 瞭解t分佈隨機鄰近嵌入（t-Distributed Stochastic Neighbor Embedding (t-SNE)）原理以及如何在
Python上應用。 - 視覺化這兩種演算法的降維結果。
- 比較這兩種演算法之間的優缺點
資料降維
如果你已經有過在大型資料集（包含非常多的特徵）上進行相關工作的經歷，你應該可以深刻（fathom）的體會到想要理解或者挖掘特徵間關係是非常困難的。這一問題不僅體現在資料探索性分析（EDA）上，而且還會影響機器學習模型的效能，因此你可能會因此讓你的模型過擬合或者破壞了你對模型的某些假設，例如獨立特徵對線性迴歸模型的影響[1]。這時資料降維就顯得格外重要。在機器學習中，降維技術是減少隨機變數個數，得到一組“不相關”主變數的過程。通過減少特徵空間的維數，你僅需要考慮一小部分特徵之間的關係，而且你可以很輕易地在這部分特徵上做探索或者視覺化，同時也可以減少你的模型出現過擬合的可能。
[注1] 考慮線性相關的兩個變數以及線性無關的兩個變數對整個模型的影響。
降維通常可以通過以下方式實現：
- 特徵限制：通過限制特徵來縮減特徵空間。雖然可以達到降維目的，但其缺點就是你會丟失一部分被你刪掉的特徵所蘊含的資訊。
- 特徵選擇：應用一些統計檢驗的方法去根據它們的重要程度進行排序，然後選擇其中最主要的幾個特徵子集。這類方法依然會面臨資訊丟失的風險，而且還會因為不同的檢驗方法會導致出現不同的特徵排序方式，因此是不穩定的。
- 特徵提取：通過建立一些融合了若干舊特徵的獨立特徵。這類技術可以很好的應用線上性或非線性的降維技術中。
主成分分析（PCA）
主成分分析是一種線性的特徵提取技術，它將資料通過線性對映的方式投影到低維空間中，通過這樣的方式能夠確保原資料在低維空間中方差最大。它通過計算其特徵的協方差矩陣中的特徵向量來實現這一目的。與最大的特徵值（主成分）一一對應的特徵向量會被用於重建成新的資料，並且保證這些資料在該特徵向量方向上的方差最大。
簡單來說，PCA以特定的方式融合了你所有的輸入屬性（特徵）值，這樣你就可以在刪除不重要特徵的同時不必擔心丟失最有價值的部分。還有一個更為顯著的好處，經過PCA處理之後的每一個新特徵都是獨立於其他特徵的[2]。
[注2] 矩陣分解中的所有特徵向量都是線性無關的。
t分佈隨機鄰近插入（t-SNE）
t-SNE是一種非線性的降維技術，非常適合用於高維資料的視覺化。廣泛應用於影像處理、自然語言處理，基因資料以及語音處理。為了保證足夠淺顯易懂，這裡僅對t-SNE的工作原理做簡要介紹：
- 該演算法一開始通過計算在高維空間中的資料點的相似度概率和與其對應的低維空間中的點的相似度的概率。點的相似度計算方法是[3]：以A為中心的高斯分佈中，如果按概率密度的比例選取相鄰點，則點A將選擇點B作為其相鄰點的條件概率，以此計算點A的相似性。
- 為了更好將資料投影至低維空間中，演算法嘗試去最小化高維資料空間和低維資料空間之間的條件概率（相似度）之差。
- 為了評估
t-SNE條件概率差和的最小化，使用梯度下降的方法最小化原分佈中資料與對映分佈中的對應資料的KL散度[4]（Kullback-Leibler divergence）的總和。
[注3] 更容易理解的方式是通過歐式距離算出的AB兩點的距離轉換為條件概率以此來表達點與點之間的相似度，此概率與大，AB兩點的相似度就越高。
[注4] 即相對熵或稱資訊增益。讓其值變小也就是為了讓相似度更高的資料點聚集在一起。資訊增益小則說明區分該例項的難度大，換個角度來說就是這兩個例項非常相似。關於資訊增益相關概念可以瀏覽一下我的另一篇博文
如果有興趣更進一步研究t-SNE的同學可以檢視附錄中的論文1。
簡單來說，t-SNE主要目的在最小化兩種分佈的差異性：第一個是度量輸入例項成對相似性，另一種是嵌入在相應低維空間中的成對點的相似性。
通過上述方式，t-SNE對映多維資料到低維空間中並且試圖找到基於多個特徵的資料點的相似性來區分觀測資料群，從而發現資料中的模式。然而，經過這一過程後，輸入特徵開始變得模糊，並且你無法僅基於t-SNE的結果進行任何推斷[5]。因此這也是為什麼t-SNE主要還是用來做EDA和視覺化的原因。
[注5] 與PCA比較就可以很顯然的看出，經過PCA處理過後的結果能夠得知每一個成分的方差貢獻度（解釋方差），然後t-SNE僅僅是基於相似度進行判定，沒辦法從其結果推斷類似的資訊。
t-SNE應用 Python實現
現在你可以在Python中基於開源資料集應用t-SNE演算法並且將其降維結果視覺化。與此同時，你也要在相同資料集應用PCA演算法並視覺化結果，然後與t-SNE進行比較。
本次的資料集使用Fashion-MNIST並且您可以點選此處進行下載。
Fashion-MNIST是類似MNIST手寫影像資料集的公共資料集，由70,000條已標註為10種類別的時尚服裝資料，每一個例項都是由28x28的灰度影像組成，其中訓練資料集含有60,000條，測試資料集有10,000條，用此資料集比用原MNIST資料集能更好的對比結果。Fashion-MNIST的標籤與MNIST一樣是0-9，但是不同的是這是個數字代表的是對應的時尚服裝產品，下面對每一個數字對應的含義進行解釋說明：
| 標註編號 | 描述 |
|---|---|
| 0 | T-shirt/top（T恤） |
| 1 | Trouser（褲子） |
| 2 | Pullover（套衫） |
| 3 | Dress（裙子） |
| 4 | Coat（外套） |
| 5 | Sandal（涼鞋） |
| 6 | Shirt（汗衫） |
| 7 | Sneaker（運動鞋） |
| 8 | Bag（包） |
| 9 | Ankle boot（踝靴） |
之後你可以在Fashion-MNIST官方倉庫的utils檔案下的mnist_reader.py找到讀取該資料集的特定方法，即load_mnist()，如下：
def load_mnist(path, kind='train'):
import os
import gzip
import numpy as np
"""Load MNIST data from `path`"""
labels_path = os.path.join(path,
'%s-labels-idx1-ubyte.gz'
% kind)
images_path = os.path.join(path,
'%s-images-idx3-ubyte.gz'
% kind)
with gzip.open(labels_path, 'rb') as lbpath:
labels = np.frombuffer(lbpath.read(), dtype=np.uint8,
offset=8)
with gzip.open(images_path, 'rb') as imgpath:
images = np.frombuffer(imgpath.read(), dtype=np.uint8,
offset=16).reshape(len(labels), 784)
return images, labels
複製程式碼之後你可以通過load_mnist()方法讀取訓練和測試資料集並將其應用在你的演算法上即可，只要將資料集的.gz檔案路徑作為第一個引數以及資料型別kind作為第二個引數傳入該函式即可，如
X_train, y_train = load_mnist('.', kind='train') # 我的.gz檔案就放在當前路徑下
複製程式碼讀取資料之後，你可以檢查一下訓練資料的基本屬性，如shape，你會看到訓練資料由60,000個例項以及784個特徵組成
X_train.shape
複製程式碼(60000, 784)
複製程式碼同樣的標註資料也可以看到是由0-910個標籤組成
y_train
複製程式碼array([9, 0, 0, ..., 3, 0, 5], dtype=uint8)
複製程式碼接下來，為了保證能夠程式碼的完整性比如引入所需要的第三方庫，同時為了確保可復現，還需要設定Random State引數[6]為123。程式碼如下：
[注6] 我們知道計算機的隨機都是偽隨機，因此為了確保程式碼結果是可復現的都會設定一個隨機因子，至於這個值是多少並沒有規定，例如我本人就喜歡設定成42，原因是在《銀河系漫遊指南》中42是超級計算機得出的生命終極答案。
import time
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patheffects as PathEffects
%matplotlib inline
import seaborn as sns
sns.set_style('darkgrid')
sns.set_palette('muted')
sns.set_context("notebook", font_scale=1.5,
rc={"lines.linewidth": 2.5})
RS = 123
複製程式碼為了能夠相容兩個演算法結果的視覺化展示，需要建立一個fashion_scatter()函式。該函式接受兩個引數：引數1 x 接受一個演算法結果的2維矩陣輸入；引數2 colors接受一個1維的標籤陣列。該函式會根據colors中的標籤對x散點資料繼續染色。
fashion_scatter()定義如下：
# Utility function to visualize the outputs of PCA and t-SNE
def fashion_scatter(x, colors):
# choose a color palette with seaborn.
num_classes = len(np.unique(colors))
palette = np.array(sns.color_palette("hls", num_classes))
# create a scatter plot.
f = plt.figure(figsize=(8, 8))
ax = plt.subplot(aspect='equal')
sc = ax.scatter(x[:,0], x[:,1], lw=0, s=40, c=palette[colors.astype(np.int)])
plt.xlim(-25, 25)
plt.ylim(-25, 25)
ax.axis('off')
ax.axis('tight')
# add the labels for each digit corresponding to the label
txts = []
for i in range(num_classes):
# Position of each label at median of data points.
xtext, ytext = np.median(x[colors == i, :], axis=0)
txt = ax.text(xtext, ytext, str(i), fontsize=24)
txt.set_path_effects([
PathEffects.Stroke(linewidth=5, foreground="w"),
PathEffects.Normal()])
txts.append(txt)
return f, ax, sc, txts
複製程式碼為了不讓你的機器承擔過重的記憶體與執行時間壓力，在本次實踐中我們只取訓練集中前20,000條作為訓練資料，同時你要確保這20,000條資料一定要涵蓋10個標籤。
x_subset = X_train[0:20000]
y_subset = y_train[0:20000]
np.unique(y_subset)
複製程式碼array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=uint8)
複製程式碼現在可以在訓練子集上使用PCA演算法並視覺化其結果了，為了方便我們直接使用scikit-learn提供的PCA演算法，將n_components設定為4，這個引數決定了最後輸出的資料維度。關於PCA的更多用法可以自行前往scikit-learn的官網檢視
from sklearn.decomposition import PCA
time_start = time.time()
pca = PCA(n_components=4)
pca_result = pca.fit_transform(x_subset)
print(f'PCA done! Time elapsed: {time.time()-time_start} seconds')
複製程式碼PCA done! Time elapsed: 0.6091551780700684 seconds
複製程式碼現在可以將pca_reuslt存入DataFrame中，之後我們可以檢查這四個主成分的各自的資料方差。
pca_df = pd.DataFrame(columns = ['pca1','pca2','pca3','pca4'])
pca_df['pca1'] = pca_result[:,0]
pca_df['pca2'] = pca_result[:,1]
pca_df['pca3'] = pca_result[:,2]
pca_df['pca4'] = pca_result[:,3]
print(f'Variance explained per principal component: {pca.explained_variance_ratio_}')
複製程式碼Variance explained per principal component: [0.29021329 0.1778743 0.06015076 0.04975864]
複製程式碼注意：訓練子集中的第一和第二主成分的解釋方差之和幾乎達到了48%。因此我們僅對這兩個主成分的資料進行視覺化即可。
top_two_comp = pca_df[['pca1','pca2']] # taking first and second principal component
f, ax, sc, txts = fashion_scatter(top_two_comp.values,y_subset) # Visualizing the PCA output
f.show()
複製程式碼如圖中所示，PCA演算法試圖將不同點區分且將相同的點聚集起來。因此，這個圖也可以用來做資料探索性分析。當然主成分1與主成分2還可以在分類和聚類演算法上直接充當特徵。
現在我們使用t-SNE來做相同的事情，當然t-SNE也在被scikit-learn實現了，我們直接使用即可，在此我們提供一些t-SNE常用的引數說明，更多關於t-SNE的說明可以參考scikit-learn的官方文件
- n_components (預設為2): 嵌入空間的維數。
- perplexity (預設為30): 複雜度的含義與其他流行學習方法[7]中的最鄰近個數的含義相似，通常在5-50之間考慮。
- early_exaggeration (預設為12.0): 控制原始空間在嵌入空間中的密度集大小。
- learning_rate (預設為200.0): 學習率的常見範圍在10.0~1000.0之間。
- n_iter (預設為1000): 演算法優化的的最大迭代次數，至少應該設定為250次。
- method (預設為 ‘barnes_hut’): 學習方法，Barnes-Hut方法執行時間為O(NlogN)。method='exact'會比較慢，演算法執行時間為O(N2)。
[注7] 流行學習方法是指從高維取樣資料中恢復低維流形結構，即找到高維空間中的低維流形。引自簡書《流行學習-實現高維資料的降維與視覺化》。
在此我們先以預設引數執行t-SNE演算法：
from sklearn.manifold import TSNE
import time
time_start = time.time()
fashion_tsne = TSNE(random_state=RS, n_jobs=-1).fit_transform(x_subset)
print(f't-SNE done! Time elapsed: {time.time()-time_start} seconds')
複製程式碼t-SNE done! Time elapsed: 882.41050598 seconds
複製程式碼很顯然，t-SNE在相同資料集上的執行時間要比PCA長太多[8]。
[注8] 在sklearn 0.22版本更新了n_jobs引數，可以使用多個CPU平行計算，會對執行時間有所改觀。程式碼如下：
fashion_tsne = TSNE(random_state=RS, n_jobs=-1).fit_transform(x_subset)
複製程式碼因機器而異，在我本地跑完整個演算法耗時223秒。
現在視覺化t-SNE的結果：
f, ax, sc, txts = fashion_scatter(fashion_tsne, y_subset)
f.show()
複製程式碼正如你所見，這次的降維結果比剛才PCA的結果好了許多，基本可以清晰的看到每個標籤都被聚集在它們各自的簇中，如果你將此資料應用於聚類演算法中，應該會更精確的描繪每一個簇。
在sklearn關於t-SNE的文件中有提到：
It is highly recommended to use another dimensionality reduction method (e.g., PCA for dense data or TruncatedSVD for sparse data) to reduce the number of dimensions to a reasonable amount (e.g., 50) if the number of features is very high. This will suppress some noise and speed up the computation of pairwise distances between samples.
當如果維數非常高時，強烈建議使用其他的降維方式（如稠密矩陣使用PCA或者稀疏矩陣使用TruncatedSVD）去降低至一個合理的維度（如50），這樣做會極大的抑制資料噪聲，加快樣本間歐氏距離的計算時間。
現在我們可以根據官方給出的這個建議，先使用PCA演算法[9]將資料維數降至50維，然後再在t-SNE上使用此資料，並視覺化其結果。
[注9] 個人感覺這裡分別嘗試PCA與TruncatedSVD這兩中方法會更加合理一些。
time_start = time.time()
pca_50 = PCA(n_components=50)
pca_result_50 = pca_50.fit_transform(x_subset)
print(f'PCA with 50 components done! Time elapsed: {time.time()-time_start} seconds')
print(f'Cumulative variance explained by 50 principal components: {np.sum(pca_50.explained_variance_ratio_)}')
複製程式碼PCA with 50 components done! Time elapsed: 1.1455249786376953 seconds
Cumulative variance explained by 50 principal components: 0.8624682420611026
複製程式碼現在將pca_result_50應用在t-SNE上：
import time
time_start = time.time()
# 這裡直接應用 n_jobs=-1 引數，啟用所有cpu進行計算
fashion_pca_tsne = TSNE(random_state=RS, n_jobs=-1).fit_transform(pca_result_50)
print(f't-SNE done! Time elapsed: {time.time()-time_start} seconds')
複製程式碼t-SNE done! Time elapsed: 93.07448196411133 seconds
複製程式碼可以看出經過以上步驟後的執行時間大幅度縮減。
從圖上可以看出來，幾乎非常接近於直接運用t-SNE的圖，但是還是有一些地方不同，在0標籤或者說在T-shirt/top這一簇被更加緊密的聚集在一起了。
PCA 與 t-SNE 的區別
儘管PCA與t-SNE各有優劣，但是它們之間還是有一些非常重要的區別：
t-SNE的計算成本非常高，相同數量級的資料集在PCA上的計算速度會比t-SNE要快許多。- 從理論上來說，
PCA是一種矩陣分解技術，而t-SNE是一種概率方法。 - 在類似
PCA一樣的線性降維演算法中，會將不同的資料點置於距離較遠的低維空間中。但是，為了在低維非線性流行上表示高維資料，必須將相似的資料點緊密的表示在一起，這也是t-SNE與PCA應用場景不同之處。 - 在使用
t-SNE的時候，即使是相同的超引數但是由於在不同時期執行的結果可能不盡相同，因此在使用t-SNE時必須觀察許多圖，而PCA則是穩定的。 - 由於
PCA是一種線性的演算法，它無法解釋特徵之間的複雜多項式關係也即非線性關係，而t-SNE可以獲知這些資訊。
結論
恭喜！你通過本教程已經對降維技術有了一定的瞭解，並且知道如何使用PCA和t-SNE這兩種主流的降維技術，並且瞭解瞭如何建立一些漂亮的圖來比較它們之間的結果。但是，這些僅僅是非常基礎的知識，在t-SNE的原理中還有許多知識可以去深究，希望通過本教程可以讓您在您的日常工作中更好的使用t-SNE。