降維打擊-觀察高維世界

01-24

10-03

為了避免歧義，將第二小節「互相距離相同」改為「兩兩距離相同」

較早前寫的一篇文章，這裡完善了一下發出來，希望能夠有所啟發。

我們知道1維的世界最多有2個點兩兩距離相同，2維的世界最多有3個點，3維有4個。那麼4維呢？數學上的直觀可以告訴我們是5個，並且可以證明出來。但是，你能想像這種情況嗎？–不能，因為我們天生如此。

你可能覺得這個沒有什麼意思——各種數學工具做的足夠好了。但是，對於一些複雜而又不平凡的情況，不是用簡單的數學可以解決的。比如如果我們把手寫數字的圖片想像成一個向量–向量每一項代表圖像中每一點的像素值，那麼即使對於mnist這種最簡單手寫數字圖片的數據集之一，向量的維度也可以達到 28*28 = 784。

這是MNIST數據集的幾個例子：

我們這樣表示向量：

先把它表示為一個矩陣（應該28*28的，沒有畫全）

然後攤平即可

整個數據集構成了784維空間中向量（或者點）的集合，並且難以用一個數學公式描述出它的特性。

這個時候，為了直觀表示出高維中數據的「樣子」，我們需要降維——即把高維空間投射到3維以內。

降維就其本身來說是非常平凡的東西。最簡單的方法是「投影「，就是沿著一個軸把一個方向的數據的變化消除（比如把3維的一個方向上的東西壓平到平面上成為2維），或者說將向量的某個分量全部置0。這樣每次投影可以減少一個緯度。但是對於784維，至少需要781次投影，這樣幾乎確信自己很難發現什麼（一般來說4維以上的進行投影就基本無效了）。更麻煩的問題是，我們不知道沿著哪個方向投影更好。

我們需要更加漂亮的方式來降低維度，同時維持高維度的」感覺「。

首先我們有一種數學方法。這種方式叫PCA（主成分分析）。PCA的數學定義是：一個正交化線性變換，把數據變換到一個新的坐標系統中，使得這一數據的任何投影的第一大方差在第一個坐標（稱為第一主成分）上，第二大方差在第二個坐標（第二主成分）上，依次類推。從資訊理論角度考慮主成分分析在降維上是最優的（證明略，主要通過信息熵），它解決了投影問題。但是實際過程中我們可以發現它做的並不總是令人滿意，比如下圖是它對MNIST的784維空間的一個投影(選了一部分數據，更多數據就糟到看不清了)，不同顏色代表不同數字。我們可以看到很多數字「糾纏」到一起（除了左下角的1和右下角的0），和我們對數字容易區分的認知相差很遠。

你可能感到不理解：為什麼資訊理論上面可以證明最優的在實際中並不是最優？這是因為，「最優」的前提是對所有情況發生概率一致的基礎上得出的，對有先驗知識的東西並不成立。這個再次凸顯了「A more general theory is more useless」的道理。

比如，PCA 處理下面雜訊圖片的方式和處理上面的沒有什麼不同。它不關心數據本身，而僅僅是一個方向上面的方差，所以是一個「過大」的理論（但對於方差顯著的數據很有效）。

那麼我們有沒有其他方法呢？

註：

PCA 常被稱為線性降維方法，因為它對數據進行了線性變換。另外一種常用的線性降維方法是 LDA (Linear Discriminant Analysis)，下面正文要談的都是非線性的。PCA也可以使用 kernel trick 變成非線性的。
Laplacian Eigenmaps（由求圖的 Laplace 矩陣的特徵值而得名）也是一種線性降維方法。
LLE (Locally Linear Embedding)，局部線性嵌入演算法，是一種非線性演算法。LLE演算法認為每一個數據點都可以由其近鄰點的線性加權組合構造得到。
PCA, LDA, Laplacian Eigenmaps, LLE 都屬於譜（spectral）方法（涉及矩陣特徵向量的方法）。
線性方法計算效率是非常高的，非線性演算法對於大規模數據需要考慮隨機遊走(random walk)等優化，否則計算代價非常高。

距離表示

仔細想想我們會發現，表示出數據的空間結構並不是一件完全沒有規律可遵循的事情。首先，旋轉、平移和翻轉不會改變數據本身的結構，對於高維空間也是如此。並且結構還明顯受點的個數的影響–3個點對於高於2維的空間是沒有特殊意義的，我們需要一種東西可以保證不受維數增加的影響。

