|
UMAP應該說是目前最好的降維算法了，能最大程度的保留原始數據的特征同時大幅度的降低特征維數。 
這是《生命科學的數理統計和機器學習》的相關探討，我試圖介紹生物信息學、生物醫學、遺傳學等常見的分析技術。今天，我們將深入探討一種令人興奮的降維技術，稱為UMAP，它主導了當今的單細胞基因組學。在這里，我將嘗試質疑UMAP作為一種過于數學化的方法的神話，并使用簡單的語言來解釋它。在本文章中，我將展示如何在Python中從頭開始編寫UMAP，以及(更高級的)如何創建一種比UMAP提供更好的可視化效果的降維技術。然而，現在我們將從UMAP背后的故事開始，并強調tSNE和UMAP之間的關鍵區別。
如果你不知道tSNE是什么，它是如何工作的，也沒有讀過2008年的革命性的van der Maaten & Hinton原稿，你可能不需要知道，因為tSNE現在基本上已經死了。盡管tSNE對一般的單細胞基因組學和數據科學產生了巨大的影響，但人們普遍認為它有一些缺點，這些缺點很快將得到解決。 
究竟是什么讓我們想放棄使用tSNE單細胞基因組學?這里我用簡短的評論總結了幾點: 在scRNAseq中，tSNE并不能很好地擴展快速增加的樣本量。試圖用FItSNE來加速它會導致大量的內存消耗，使得在計算機集群之外無法進行分析。（集群之痛） tSNE沒有保留全局數據結構，這意味著只有在集群距離內才有意義，而集群之間的相似性并不能保證，因此使用tSNE作為集群不是一個很好的主意。 
tSNE實際上只能嵌入到2維或3維中，即僅用于可視化目的，因此很難將tSNE作為一種通用的降維技術來生產10個或50個組件。 請注意，這對于更現代的FItSNE算法來說仍然是一個問題。 tSNE執行從高維到低維的非參數映射，這意味著它不利用驅動觀察到的集群的特性(也稱為PCA加載)。 tSNE不能直接處理高維數據，在將其插入tSNE之前，通常使用自動編碼器或PCA執行預降維 - tSNE的計算占用了大量的內存，在使用大型perplexity超參數時，這一點變得尤為明顯，因為k近鄰的初始步長(就像Barnes-Hut過程)變得低效，而且對時間減少也很重要。更現代的FItSNE算法也不能解決這個問題。
tSNE是一種比較簡單的機器學習算法，可以用以下四個公式表示:
(1)定義了高維空間中任意兩點之間觀測距離的高斯概率，滿足對稱性規則。Eq。 (2)引入了困惑的概念作為一個約束,確定最優σ為每個樣本。 (3)聲明了低維嵌入中點對之間距離的學生t分布。學生t分布是為了克服嵌入低維時的擁擠問題。 (4)給出了Kullback-Leibler散度損失函數，將高維概率投影到低維概率上，并給出了梯度下降優化中使用的梯度的解析形式。 
看看上面的圖，我想說的是 t分布應該提供全局距離信息，因為它們將高維空間中相距較遠的點推到低維空間中更遠的點。 然而，這種良好的意愿被成本函數(KL-divergence)的選擇所扼殺，我們將在后面看到其原因。 當了解UMAP時，我的第一印象是這是一種全新的、有趣的降維技術，它基于堅實的數學原理，因此與純機器學習半經驗算法tSNE非常不同。我的生物學同事告訴我，最初的UMAP論文“太數學化了”，看著論文的第二部分，我很高興地看到，嚴格而準確的數學終于進入生命和數據科學。然而，在閱讀UMAP文檔和觀看Leland McInnes在2018年SciPy大會上的演講時，我感到困惑，覺得UMAP是另一種鄰居圖技術，它與tSNE如此相似，以至于我很難理解UMAP與tSNE究竟有什么不同。 從UMAP論文中，雖然Leland McInnes試圖在附錄c中總結UMAP和tSNE之間的差異，但是它們之間的差異并不是很明顯。在這里，我將首先總結我所注意到的UMAP和tSNE之間的不同之處，然后嘗試解釋為什么這些不同是重要的，并找出它們的影響有多大。 UMAP在高維中使用指數概率分布，但不一定是像tSNE那樣的歐氏距離，而是任何距離都可以代入。另外，概率沒有歸一化:
ρ是一個重要的參數,它代表了每個i數據點距離第一個最近鄰。這保證了流形的局部連接性。換句話說，這為每個數據點提供了一個局部自適應指數核，因此距離度量在點與點之間變化。 ρ參數是唯一在UMAP部分2和3之間的橋梁。否則,我也看不出什么模糊單形建筑,即從第二節的拓撲數據分析,與從第三節UMAP的算法實現,似乎在一天結束的時候模糊單形集導致最近鄰圖施工。 UMAP對高維或低維概率都不應用標準化，這與tSNE非常不同，感覺很奇怪。然而,高或低維函數形式的概率可以看到他們已經擴展[0,1],事實證明,缺乏規范化的情況下,類似于情商分母。
(1),可以顯著降低計算時間高維圖像由于求和或集成是一個代價高昂的計算過程。想想馬爾可夫鏈蒙特卡羅(MCMC)它基本上是試圖近似地計算在貝葉斯規則的分母上的積分（UMAP使用最近鄰居的數量而不是perplexity） (2)定義perplexity, UMAP則定義了沒有log2函數的最近鄰居k的個數，即: 
symmterization是必要的因為UMAP融合在一起的點與本地不同的指標(通過參數ρ),它可能發生圖A和B節點之間的重量不等于B之間的權重和節點。為什么UMAP使用這種對稱而不是tSNE使用的對稱還不清楚。我將在下一篇文章(從頭開始編寫UMAP)中展示我對不同的對稱化規則的實驗，這并沒有使我相信這是如此重要的一步，因為它對最終的低維嵌入式產生了很小的影響。UMAP使用曲線族1 / (1+a*y^(2b))在低維中建模距離概率，不是完全的學生t分布，但非常非常相似，請注意再次沒有應用標準化:
其中，對于默認UMAP超參數a≈1.93,b≈0.79(實際上，對于min_dist = 0.001)。在實踐中，UMAP從非線性最小二乘擬合到帶有min_dist超參數的分段函數中找到a和b: 
為了更好地理解曲線族1 / (1+a*y^(2b))的行為，讓我們畫出不同a和b的曲線:plt.figure(figsize=(20, 15))y = np.linspace(0, 10, 1000)
my_prob = lambda y, a, b: np.power(1 + a*y**(2*b), -1)plt.plot(y, my_prob(y, a = 1, b = 1))plt.plot(y, my_prob(y, a = 1.93, b = 0.79))plt.plot(y, my_prob(y, a = 1.93, b = 5))
plt.gca().legend(('a = 1, b = 1', 'a = 1.93, b = 0.79', 'a = 1.93, b = 5'), fontsize = 20)plt.title('Low-Dimensional Probability of Pairwise Euclidean Distances', fontsize = 20)plt.xlabel('Y', fontsize = 20); plt.ylabel('Q(Y)', fontsize = 20)plt.show()
我們可以看到曲線族對參數b非常敏感，在大的參數b處，在小的參數y處，曲線族形成了一種高峰。這意味著在UMAP超參數min_dist之下，所有的數據點都是同樣緊密相連的。由于Q(Y)函數的行為幾乎像一個Heaviside階躍函數，這意味著UMAP為所有在低維空間中相互靠近的點分配了幾乎相同的低維坐標。min_dist正是導致在UMAP維數降低圖中經常觀察到的超密集集群的原因。 為了演示如何準確地找到a和b參數，讓我們展示一個簡單的分段函數(其中高峰部分是通過min_dist參數定義的)，并使用函數族1 / (1+a*y^(2b))通過優化來擬合它。curve_fit來自Scipy Python庫。作為擬合的結果，我們得到了函數1 / (1+a*y^(2b))的初值a和初值b。 from scipy import optimizeimport matplotlib.pyplot as pltimport numpy as np
MIN_DIST = 1
x = np.linspace(0, 10, 300)
def f(x, min_dist): y = [] for i in range(len(x)): if(x[i] <= min_dist): y.append(1) else: y.append(np.exp(- x[i] + min_dist)) return y
dist_low_dim = lambda x, a, b: 1 / (1 + a*x**(2*b)) p , _ = optimize.curve_fit(dist_low_dim, x, f(x, MIN_DIST))print(p)
plt.figure(figsize=(20,15))plt.plot(x, f(x, MIN_DIST), 'o')plt.plot(x, dist_low_dim(x, p[0], p[1]), c = 'red')plt.title('Non-Linear Least Square Fit of Piecewise Function', fontsize = 20)plt.gca().legend(('Original', 'Fit'), fontsize = 20)plt.xlabel('X', fontsize = 20)plt.ylabel('Y', fontsize = 20)plt.show()
UMAP使用二元交叉熵(CE)作為成本函數，而不是像tSNE那樣使用kl -散度。
在下一節中，我們將展示CE成本函數中的這個附加(第二個)項使UMAP能夠捕獲全局數據結構，而tSNE只能在適當的perplexity值上對局部結構建模。由于我們需要知道交叉熵的梯度，以便以后實現梯度下降，讓我們快速計算它。忽略只包含p(X)的常數項，我們可以將交叉熵重新寫一下，并將其微分如下:
UMAP使用圖拉普拉斯變換分配初始的低維坐標，與tSNE使用的隨機正常初始化形成對比。 然而，這應該會對最終的低維表示產生較小的影響，至少在tSNE中是這樣的。 但是，這應該會減少UMAP每次運行時的更改，因為它不再是一個隨機初始化。 Linderman和Steinerberger提出了一個有趣的假設，即在tSNE的初始階段最小化kl散度，并在早期進行放大，這一假設的動機是通過圖拉普拉斯變換來進行初始化。
圖拉普拉斯、譜聚類、拉普拉斯Eignemaps、擴散圖、譜嵌入等，實際上是指將矩陣分解和鄰接圖方法結合起來解決降維問題的同一種有趣的方法。在這種方法中，我們首先構造一個圖(或knn圖)，然后通過構造拉普拉斯矩陣用矩陣代數(鄰接矩陣和度矩陣)將其形式化，最后分解拉普拉斯矩陣，即求解特征值分解問題。
我們可以使用scikit-learn Python庫，并使用spectralembedded函數在演示數據集(即與癌癥相關的成纖維細胞(CAFs) scRNAseq數據)上輕松地顯示初始的低維坐標: from sklearn.manifold import SpectralEmbeddingmodel = SpectralEmbedding(n_components = 2, n_neighbors = 50)se = model.fit_transform(np.log(X_train + 1))plt.figure(figsize=(20,15))plt.scatter(se[:, 0], se[:, 1], c = y_train.astype(int), cmap = 'tab10', s = 50)plt.title('Laplacian Eigenmap', fontsize = 20)plt.xlabel('LAP1', fontsize = 20)plt.ylabel('LAP2', fontsize = 20)plt.show()
最后，UMAP使用隨機梯度下降(SGD)代替常規梯度下降(GD)，如tSNE / FItSNE，這既加快了計算速度，又減少了內存消耗。 現在讓我們簡要地討論一下為什么他們說tSNE只保留數據的局部結構。可以從不同的角度來理解tSNE的局部性。首先,我們有σ參數Eq。(1)本地數據點集這樣互相“感覺”。因為成對歐幾里得距離衰減指數的概率,在小的σ值,它基本上是零遙遠的點(大型X)和快速增長僅為最近的鄰居(小X)。相比之下,在大的σ,遙遠而近點的概率成為限制可比和σ→∞,概率就等于1為所有任何一對點之間的距離,即成為等距點。plt.figure(figsize=(20, 15))x = np.linspace(0, 10, 1000)sigma_list = [0.1, 1, 5, 10]
for sigma in sigma_list: my_prob = lambda x: np.exp(-x**2 / (2*sigma**2)) plt.plot(x, my_prob(x))plt.gca().legend(('$\sigma$ = 0.1','$\sigma$ = 1', '$\sigma$ = 5', '$\sigma$ = 10'), fontsize = 20)plt.title('High-Dimensional Probability of Pairwise Euclidean Distances', fontsize = 20)plt.xlabel('X', fontsize = 20); plt.ylabel('P(X)', fontsize = 20)plt.show()
有趣的是,如果我們擴大成對歐幾里得距離的概率高維度成泰勒級數在σ→∞,我們會在第二近似冪律:
關于兩兩歐幾里得距離的冪律類似于多維定標法(MDS)的成本函數，MDS是通過保存每對點之間的距離來保存全局距離，而不管它們是相距很遠還是很近。一個可以解釋這個大的σtSNE遠程數據點之間的相互作用,所以是不完全正確的說tSNE只能處理當地的距離。然而，我們通常會受到perplexity有限值的限制，Laurens van der Maaten建議perplexity的取值范圍在5到50之間，盡管在局部信息和全局信息之間可能會有一個很好的折衷，那就是使用平方根≈N^(1/2)來選擇perplexity，其中N為樣本量。相反的極限,σ→0,我們最終的極端“局部性”高維概率的行為類似于狄拉克δ函數的行為。
另一種理解tSNE“局部性”的方法是考慮KL-divergence函數。假設X是高維空間中點之間的距離Y是低維空間中點之間的距離
根據kl -散度的定義：
方程(9)的第一項對于X的大小都是趨近于0的，對于X的大小也是趨近于0的，因為指數趨近于1，而log(1)=0。對于大X，這一項仍然趨近于0因為指數前因子趨近于0的速度快于對數趨近于負無窮。因此，為了直觀地理解kl散度，只考慮第二項就足夠了:
這是一個看起來很奇怪的函數，讓我們畫出KL(X, Y) import numpy as npimport matplotlib.pyplot as pltfrom matplotlib import cmfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure(figsize=(20, 15))ax = fig.gca(projection = '3d')# Set rotation angle to 30 degreesax.view_init(azim=30)
X = np.arange(0, 3, 0.1)Y = np.arange(0, 3, 0.1)X, Y = np.meshgrid(X, Y)Z = np.exp(-X**2)*np.log(1 + Y**2)
# Plot the surface.surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,linewidth=0, antialiased=False)ax.set_xlabel('X', fontsize = 20)ax.set_ylabel('Y', fontsize = 20)ax.set_zlabel('KL (X, Y)', fontsize = 20)ax.set_zlim(0, 2.2)ax.zaxis.set_major_locator(LinearLocator(10))ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))fig.colorbar(surf, shrink=0.5, aspect=5)plt.show()
這個函數的形狀非常不對稱。如果點在高維度X之間的距離很小,指數因子變成1和對數項行為日志(1 + Y ^ 2)這意味著如果Y是在低維空間的距離大,將會有一個大的懲罰,因此tSNE試圖減少Y在小X為了減少罰款。相反，對于高維空間中的長距離X, Y基本上可以是0到∞之間的任何值，因為指數項趨于0，并且總是勝過對數項。因此，在高維空間中相距遙遠的點，在低維空間中可能會相互靠近。因此，換句話說，tSNE并不能保證高維空間中相距較遠的點在低維空間中會保持較遠的距離。然而，它確實保證了在高維空間中相鄰的點在低維空間中保持相鄰。所以tSNE不是很擅長遠距離投射至低維,所以它只保留本地數據結構提供了σ不去∞。 與tSNE不同，UMAP使用交叉熵(CE)作為成本函數，而不是KL-divergence
這導致了地方-全球結構保護平衡的巨大變化。在X的小值處，我們得到了與tSNE相同的極限，因為第二項由于前因子和對數函數比多項式函數慢的事實而消失:
因此，為了使懲罰規則最小化，Y坐標必須非常小，即Y→0。這與tSNE的行為完全一樣。但是，在大X的相反極限，即X→∞時，第一項消失，第二項的前因子為1，得到:
這里，如果Y很小，我們會得到一個很大的懲罰，因為Y在對數的分母上，因此，我們鼓勵Y很大，這樣，對數下的比率就變成了1，我們得到零懲罰。因此，我們在X→∞處得到Y→∞，所以從高維空間到低維空間的整體距離保持不變，這正是我們想要的。為了說明這一點，讓我們繪制UMAP CE成本函數: import numpy as npimport matplotlib.pyplot as pltfrom matplotlib import cmfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure(figsize=(20, 15))ax = fig.gca(projection = '3d') # Set rotation angle to 30 degreesax.view_init(azim=30)
X = np.arange(0, 3, 0.001)Y = np.arange(0, 3, 0.001)X, Y = np.meshgrid(X, Y)Z = np.exp(-X**2)*np.log(1 + Y**2) + (1 - np.exp(-X**2))*np.log((1 + Y**2) / (Y**2+0.01))
# Plot the surface.surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,linewidth=0, antialiased=False)ax.set_xlabel('X', fontsize = 20)ax.set_ylabel('Y', fontsize = 20)ax.set_zlabel('CE (X, Y)', fontsize = 20)ax.set_zlim(0, 4.3)ax.zaxis.set_major_locator(LinearLocator(10))ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
# Add a color bar which maps values to colors.fig.colorbar(surf, shrink=0.5, aspect=5)plt.show()
在這里，我們可以看到圖的“右”部分看起來與上面的kl散度曲面非常相似。這意味著在X低的時候，為了減少損失，我們仍然想要Y低。然而，當X很大時，Y的距離也要很大，因為如果它很小，CE (X, Y)的損失將是巨大的。記住，之前，對于KL (X, Y)曲面，在X很大的情況下，我們在高Y值和低Y值之間沒有差別，這就是為什么CE (X, Y)代價函數能夠保持全局距離和局部距離。 我們知道UMAP是速度比tSNE擔憂)時大量的數據點,b)嵌入維數大于2或3,c)大量環境維度的數據集。在這里,讓我們試著了解UMAP要優于tSNE來自于數學和算法實現。 tSNE和UMAP基本上都包含兩個步驟： 建立一個圖在高維度的帶寬和計算指數概率,σ,使用二進制搜索和固定數量的最近的鄰居需要考慮。 通過梯度下降優化低維表示。第二步是算法的瓶頸，它是連續的，不能是多線程的。由于tSNE和UMAP都執行第二步，所以并不清楚為什么UMAP比tSNE更有效。
然而，我注意到UMAP的第一步比tSNE快得多。這有兩個原因： 首先，我們去掉了最近鄰數量定義中的log部分，即不像tSNE那樣使用全熵:
由于在算法上，對數函數是通過泰勒級數展開來計算的，而且在線性項前面加上一個對數前因子實際上并沒有增加多少，因為對數函數比線性函數慢，所以最好完全跳過這一步。 第二個原因是我們忽略了高維概率的歸一化，即式(1)中用于tSNE的高維概率歸一化。可以說，這一小步實際上對演出產生了戲劇性的影響。這是因為求和或積分是一個計算開銷很大的步驟。
接下來，UMAP實際上在第二步中也變得更快了。這種改善也有幾個原因: 采用隨機梯度下降法(SGD)代替常規梯度下降法(GD)，如tSNE或FItSNE。這提高了速度，因為對于SGD，您可以從一個隨機樣本子集計算梯度，而不是像常規GD那樣使用所有樣本子集。除了速度之外，這還減少了內存消耗，因為您不再需要為內存中的所有樣本保持梯度，而只需要為一個子集保持梯度。 我們不僅跳過了高維和低維概率的標準化。這也省去了第二階段(優化低維嵌入)的昂貴求和。 由于標準的tSNE使用基于樹的算法進行最近鄰搜索，由于基于樹的算法隨著維數的增加呈指數增長，所以生成超過2-3個嵌入維數的速度太慢。這個問題在UMAP中通過去掉高維和低維概率的標準化得到了解決。 增加原始數據集的維數，我們引入稀疏性，即我們得到越來越多的碎片化流形，即有時存在稠密區域，有時存在孤立點(局部破碎流形)。UMAP解決這個問題通過引入本地連接ρ參數融合在一起(某種程度上)通過引入自適應稀疏區域指數內核,考慮了本地數據連接。這就是為什么UMAP(理論上)可以處理任意數量的維度，并且在將其插入主維度縮減過程之前不需要預降維步驟(自動編碼器、PCA)的原因。
在這篇文章中，我們了解到盡管tSNE多年來一直服務于單細胞研究領域，但它有太多的缺點，如速度快、缺乏全球距離保存。UMAP總體上遵循了tSNE的哲學，但是引入了一些改進，例如另一個成本函數和缺少高維和低維概率的標準化。
|