因此，我們開始考慮用點之間的距離（這裡使用歐氏距離）代替點的坐標 – 距離不受旋轉和平移的影響，並且只和點數有關。

但是，距離表示會不會丟失某些結構信息呢？不會，我們以2維為例構造式的說明一下：

只有一個距離的時候當然隨便什麼角度和位置排放點都不會影響結構：

假設不是退化到一維的情況（即存在不共線的點），它們三個點的三個邊一定構成唯一形狀的三角形

下面任何一個點將和這三個點構成3條邊，也就是3個方程，足以確定2維坐標的兩個分量。否則那個點不在二維平面上。

其他維度類似。所以各點間的距離確實完整保留了空間結構，同時又去除了不必要的轉動和平移翻轉。同時，距離本身不直接包含維度信息，不會受到冗餘維度的影響。

但是，正是因為距離本身不直接包含維度，所以我們似乎沒有什麼特別明確的演算法。因此我們退一步，把問題轉化為優化問題的形式，以備我們有更靈活的選擇：

$min_{D_Y} E(D_X, D_Y)$

其中 $D_X, D_Y$ 分別是原空間和當前空間中的邊。E（能量函數，損失函數）是一種評估它們「結構相似」的函數。

然後我們試圖尋找一個策略，使得它們「最為相似」。這樣比直接求 $D_Y$ 有更多選擇的餘地。

但困難之處還在於，你必須保持低維度空間的距離是合法的（比如三維中一共5點10條邊，長度都相等是無解的）。所以我們需要間接的方法：

設 $X$ 是原空間的點（向量）， $Y$ 是低維空間的點， $D_X, D_Y$ 分別是它們對應的邊。我們尋求一個 $Y$ 最小化 $E$ ：

$min_Y E(D_X, D_Y)$

這樣一來在於，我們將難以直接求解 $Y$ ，因為 $Y$ 的距離是關於 $Y$ 的相當複雜的函數（帶 $n*dim$ 個自變數， $n*(n-1)/2$ 個因變數，還有平方，根號），再加上 $E$ 的與 $D_X$ 相關的複雜運算，使得我們基本不可能有個可以接受的公式直接根據 $E$ 求得 $Y$ 。

不過不用害怕。這類問題的解法早已熟知，那就是採用各種數值優化演算法，對於多數使用梯度下降就足夠了，這裡暫時不做更多介紹。

所以下面的關建是如何選擇 $E$ （能量函數，損失函數）

基於距離相似的 E

MDS (Multidimensional scaling)

首先我們想到的是，顯然一般情況下我們沒有辦法讓 $D_Y$ 和 $D_X$ 完全相同，否則低維空間無法容納這樣的 $D_Y$ 。這讓我們想到了直線擬合：對於大部分實際數據都無法真正過一條直線，因此我們有最小二乘法，即使得直線和數據點之間的平方誤差最小（這是少數我們能夠有解析解的情況，對於不能轉化為直線的一些複雜的情況我們只能採用優化的方式）。按照最小二乘法的思想，我們很容易設計以下的 $E$ （按照求和約定省略了部分標記，但保留了求和符號。下同）：

$E_{MDS} = frac 1{sum D_{X_i}^2}sum (D_{X_i}-D_{Y_i})^2$

前面的係數是正則化係數，目的是為了保證 EE 不受距離加倍的影響（同時把所有距離乘上一個正因子不應該改變結果）

這種方法屬於所謂 Multidimensional scaling

結果如下：

可以看到一些團簇，但是仍然不能讓人滿意。

我們注意到，問題在於一些理應分開的點變的像一團粥一樣，Multidimensional scaling 同等的對待所有的距離，這是和PCA一樣的問題。因此我們希望能夠更加強調一些局部結構，讓它們更加緊密一些。

Sammon Mapping

我們稍稍修改一下上面的 $E$ ，把一部分 $D_X$ 移入到後面參加運算。

$E_{sammon} = frac 1{sum D_{X_i}}sum frac{(D_{X_i}-D_{Y_i})^2}{D_{X_i}}$

這就是所謂 Sammon Mapping。它的效果是強調更短的邊（在分母上，越小就意味著上面的平方項會帶來更大的損失，從而約束上面的平方項更加小）

然而它的效果似乎並不是那麼明顯：

然而大量的事實確實證實 Sammon Mapping 在比較簡單的情況下很有效。為什麼這裡就不行了呢？

因為，這是784維的真實數據，這麼高的維數幾乎是種詛咒。

Isomap

之前我們採用距離是歐幾里德距離，它是平坦空間的最短距離。

Isomap 把距離推廣為了「測地距離」。「測地距離」是兩點之間最短路徑長度（離散空間）或者某種加權的測度（連續的空間）。一般Isomap會先把問題轉化為離散的處理，這裡不舉例子了。

附：維數災難

我們來做一個有趣的實驗：

對於k維，立方體的頂點數是 $2^k$ ，原因是在單位正交坐標系中，k維立方體的頂點可以用k個0/1構成的坐標描述，比如一個正方形頂點為 $(0,0),(0,1),(1,0),(1,1)$ 。這種廣義的正方體有很多奇妙的性質，比如相鄰的點僅有一個分量不同，可以從一個點出發不重複沿著邊走遍所有點等等。

然後我們計算每個點到其他點的距離的方差（寫個程序就行了，注意一些技巧）。結論是方差隨著維數增多緩慢減小，並穩定在 0.125 左右。

這意味著什麼？我們知道，一個分布的方差是由數據分布的「形狀」決定的，越彌散的分布有越大的方差。方差不變意味著分布形狀基本不變，但是引入了過多的數據，這樣一來一個分布局部區間必然包含了更多的數據，也就是更多距離相似的邊，雖然立方體這個數學上的模型沒有變化。

比如，對於2維數據，對角線和邊長相差 $sqrt{2}-1approx 0.414$ ，但是對於1000維的數據，存在 499,500 個長 $sqrt{998}$ 和 166,167,000 個長 $sqrt{997}$ 的邊，而它們的長度僅僅相差 0.016，而相對相差只有 0.0005。

這個簡單的例子就是所謂「維數災難」的一種形式，即隨著維數增大，邊之間的「解析度」會降低。這可能就是 Sammon Mapping 沒有取得很好結果的原因，即它除的分母太相似了，沒有過多的區分度。

基於物理體系模擬的 E

之前我們公式中的主項， $sum (D_{X_i}-D_{Y_i})^2$ 是否有什麼物理意義呢？

如果把現有空間中的邊看作彈簧，那麼上面的公式就是勁度係數為1的彈簧體系的彈性勢能。其中對於彈簧 $D_{Y_i}$ ， $D_{X_i}$ 是它的平衡位置。

我們觀察到 Sammon Mapping 試圖靠攏相鄰點的做法失敗了，那麼我們是不是可以考慮推遠不相鄰的點？

做法很簡單，就是引入一個勢能場，比如讓所有點都帶上單位同種電荷。

這樣一來能量函數（名副其實）就是

$E_{ForceField} = sumfrac1{D_{Y_i}} + sum k(D_{X_i}-D_{Y_i})^2$

其中k是勁度係數。(註：靜電勢能選擇反比形式，是因為考慮到三維中平方受力的穩定性較好。對於不同維度，維持靜電場為有源無旋場根據高斯定理會導致不同的勢能形式。)

效果如下：

嗯，好像沒有什麼改進。原因同 Sammon Mapping，前面的靜電勢能項因為大家距離都差不多而沒有起到實質作用。

那麼我們有沒有改進方案呢？在我們腦海中，一個彈簧和斥力的體系會自己根據受力平衡「布局」到一個好的位置，這也是布圖的一種常見演算法。

它的能量函數是 (k是是勁度係數， $D_0$ 是彈簧的平衡位置)：

$E_{ForceField} = sumfrac1{D_{Y_i}} + sum_{Y_i^* in edges} k(D_{Y_i^*}-D_0)^2$

你可能感到不可思議的是， $D_X$ 到哪裡去了？我們怎麼可能不需要原來的距離就得到結果呢？

答案是， $D_X$ 用來確定edges，也就是保留哪些邊。

注意，之前我們是把所有的邊都納入到能量體系當中，而現在我們需要根據 $D_X$ 選擇一些作為「彈簧」。

選擇的方法是使用古老的 KNN (k-nestest neighbors)演算法，也就是選擇每個點最近的k個邊。

這樣我們就有

$E_{KNN-ForceField} = sumfrac1{D_{Y_i}} + sum_{Y_i^* in knnedges (D_X)} k(D_{Y_i^*}-D_0)^2$

knn 的 k = 4 的時候，

k = 3

嗯，總體上要好很多。主要貢獻應該歸於用KNN預處理了 DXDX，KNN僅僅選擇最短的幾個邊，即使各個邊長度差距不大，這樣大大增加了區分度。之前的方法主要通過把 $D_{X_i}$ 放在分母，區分度是不夠的。

但是引入KNN的方法把問題變成了黑盒：k 對數據本身的意義是什麼？這個不是很清楚。另外邊要不選中，要不不選中，似乎也是一種很極端的做法。

基於概率分布的 E

我們希望有一種演算法，能夠一定程度上克服使用 KNN 的缺陷。

一種方法是用概率分布取代「是否」這種決定式的方法。這類方法稱為 SNE （隨機鄰居嵌入）。

首先，它將距離轉換為概率分布。概率分布中概率比較大的表示更有可能「靠在一起」。

附：將有限個連續數量值轉換為概率分布

首先，由於數量值是連續的（比如實數），所以我們需要一個連續的概率分布的概率密度函數(pdf)。

我們直接先把概率密度函數作用上去

$p_1,p_2,p_3...p_n = f(d_1),f(d_2),f(d_3)...f(d_n)$

這個時候得到了一些離散的概率值，因為數量值的個數是有限的。

接下來我們歸一化這些概率值：

$p_i^* = frac{p_i}{sum p_j}$ 現在 p?p? 變成了一個標準的離散的概率分布，它是連續分布的一種近似。

$square$

然後我們評測原來和現在概率的「相似性」。

對於有重疊的概率分布，一種非常有效而且有很強理論支持的用於表徵相似的測度叫 KL divergence （KL散度），一般用記號 ${ m KL}left(cdot|cdot ight)$ 表示。它也叫相對熵。

所以形式表示為

$E_{SNE} = { m KL}left(K_X(D_X)|K_Y(D_Y) ight)$

符號意義同上文。其中K表示某種映射，將邊集合轉換為概率分布（將有限個連續數量值轉換為概率分布）。

SNE - Part I 將原始高維空間中的每個點相關的距離（鄰居距離）轉化為概率分布

我們注意到，先前 Sammon Mapping 這些效果不好的原因是對於距離的區分度不是特別好。因此我們考慮通過控制分布的形狀來控制區分度。

我們選擇的是均值為0的高斯分布（正態分布）。

形式如下：

$N(0, sigma) = frac1{sqrt{2pi}sigma}expleft(-frac{x^2}{2sigma^2} ight)$

高斯分布是在方差相等的分布中信息熵最大的，其信息熵為 $lnleft(sigmasqrt{2pi e} ight)$ 。在分布未知的情況下根據中心極限定理我們非常有理由選擇 Gaussian distribution，因為它是對未知分布的良好近似。

然後我們把這種分布應用到每個點到其他點的距離，得到一組歸一化分布 $g^{(i)}_sigma$ ，其中 i 代表第 i 個點， $sigma$ 是所選的 gaussian 分布的方差。（注意，這裡是每個點而言的分布，而不是對所有邊的分布）

現在問題是，我們如何選擇 $sigma$ ？

我們用信息熵評價分布的區分度（即 ${ m H}(g^{(i)}_sigma)$ ），熵越高則區分度越低，來控制 $sigma$ 的取值，則可以證明 ${ m H}(g^{(i)}_sigma)$ 是關於 $sigma$ 的單調遞增函數。

這個結論並不奇怪。 $sigma$ 代表方差，方差越小則分布越「尖銳」，區分度越高；否則分布越平坦，區分度越低。

但是直接指定信息熵看上去不像是一個良好的做法。因為對於連續分布本身，

${ m H}left(N(0,sigma) ight) = lnleft(sigmasqrt{2pi e} ight)$

即

$sigma propto e^{H(0,sigma)}$

這意味著 σσ 對信息量過於敏感，不便於我們調節。

我們定義 $Perplexity = 2^{H}$ ，則 $sigma propto Perplexity$ ，相應的離散距離分布也大致滿足這個關係。

（這裡我們不區分自然對數，2對數和常用對數，因為僅僅影響熵的「單位」而已）

所以，總結一下，就是

$Find sigma, s.t. { m H}(g^{(i)}_sigma) = ln (Perplexity)$

尋找我們使用二分法即可（由於 Perplexity 隨著 $sigma$ 單調增），因為函數太過複雜不適合使用牛頓迭代法等需要導數的方式。

SNE - Part II 將原始高維空間中鄰居距離轉化為所有邊距離的概率分布

到這裡仍然沒有結束。因為我們僅僅求得關於每個點的鄰居距離的概率分布（對有n個點的數據，就有n個分布），目的是控制距離的顯著程度。

我們把這些分布排列成矩陣，其中每行都是一個歸一化的概率分布。

$left[ egin{matrix} g^{(1)} \ g^{(2)} \ vdots \ g^{(n)} end{matrix} ight] = left[ egin{matrix} g^{(1,1)} & g^{(1,2)} & cdots & g^{(1,n)}\ g^{(2,1)} & g^{(2,2)} & cdots & g^{(2,n)}\ vdots & vdots & ddots & vdots\ g^{(n,1)} & g^{(n,2)} & cdots & g^{(n,n)} end{matrix} ight] % ag{3}$

下面我們需要將這n個獨立的分轉化為關於所有邊的分布。這樣一來 $g^{(i,j)}+g^{(j,i)}$ 代表點 $i$ 和 $j$ 的距離相應的概率。

首先我們進行歸一化。由於有n個已經歸一化的分布，所以為了使得總體分布歸一化，我們把矩陣中每一項除以 $n$ 即可。

然後進行對稱化。由於 $g^{(i,j)}, g^{(j,i)}$ 對應同一個距離，不應該有差異，所以我們讓它們各取一半，即

$p^{(i,j)} = p^{(j,i)} = frac12left(g^{(i,j)}+g^{(j,i)} ight)$

最終我們把距離分布轉化為了概率分布，即

$left[ egin{matrix} d^{(1,1)} & d^{(1,2)} & cdots & d^{(1,n)}\ d^{(2,1)} & d^{(2,2)} & cdots & d^{(2,n)}\ vdots & vdots & ddots & vdots\ d^{(n,1)} & d^{(n,2)} & cdots & d^{(n,n)} end{matrix} ight] ightarrow left[ egin{matrix} p^{(1,1)} & p^{(1,2)} & cdots & p^{(1,n)}\ p^{(2,1)} & p^{(2,2)} & cdots & p^{(2,n)}\ vdots & vdots & ddots & vdots\ p^{(n,1)} & p^{(n,2)} & cdots & p^{(n,n)} end{matrix} ight]$

我們將最後得到的概率矩陣記為 $P_X$

SNE - Part IV 將低維度維空間中距離轉化為概率分布

低維度中我們不應採用先前那樣的帶參數的概率密度函數，原因是如果使用就代表我們「先驗」地定義結果的性質而不是輸入的性質，這是一種「自欺欺人」的做法，即我們不應該為了更好的結果而改變結果的形式。

傳統的 SNE 使用的是 $N(0, 1)$ ，即標準正態分布。但是事實證明使用自由度為1的（0均值）t 分布效果更好，對應的方法稱為 t-SNE。

t 分布是統計學三大分布（ $chi^2$ 分布，t 分布，F 分布）之一，又稱為學生分布 (和一段署名為「學生」的匿名發表軼事有關）。自由度為1的密度函數為

$t(X) = frac1{pi(1+x^2)}$ 它比標準正態分布更加平緩（紅色為t分布）：

為什麼 t 分布效果更好？一個可能的原因在於 t 分布本身的性質。正態分布用於估計方差為σσ 的數據的均值，而 t 分布用於估計方差未知的數據的均值，也就是 t 分布更加低維數據沒有先驗方差的情況，而標準正態分布假定了方差為1。

我們使用相同的方法將 t 分布應用於 $D_Y$ 並且歸一化得到所有邊的概率分布，記為 $Q_Y$ 。

SNE - Part V 優化

最終的E為

$E_{t-SNE} = { m KL}(P_X|Q_Y), with certain perplexity$

優化目標為

$min_Y { m KL}(P_X|Q_Y), with certain perplexity$

t-SNE 的優化目標相對之前介紹的內容複雜很多，而且是非凸的（即存在大量的局部最優值），往往需要加上絕熱退火、動量下降等方法。

下面是結果：

可以看出是相當不錯的。

下面的圖是 t-SNE 論文中的對比結果，t-SNE是碾壓的勝利。

t-SNE 也有一些問題，比如速度慢，結果因為有局部最優所以有點隨機，對於非常一般的一些數據處理不好等等。但是對於大型的真實數據集和對深度神經網路的分析上面 t-SNE 優勢很大，所以常常用於深度學習的結果分析中。

t-SNE 的 perplexity 是必要的參數嗎？

相信 t-SNE 讓人不滿意的地方之一在於有個 perplexity 參數。這個參數是必要的嗎？

幾乎可以說是必要的，而且很有意義，甚至是 t-SNE 優勢之一。

我們考慮下面的情況：

我們把一條直線上面的點繞折構成一個「平面」：

然後把平面像壽司卷那樣捲成壽司棒（下面的圖切開了 T T）

側視圖

然後我們延長這個「棒子」，把它作為「直線」，繼續如上的構造。

問題是，在二維空間它應該呈現出什麼樣子才能體現它的性質？一條線？一個面？一個圓餅？似乎都有道理，因為在於尺度的大小。上述構造的實質是不斷把各個維度的「捲曲」作為信息加入到數據中，這些信息對應了多級結構，我們需要獲取這些信息的時候就必須考慮到不同尺度的結構。

t-SNE 的 perplexity 恰恰控制的是尺度，原因如下：

由於採用了 gaussian 分布，所以我們知道大約68%的數據集中在一個標準差 $sigma$ 內，也就是大約在一個 perplexity 的範圍內。

perplexity 非常小的時候，在計算鄰居距離的過程中，如下圖所示，剛好能夠覆蓋最小尺度的「鄰居」。

在全局的歸一化和對稱化完成後，實質上形成了這樣的分布：

這樣 t-SNE 會展現「直線」的姿態（一級結構）。

如果perplexity更大一點就會像下面一樣

但是側面覆蓋範圍還不夠大

所以最終會得到展開的平面（二級結構）

再大就足以覆蓋三級結構(棒狀)

所以說 perplexity 一定程度上決定了我們看到的結構尺度。

Big Data talks

中心極限定理告訴我們，對於獨立同分布的數據，數據樣本越多則能越精確地估計數據的均值。統計學三大分布也告訴我們這個結論對估計方差比、均值差等也成立。

我們可以直覺上感到，隨著數據增多，對於數據其他特性的估計也會變的更好，包括數據的「形狀」。

之前採用的例子都是 1,000 個 MNIST 數字圖像。我們試試用 2,000 個：

使用2000個圖像的 Sammon mapping，改進有限。

使用2000個圖像的 KNN-ForceField（k=3），進步挺大

使用2000個圖像的 t-SNE (perplexity = 40)，進步更大:

t-SNE 處理的 MNIST 全體大合影（60,000個），效果非常好

可以看出 t-SNE 甚至保留了很多規律，比如右側的1的位置和它的形態明顯相關

另外圖中也暗示了某些數字難以區分。

維度

之前我們都是把高維度投射到2維的空間中，如果是3維那麼如何？下面我們先用 1000 張圖片作測試。

Sammon Mapping，有提升

KNN-ForceField (k=3)

t-SNE 效果也不錯

將數目增加到 2000個

KNN-ForceField (k=3) 側視圖

t-SNE 側視圖1

t-SNE 側視圖2

效果是不錯的

一般情況下，降維到三維要優於到二維，這個也符合常理。但是三維需要「轉著看」，所以論文和著作中更偏向二維的。

章節附註

可視化使用了 d3js 和 threejs，其中 d3js 用於控制顏色和選擇控制項，threejs 是個使用 WebGL 進行渲染的庫。

另外參見了 http://colah.github.io/ 和 Andrej Karpathy 的一些blog和項目。很高興我回過來看的時候發現很多想法是和他們想到一起了。

由於快登機了，就先寫這麼多~