漫談 Clustering

http://blog.pluskid.org/?p=17

1. k-means

   好久沒有寫 blog 了，一來是 blog 下線一段時間，而租 DreamHost 的事情又一直沒弄好；二來是沒有太多時間，天天都跑去實驗室。現在主要折騰 Machine Learning 相關的東西，因為很多東西都不懂，所以平時也找一些資料來看。按照我以前的更新速度的話，這么長時間不寫 blog 肯定是要被悶壞的，所以我也覺得還是不定期地整理一下自己了解到的東西，放在 blog 上，一來梳理總是有助於加深理解的，二來也算共享一下知識了。那么，還是從 clustering 說起吧。

   Clustering 中文翻譯作"聚類"，簡單地說就是把相似的東西分到一組，同 Classification（分類）不同，對於一個classifier ，通常需要你告訴它"這個東西被分為某某類"這樣一些例子，理想情況下，一個 classifier 會從它得到的訓練集中進行"學習"，從而具備對未知數據進行分類的能力，這種提供訓練數據的過程通常叫做 supervised learning (監督學習)，而在聚類的時候，我們並不關心某一類是什么，我們需要實現的目標只是把相似的東西聚到一起，因此，一個聚類算法通常只需要知道如何計算相似度就可以開始工作了，因此 clustering 通常並不需要使用訓練數據進行學習，這在 Machine Learning 中被稱作 unsupervised learning (無監督學習)。

   舉一個簡單的例子：現在有一群小學生，你要把他們分成幾組，讓組內的成員之間盡量相似一些，而組之間則差別大一些。最后分出怎樣的結果，就取決於你對於"相似"的定義了，比如，你決定男生和男生是相似的，女生和女生也是相似的，而男生和女生之間則差別很大，這樣，你實際上是用一個可能取兩個值"男"和"女"的離散變量來代表了原來的一個小學生，我們通常把這樣的變量叫做"特征"。實際上，在這種情況下，所有的小學生都被映射到了兩個點的其中一個上，已經很自然地形成了兩個組，不需要專門再做聚類了。另一種可能是使用"身高"這個特征。我在讀小學候，每周五在操場開會訓話的時候會按照大家住的地方的地域和距離遠近來列隊，這樣結束之后就可以結隊回家了。除了讓事物映射到一個單獨的特征之外，一種常見的做法是同時提取 N 種特征，將它們放在一起組成一個 N 維向量，從而得到一個從原始數據集合到 N 維向量空間的映射——你總是需要顯式地或者隱式地完成這樣一個過程，因為許多機器學習的算法都需要工作在一個向量空間中。

   那么讓我們再回到 clustering 的問題上，暫且拋開原始數據是什么形式，假設我們已經將其映射到了一個歐幾里德空間上，為了方便展示，就使用二維空間吧，如下圖所示：

   

   從數據點的大致形狀可以看出它們大致聚為三個 cluster ，其中兩個緊湊一些，剩下那個松散一些。我們的目的是為這些數據分組，以便能區分出屬於不同的簇的數據，如果按照分組給它們標上不同的顏色，就是這個樣子：

   

   那么計算機要如何來完成這個任務呢？當然，計算機還沒有高級到能夠"通過形狀大致看出來"，不過，對於這樣的 N 維歐氏空間中的點進行聚類，有一個非常簡單的經典算法，也就是本文標題中提到的k-means 。在介紹 k-means 的具體步驟之前，讓我們先來看看它對於需要進行聚類的數據的一個基本假設吧：對於每一個cluster ，我們可以選出一個中心點 (center) ，使得該 cluster 中的所有的點到該中心點的距離小於到其他 cluster 的中心的距離。雖然實際情況中得到的數據並不能保證總是滿足這樣的約束，但這通常已經是我們所能達到的最好的結果，而那些誤差通常是固有存在的或者問題本身的不可分性造成的。例如下圖所示的兩個高斯分布，從兩個分布中隨機地抽取一些數據點出來，混雜到一起，現在要讓你將這些混雜在一起的數據點按照它們被生成的那個分布分開來：

   

   由於這兩個分布本身有很大一部分重疊在一起了，例如，對於數據點2.5來說，它由兩個分布產生的概率都是相等的，你所做的只能是一個猜測；稍微好一點的情況是2，通常我們會將它歸類為左邊的那個分布，因為概率大一些，然而此時它由右邊的分布生成的概率仍然是比較大的，我們仍然有不小的幾率會猜錯。而整個陰影部分是我們所能達到的最小的猜錯的概率，這來自於問題本身的不可分性，無法避免。因此，我們將 k-means 所依賴的這個假設看作是合理的。

   基於這樣一個假設，我們再來導出 k-means 所要優化的目標函數：設我們一共有 N 個數據點需要分為 K 個 cluster，k-means要做的就是最小化

   這個函數，其中 在數據點n被歸類到cluster k的時候為1，否則為0。直接尋找  和  來最小化  並不容易，不過我們可以采取迭代的辦法：先固定  ，選擇最優的  ，很容易看出，只要將數據點歸類到離他最近的那個中心就能保證  最小。下一步則固定 ，再求最優的 。將  對  求導並令導數等於零，很容易得到  最小的時候  應該滿足：

   亦即 的值應當是所有cluster k中的數據點的平均值。由於每一次迭代都是取到  的最小值，因此  只會不斷地減小（或者不變），而不會增加，這保證了k-means 最終會到達一個極小值。雖然k-means並不能保證總是能得到全局最優解，但是對於這樣的問題，像k-means 這種復雜度的算法，這樣的結果已經是很不錯的了。

   下面我們來總結一下k-means 算法的具體步驟：

2. 選定 K個中心   的初值。這個過程通常是針對具體的問題有一些啟發式的選取方法，或者大多數情況下采用隨機選取的辦法。因為前面說過 k-means 並不能保證全局最優，而是否能收斂到全局最優解其實和初值的選取有很大的關系，所以有時候我們會多次選取初值跑 k-means ，並取其中最好的一次結果。
3. 將每個數據點歸類到離它最近的那個中心點所代表的 cluster 中。
4. 用公式  計算出每個cluster 的新的中心點。
5. 重復第二步，一直到迭代了最大的步數或者前后的  的值相差小於一個閾值為止。

   按照這個步驟寫一個 k-means 實現其實相當容易了，在 SciPy 或者 Matlab 中都已經包含了內置的 k-means 實現，不過為了看看 k-means 每次迭代的具體效果，我們不妨自己來實現一下，代碼如下（需要安裝 SciPy 和 matplotlib） ：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91

   #!/usr/bin/python     from __future__ import with_statement import cPickle as pickle from matplotlib import pyplot from numpy import zeros, array, tile from scipy.linalg import norm import numpy.matlib as ml import random     def kmeans(X, k, observer=None, threshold=1e-15, maxiter=300): N = len(X) labels = zeros(N, dtype=int) centers = array(random.sample(X, k)) iter = 0     def calc_J(): sum = 0 for i in xrange(N): sum += norm(X[i]-centers[labels[i]]) return sum     def distmat(X, Y): n = len(X) m = len(Y) xx = ml.sum(X*X, axis=1) yy = ml.sum(Y*Y, axis=1) xy = ml.dot(X, Y.T)     return tile(xx, (m, 1)).T+tile(yy, (n, 1)) - 2*xy     Jprev = calc_J() while True: # notify the observer if observer is not None: observer(iter, labels, centers)     # calculate distance from x to each center # distance_matrix is only available in scipy newer than 0.7 # dist = distance_matrix(X, centers) dist = distmat(X, centers) # assign x to nearst center labels = dist.argmin(axis=1) # re-calculate each center for j in range(k): idx_j = (labels == j).nonzero() centers[j] = X[idx_j].mean(axis=0)     J = calc_J() iter += 1     if Jprev-J < threshold: break Jprev = J if iter >= maxiter: break     # final notification if observer is not None: observer(iter, labels, centers)     if __name__ == '__main__': # load previously generated points with open('cluster.pkl') as inf: samples = pickle.load(inf) N = 0 for smp in samples: N += len(smp[0]) X = zeros((N, 2)) idxfrm = 0 for i in range(len(samples)): idxto = idxfrm + len(samples[i][0]) X[idxfrm:idxto, 0] = samples[i][0] X[idxfrm:idxto, 1] = samples[i][1] idxfrm = idxto     def observer(iter, labels, centers): print "iter %d." % iter colors = array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) pyplot.plot(hold=False) # clear previous plot pyplot.hold(True)     # draw points data_colors=[colors[lbl] for lbl in labels] pyplot.scatter(X[:, 0], X[:, 1], c=data_colors, alpha=0.5) # draw centers pyplot.scatter(centers[:, 0], centers[:, 1], s=200, c=colors)     pyplot.savefig('kmeans/iter_%02d.png' % iter, format='png')     kmeans(X, 3, observer=observer)

   代碼有些長，不過因為用 Python 來做這個事情確實不如 Matlab 方便，實際的 k-means 的代碼只是 41 到 47 行。首先 3 個中心點被隨機初始化，所有的數據點都還沒有進行聚類，默認全部都標記為紅色，如下圖所示：

   

   然后進入第一次迭代：按照初始的中心點位置為每個數據點着上顏色，這是代碼中第 41 到 43 行所做的工作，然后 45 到 47 行重新計算 3 個中心點，結果如下圖所示：

   

   可以看到，由於初始的中心點是隨機選的，這樣得出來的結果並不是很好，接下來是下一次迭代的結果：

   

   可以看到大致形狀已經出來了。再經過兩次迭代之后，基本上就收斂了，最終結果如下：

   

   不過正如前面所說的那樣 k-means 也並不是萬能的，雖然許多時候都能收斂到一個比較好的結果，但是也有運氣不好的時候會收斂到一個讓人不滿意的局部最優解，例如選用下面這幾個初始中心點：

   

   最終會收斂到這樣的結果：

   

   不得不承認這並不是很好的結果。不過其實大多數情況下 k-means 給出的結果都還是很令人滿意的，算是一種簡單高效應用廣泛的 clustering 方法。

   Update 2010.04.25: 很多人都問我要 cluster.pkl ，我干脆把它上傳上來吧，其實是很容易自己生成的，點擊這里下載

6. k-medoids

   上一次我們了解了一個最基本的 clustering 辦法 k-means ，這次要說的 k-medoids 算法，其實從名字上就可以看出來，和 k-means 肯定是非常相似的。事實也確實如此，k-medoids 可以算是 k-means 的一個變種。

   k-medoids 和 k-means 不一樣的地方在於中心點的選取，在 k-means 中，我們將中心點取為當前 cluster 中所有數據點的平均值：

   並且我們已經證明在固定了各個數據點的 assignment 的情況下，這樣選取的中心點能夠把目標函數  最小化。然而在 k-medoids 中，我們將中心點的選取限制在當前 cluster 所包含的數據點的集合中。換句話說，在 k-medoids 算法中，我們將從當前 cluster 中選取這樣一個點——它到其他所有（當前 cluster 中的）點的距離之和最小——作為中心點。k-means 和 k-medoids 之間的差異就類似於一個數據樣本的均值 (mean) 和中位數 (median) 之間的差異：前者的取值范圍可以是連續空間中的任意值，而后者只能在給樣本給定的那些點里面選。那么，這樣做的好處是什么呢？

   一個最直接的理由就是k-means 對數據的要求太高了，它使用歐氏距離描述數據點之間的差異 (dissimilarity) ，從而可以直接通過求均值來計算中心點。這要求數據點處在一個歐氏空間之中。

   然而並不是所有的數據都能滿足這樣的要求，對於數值類型的特征，比如身高，可以很自然地用這樣的方式來處理，但是類別 (categorical) 類型的特征就不行了。舉一個簡單的例子，如果我現在要對犬進行聚類，並且希望直接在所有犬組成的空間中進行，k-means 就無能為力了，因為歐氏距離在這里不能用了：一只Samoyed 減去一只 Rough Collie 然后在平方一下？天知道那是什么！再加上一只 German Shepherd Dog 然后求一下平均值？根本沒法算，k-means 在這里寸步難行！

   在k-medoids 中，我們把原來的目標函數中的歐氏距離改為一個任意的dissimilarity measure 函數 ：

   

   Samoyed Rough Collie

   最常見的方式是構造一個 dissimilarity matrix  來代表 ，其中的元素  表示第i只狗和第 j只狗之間的差異程度，例如，兩只 Samoyed 之間的差異可以設為 0 ，一只 German Shepherd Dog 和一只 Rough Collie 之間的差異是 0.7，和一只 Miniature Schnauzer 之間的差異是 1 ，等等。

   除此之外，由於中心點是在已有的數據點里面選取的，因此相對於 k-means 來說，不容易受到那些由於誤差之類的原因產生的 Outlier 的影響，更加 robust 一些。

   扯了這么多，還是直接來看看 k-medoids 的效果好了，由於 k-medoids 對數據的要求比 k-means 要低，所以 k-means 能處理的情況自然 k-medoids 也能處理，為了能先睹為快，我們偷一下懶，直接在上一篇文章中的 k-means 代碼的基礎上稍作一點修改，還用同樣的例子。將代碼的 45 到 47 行改成下面這樣：

   45 46 47 48 49 50

   for j in range(k): idx_j = (labels == j).nonzero() distj = distmat(X[idx_j], X[idx_j]) distsum = ml.sum(distj, axis=1) icenter = distsum.argmin() centers[j] = X[idx_j[0][icenter]]

   可以看到 k-medoids 在這個例子上也能得到很好的結果：

   

   而且，同 k-means 一樣，運氣不好的時候也會陷入局部最優解中：

   

   如果仔細看上面那段代碼的話，就會發現，從 k-means 變到 k-medoids ，時間復雜度陡然增加了許多：在 k-means 中只要求一個平均值  即可，而在 k-medoids 中則需要枚舉每個點，並求出它到所有其他點的距離之和，復雜度為  。

   看完了直觀的例子，讓我們再來看一個稍微實際一點的例子好了：Document Clustering ——那個永恆不變的主題，不過我們這里要做的聚類並不是針對文檔的主題，而是針對文檔的語言。實驗數據是從 Europarl 下載的包含 Danish、German、Greek、English、Spanish、Finnish、French、Italian、Dutch、Portuguese 和 Swedish 這些語言的文本數據集。

   在 N-gram-based text categorization 這篇 paper 中描述了一種計算由不同語言寫成的文檔的相似度的方法。一個（以字符為單位的） N-gram 就相當於長度為 N 的一系列連續子串。例如，由 hello 產生的 3-gram 為：hel、ell 和 llo ，有時候還會在划分 N-gram 之前在開頭和末尾加上空格（這里用下划線表示）：_he、hel、ell、llo、lo_ 和 o__ 。按照 Zipf's law ：

   The nth most common word in a human language text occurs with a frequency inversely proportional to n.

   這里我們用 N-gram 來代替 word 。這樣，我們從一個文檔中可以得到一個 N-gram 的頻率分布，按照頻率排序一下，只保留頻率最高的前 k 個（比如，300）N-gram，我們把叫做一個"Profile"。正常情況下，某一種語言（至少是西方國家的那些類英語的語言）寫成的文檔，不論主題或長短，通常得出來的 Profile 都差不多，亦即按照出現的頻率排序所得到的各個 N-gram 的序號不會變化太大。這是非常好的一個性質：通常我們只要各個語言選取一篇（比較正常的，也不需要很長）文檔構建出一個 Profile ，在拿到一篇未知文檔的時候，只要和各個 Profile 比較一下，差異最小的那個 Profile 所對應的語言就可以認定是這篇未知文檔的語言了——准確率很高，更可貴的是，所需要的訓練數據非常少而且容易獲得，訓練出來的模型也是非常小的。

   不過，我們這里且撇開分類（Classification）的問題，回到聚類（Clustering）上，按照前面的說法，在 k-medoids 聚類中，只需要定義好兩個東西之間的距離（或者 dissimilarity ）就可以了，對於兩個 Profile ，它們之間的 dissimilarity 可以很自然地定義為對應的 N-gram 的序號之差的絕對值，在 Python 中用下面這樣一個類來表示：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

   class Profile(object): def __init__(self, path, N=3, psize=400): self.N = N self.psize = psize self.build_profile(path)     sep = re.compile(r'\W+') def build_profile(self, path): grams = {} with open(path) as inf: for line in inf: for tok in self.sep.split(line): for n in range(self.N): self.feed_ngram(grams, tok, n+1) self.create_profile(grams.items())     def create_profile(self, grams): # keep only the top most psize items grams.sort(key=itemgetter(1), reverse=True) grams = grams[:self.psize]     self.profile = dict() for i in range(len(grams)): self.profile[grams[i][0]] = i     def __getitem__(self, key): idx = self.profile.get(key) if idx is None: return len(self.profile) return idx     def dissimilarity(self, o): dis = 0 for tok in self.profile.keys(): dis += abs(self[tok]-o[tok]) for tok in o.profile.keys(): dis += abs(self[tok]-o[tok]) return dis     def feed_ngram(self, grams, tok, n): if n != 0: tok = '_' + tok tok = tok + '_' * (n-1) for i in range(len(tok)-n+1): gram = tok[i:i+n] grams.setdefault(gram, 0) grams[gram] += 1

   europarl 數據集共有 11 種語言的文檔，每種語言包括大約 600 多個文檔。我為這七千多個文檔建立了 Profile 並構造出一個 7038×7038 的 dissimilarity matrix ，最后在這上面用 k-medoids 進行聚類。構造 dissimilarity matrix 的過程很慢，在我這里花了將近 10 個小時。相比之下，k-medoids 的過程在內存允許的情況下，采用向量化的方法來做實際上還是很快的，並且通常只要數次迭代就能收斂了。實際的 k-medoids 實現可以在 mltk 中找到，今后如果有時間的話，我會陸續地把一些相關的比較通用的代碼放到那里面。

   最后，然我們來看看聚類的結果，由於聚類和分類不同，只是得到一些 cluster ，而並不知道這些 cluster 應該被打上什么標簽，或者說，在這個問題里，包含了哪種語言的文檔。由於我們的目的是衡量聚類算法的 performance ，因此直接假定這一步能實現最優的對應關系，將每個 cluster 對應到一種語言上去。一種辦法是枚舉所有可能的情況並選出最優解，另外，對於這樣的問題，我們還可以用 Hungarian algorithm 來求解。

   我們這里有 11 種語言，全排列有 11! = 39916800 種情況， 對於每一種排列，我們需要遍歷一次 label list ，並數出真正的 label （語言）與聚類得出的結果相同的文檔的個數，再除以總的文檔個數，得到 accuracy 。假設每次遍歷並求出 accuracy 只需要 1 毫秒的時間的話，總共也需要 11 個小時才能得到結果。看上去好像也不是特別恐怖，不過相比起來，用 Hungarian algorithm 的話，我們可以幾乎瞬間得到結果。由於文章的篇幅已經很長了，就不在這里介紹具體的算法了，感興趣的同學可以參考 Wikipedia ，這里我直接使用了一個現有的 Python 實現。

   雖然這個實驗非常折騰，不過最后的結果其實就是一個數字：accuracy ——在我這里達到了 88.97% ，證明 k-medoids 聚類和 N-gram Profile 識別語言這兩種方法都是挺不錯的。最后，如果有感興趣的同學，代碼可以從這里下載。需要最新版的 scipy， munkres.py 和 mltk 以及 Python 2.6 。

7. (番外篇): Vector Quantization

   在接下去說其他的聚類算法之前，讓我們先插進來說一說一個有點跑題的東西：Vector Quantization 。這項技術廣泛地用在信號處理以及數據壓縮等領域。事實上，在 JPEG 和 MPEG-4 等多媒體壓縮格式里都有 VQ 這一步。

   Vector Quantization 這個名字聽起來有些玄乎，其實它本身並沒有這么高深。大家都知道，模擬信號是連續的值，而計算機只能處理離散的數字信號，在將模擬信號轉換為數字信號的時候，我們可以用區間內的某一個值去代替着一個區間，比如，[0, 1) 上的所有值變為 0 ，[1, 2) 上的所有值變成 1 ，如此類推。其這就是一個 VQ 的過程。一個比較正式一點的定義是：VQ 是將一個向量空間中的點用其中的一個有限子集來進行編碼的過程。

   一個典型的例子就是圖像的編碼。最簡單的情況，考慮一個灰度圖片，0 為黑色，1 為白色，每個像素的值為 [0, 1] 上的一個實數。現在要把它編碼為 256 階的灰階圖片，一個最簡單的做法就是將每一個像素值 x 映射為一個整數 floor(x*255) 。當然，原始的數據空間也並不以一定要是連續的。比如，你現在想要把壓縮這個圖片，每個像素只使用 4 bit （而不是原來的 8 bit）來存儲，因此，要將原來的 [0, 255] 區間上的整數值用 [0, 15] 上的整數值來進行編碼，一個簡單的映射方案是 x*15/255 。

   

   Rechard Stallman         VQ2

   

   VQ100 VQ10

   不過這樣的映射方案頗有些 Naive ，雖然能減少顏色數量起到壓縮的效果，但是如果原來的顏色並不是均勻分布的，那么的出來的圖片質量可能並不是很好。例如，如果一個 256 階灰階圖片完全由 0 和 13 兩種顏色組成，那么通過上面的映射就會得到一個全黑的圖片，因為兩個顏色全都被映射到 0 了。一個更好的做法是結合聚類來選取代表性的點。

   實際做法就是：將每個像素點當作一個數據，跑一下 K-means ，得到 k 個 centroids ，然后用這些 centroids 的像素值來代替對應的 cluster 里的所有點的像素值。對於彩色圖片來說，也可以用同樣的方法來做，例如 RGB 三色的圖片，每一個像素被當作是一個 3 維向量空間中的點。

   用本文開頭那張 Rechard Stallman 大神的照片來做一下實驗好了，VQ 2、VQ 10 和 VQ 100 三張圖片分別顯示聚類數目為 2 、10 和 100 時得到的結果，可以看到 VQ 100 已經和原圖非常接近了。把原來的許多顏色值用 centroids 代替之后，總的顏色數量減少了，重復的顏色增加了，這種冗余正是壓縮算法最喜歡的。考慮一種最簡單的壓縮辦法：單獨存儲（比如 100 個）centroids 的顏色信息，然后每個像素點存儲 centroid 的索引而不是顏色信息值，如果一個 RGB 顏色值需要 24 bits 來存放的話，每個（128 以內的）索引值只需要 7 bits 來存放，這樣就起到了壓縮的效果。

   實現代碼很簡單，直接使用了 SciPy 提供的 kmeans 和 vq 函數，圖像讀寫用了 Python Image Library ：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

   #!/usr/bin/python     from scipy.cluster.vq import kmeans, vq from numpy import array, reshape, zeros from mltk import image     vqclst = [2, 10, 100, 256]     data = image.read('example.jpg') (height, width, channel) = data.shape     data = reshape(data, (height*width, channel)) for k in vqclst: print 'Generating vq-%d...' % k (centroids, distor) = kmeans(data, k) (code, distor) = vq(data, centroids) print 'distor: %.6f' % distor.sum() im_vq = centroids[code, :] image.write('result-%d.jpg' % k, reshape(im_vq, (height, width, channel)))

   當然，Vector Quantization 並不一定要用 K-means 來做，各種能用的聚類方法都可以用，只是 K-means 通常是最簡單的，而且通常都夠用了。

8. Gaussian Mixture Model

   上一次我們談到了用 k-means 進行聚類的方法，這次我們來說一下另一個很流行的算法：Gaussian Mixture Model (GMM)。事實上，GMM 和k-means 很像，不過GMM 是學習出一些概率密度函數來（所以GMM除了用在clustering上之外，還經常被用於density estimation），簡單地說，k-means 的結果是每個數據點被 assign 到其中某一個cluster 了，而GMM 則給出這些數據點被 assign 到每個 cluster 的概率，又稱作 soft assignment 。

   得出一個概率有很多好處，因為它的信息量比簡單的一個結果要多。比如，我可以把這個概率轉換為一個score，表示算法對自己得出的這個結果的把握。也許我可以對同一個任務，用多個方法得到結果，最后選取"把握"最大的那個結果；另一個很常見的方法是在諸如疾病診斷之類的場所，機器對於那些很容易分辨的情況（患病或者不患病的概率很高）可以自動區分，而對於那種很難分辨的情況，比如，49% 的概率患病，51% 的概率正常，如果僅僅簡單地使用 50% 的閾值將患者診斷為"正常"的話，風險是非常大的，因此，在機器對自己的結果把握很小的情況下，會"拒絕發表評論"，而把這個任務留給有經驗的醫生去解決。

   廢話說了一堆，不過，在回到 GMM 之前，我們再稍微扯幾句。我們知道，不管是機器還是人，學習的過程都可以看作是一種"歸納"的過程，在歸納的時候你需要有一些假設的前提條件，例如，當你被告知水里游的那個家伙是魚之后，你使用"在同樣的地方生活的是同一種東西"這類似的假設，歸納出"在水里游的都是魚"這樣一個結論。當然這個過程是完全"本能"的，如果不仔細去想，你也不會了解自己是怎樣"認識魚"的。另一個值得注意的地方是這樣的假設並不總是完全正確的，甚至可以說總是會有這樣那樣的缺陷的，因此你有可能會把蝦、龜、甚至是潛水員當做魚。也許你覺得可以通過修改前提假設來解決這個問題，例如，基於"生活在同樣的地方並且穿着同樣衣服的是同一種東西"這個假設，你得出結論：在水里有並且身上長有鱗片的是魚。可是這樣還是有問題，因為有些沒有長鱗片的魚現在又被你排除在外了。

   在這個問題上，機器學習面臨着和人一樣的問題，在機器學習中，一個學習算法也會有一個前提假設，這里被稱作"歸納偏執 (bias)"（bias 這個英文詞在機器學習和統計里還有其他許多的意思）。例如線性回歸，目的是要找一個函數盡可能好地擬合給定的數據點，它的歸納偏執就是"滿足要求的函數必須是線性函數"。一個沒有歸納偏執的學習算法從某種意義上來說毫無用處，就像一個完全沒有歸納能力的人一樣，在第一次看到魚的時候有人告訴他那是魚，下次看到另一條魚了，他並不知道那也是魚，因為兩條魚總有一些地方不一樣的，或者就算是同一條魚，在河里不同的地方看到，或者只是看到的時間不一樣，也會被他認為是不同的，因為他無法歸納，無法提取主要矛盾、忽略次要因素，只好要求所有的條件都完全一樣──然而哲學家已經告訴過我們了：世界上不會有任何樣東西是完全一樣的，所以這個人即使是有無比強悍的記憶力，也絕學不到任何一點知識。

   這個問題在機器學習中稱作"過擬合 (Overfitting)"，例如前面的回歸的問題，如果去掉"線性函數"這個歸納偏執，因為對於 N 個點，我們總是可以構造一個 N-1 次多項式函數，讓它完美地穿過所有的這 N 個點，或者如果我用任何大於 N-1 次的多項式函數的話，我甚至可以構造出無窮多個滿足條件的函數出來。如果假定特定領域里的問題所給定的數據個數總是有個上限的話，我可以取一個足夠大的 N ，從而得到一個（或者無窮多個）"超級函數"，能夠 fit 這個領域內所有的問題。然而這個（或者這無窮多個）"超級函數"有用嗎？只要我們注意到學習的目的（通常）不是解釋現有的事物，而是從中歸納出知識，並能應用到新的事物上，結果就顯而易見了。

   沒有歸納偏執或者歸納偏執太寬泛會導致 Overfitting ，然而另一個極端──限制過大的歸納偏執也是有問題的：如果數據本身並不是線性的，強行用線性函數去做回歸通常並不能得到好結果。難點正在於在這之間尋找一個平衡點。不過人在這里相對於（現在的）機器來說有一個很大的優勢：人通常不會孤立地用某一個獨立的系統和模型去處理問題，一個人每天都會從各個來源獲取大量的信息，並且通過各種手段進行整合處理，歸納所得的所有知識最終得以統一地存儲起來，並能有機地組合起來去解決特定的問題。這里的"有機"這個詞很有意思，搞理論的人總能提出各種各樣的模型，並且這些模型都有嚴格的理論基礎保證能達到期望的目的，然而絕大多數模型都會有那么一些"參數"（例如 K-means 中的 k），通常沒有理論來說明參數取哪個值更好，而模型實際的效果卻通常和參數是否取到最優值有很大的關系，我覺得，在這里"有機"不妨看作是所有模型的參數已經自動地取到了最優值。另外，雖然進展不大，但是人們也一直都期望在計算機領域也建立起一個統一的知識系統（例如語意網就是這樣一個嘗試）。

   廢話終於說完了，回到 GMM 。按照我們前面的討論，作為一個流行的算法，GMM 肯定有它自己的一個相當體面的歸納偏執了。其實它的假設非常簡單，顧名思義，Gaussian Mixture Model ，就是假設數據服從 Mixture Gaussian Distribution ，換句話說，數據可以看作是從數個 Gaussian Distribution 中生成出來的。實際上，我們在 K-means 和 K-medoids 兩篇文章中用到的那個例子就是由三個 Gaussian 分布從隨機選取出來的。實際上，從中心極限定理可以看出，Gaussian 分布（也叫做正態 (Normal) 分布）這個假設其實是比較合理的，除此之外，Gaussian 分布在計算上也有一些很好的性質，所以，雖然我們可以用不同的分布來隨意地構造 XX Mixture Model ，但是還是 GMM 最為流行。另外，Mixture Model 本身其實也是可以變得任意復雜的，通過增加 Model 的個數，我們可以任意地逼近任何連續的概率密分布。

   每個 GMM 由  個 Gaussian 分布組成，每個 Gaussian 稱為一個"Component"，這些 Component 線性加成在一起就組成了 GMM 的概率密度函數：

   根據上面的式子，如果我們要從GMM 的分布中隨機地取一個點的話，實際上可以分為兩步：首先隨機地在這 K個 Component 之中選一個，每個 Component 被選中的概率實際上就是它的系數  ，選中了 Component 之后，再單獨地考慮從這個 Component 的分布中選取一個點就可以了──這里已經回到了普通的 Gaussian 分布，轉化為了已知的問題。

   那么如何用 GMM 來做 clustering 呢？其實很簡單，現在我們有了數據，假定它們是由 GMM 生成出來的，那么我們只要根據數據推出 GMM 的概率分布來就可以了，然后 GMM 的 K個 Component 實際上就對應了 K 個 cluster 了。根據數據來推算概率密度通常被稱作 density estimation ，特別地，當我們在已知（或假定）了概率密度函數的形式，而要估計其中的參數的過程被稱作"參數估計"。

   現在假設我們有 N 個數據點，並假設它們服從某個分布（記作  ），現在要確定里面的一些參數的值，例如，在 GMM 中，我們就需要確定 、 和  這些參數。 我們的想法是，找到這樣一組參數，它所確定的概率分布生成這些給定的數據點的概率最大，而這個概率實際上就等於，我們把這個乘積稱作似然函數 (Likelihood Function)。通常單個點的概率都很小，許多很小的數字相乘起來在計算機里很容易造成浮點數下溢，因此我們通常會對其取對數，把乘積變成加和 ，得到 log-likelihood function 。接下來我們只要將這個函數最大化（通常的做法是求導並令導數等於零，然后解方程），亦即找到這樣一組參數值，它讓似然函數取得最大值，我們就認為這是最合適的參數，這樣就完成了參數估計的過程。

   下面讓我們來看一看 GMM 的 log-likelihood function：

   由於在對數函數里面又有加和，我們沒法直接用求導解方程的辦法直接求得最大值。為了解決這個問題，我們采取之前從 GMM 中隨機選點的辦法：分成兩步，實際上也就類似於 K-means 的兩步。

   估計數據由每個 Component 生成的概率（並不是每個 Component 被選中的概率）：對於每個數據  來說，它由第 k 個 Component 生成的概率為

   由於式子里的  和  也是需要我們估計的值，我們采用迭代法，在計算  的時候我們假定  和 均已知，我們將取上一次迭代所得的值（或者初始值）。

   估計每個 Component 的參數：現在我們假設上一步中得到的  就是正確的"數據 由 Component k 生成的概率"，亦可以當做該 Component 在生成這個數據上所做的貢獻，或者說，我們可以看作  這個值其中有 這部分是由 Component k 所生成的。集中考慮所有的數據點，現在實際上可以看作 Component 生成了  這些點。由於每個 Component 都是一個標准的 Gaussian 分布，可以很容易分布求出最大似然所對應的參數值：

   其中，並且  也順理成章地可以估計為 。

   重復迭代前面兩步，直到似然函數的值收斂為止。

   當然，上面給出的只是比較"直觀"的解釋，想看嚴格的推到過程的話，可以參考 Pattern Recognition and Machine Learning 這本書的第九章。有了實際的步驟，再實現起來就很簡單了。Matlab 代碼如下：

   （Update 2012.07.03：如果你直接把下面的代碼拿去運行了，碰到 covariance 矩陣 singular 的情況，可以參見這篇文章。）

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102

   function varargout = gmm(X, K_or_centroids) % ============================================================ % Expectation-Maximization iteration implementation of % Gaussian Mixture Model. % % PX = GMM(X, K_OR_CENTROIDS) % [PX MODEL] = GMM(X, K_OR_CENTROIDS) % % - X: N-by-D data matrix. % - K_OR_CENTROIDS: either K indicating the number of % components or a K-by-D matrix indicating the % choosing of the initial K centroids. % % - PX: N-by-K matrix indicating the probability of each % component generating each point. % - MODEL: a structure containing the parameters for a GMM: % MODEL.Miu: a K-by-D matrix. % MODEL.Sigma: a D-by-D-by-K matrix. % MODEL.Pi: a 1-by-K vector. % ============================================================     threshold = 1e-15; [N, D] = size(X);     if isscalar(K_or_centroids) K = K_or_centroids; % randomly pick centroids rndp = randperm(N); centroids = X(rndp(1:K), :); else K = size(K_or_centroids, 1); centroids = K_or_centroids; end     % initial values [pMiu pPi pSigma] = init_params();     Lprev = -inf; while true Px = calc_prob();     % new value for pGamma pGamma = Px .* repmat(pPi, N, 1); pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);     % new value for parameters of each Component Nk = sum(pGamma, 1); pMiu = diag(1./Nk) * pGamma' * X; pPi = Nk/N; for kk = 1:K Xshift = X-repmat(pMiu(kk, :), N, 1); pSigma(:, :, kk) = (Xshift' * ... (diag(pGamma(:, kk)) * Xshift)) / Nk(kk); end     % check for convergence L = sum(log(Px*pPi')); if L-Lprev < threshold break; end Lprev = L; end     if nargout == 1 varargout = {Px}; else model = []; model.Miu = pMiu; model.Sigma = pSigma; model.Pi = pPi; varargout = {Px, model}; end     function [pMiu pPi pSigma] = init_params() pMiu = centroids; pPi = zeros(1, K); pSigma = zeros(D, D, K);     % hard assign x to each centroids distmat = repmat(sum(X.*X, 2), 1, K) + ... repmat(sum(pMiu.*pMiu, 2)', N, 1) - ... 2*X*pMiu'; [dummy labels] = min(distmat, [], 2);     for k=1:K Xk = X(labels == k, :); pPi(k) = size(Xk, 1)/N; pSigma(:, :, k) = cov(Xk); end end     function Px = calc_prob() Px = zeros(N, K); for k = 1:K Xshift = X-repmat(pMiu(k, :), N, 1); inv_pSigma = inv(pSigma(:, :, k)); tmp = sum((Xshift*inv_pSigma) .* Xshift, 2); coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma)); Px(:, k) = coef * exp(-0.5*tmp); end end end

   函數返回的 Px 是一個  的矩陣，對於每一個  ，我們只要取該矩陣第 i 行中最大的那個概率值所對應的那個 Component 為 所屬的 cluster 就可以實現一個完整的聚類方法了。對於最開始的那個例子，GMM 給出的結果如下：

    

   

   相對於之前 K-means 給出的結果，這里的結果更好一些，左下角的比較稀疏的那個 cluster 有一些點跑得比較遠了。當然，因為這個問題原本就是完全有 Mixture Gaussian Distribution 生成的數據，GMM （如果能求得全局最優解的話）顯然是可以對這個問題做到的最好的建模。

   另外，從上面的分析中我們可以看到 GMM 和 K-means 的迭代求解法其實非常相似（都可以追溯到 EM 算法，下一次會詳細介紹），因此也有和 K-means 同樣的問題──並不能保證總是能取到全局最優，如果運氣比較差，取到不好的初始值，就有可能得到很差的結果。對於 K-means 的情況，我們通常是重復一定次數然后取最好的結果，不過 GMM 每一次迭代的計算量比 K-means 要大許多，一個更流行的做法是先用 K-means （已經重復並取最優值了）得到一個粗略的結果，然后將其作為初值（只要將 K-means 所得的 centroids 傳入 gmm 函數即可），再用 GMM 進行細致迭代。

   如我們最開始所討論的，GMM 所得的結果（Px）不僅僅是數據點的 label ，而包含了數據點標記為每個 label 的概率，很多時候這實際上是非常有用的信息。最后，需要指出的是，GMM 本身只是一個模型，我們這里給出的迭代的辦法並不是唯一的求解方法。感興趣的同學可以自行查找相關資料。

9. Regularized Gaussian Covariance Estimation

   我之前寫過一篇介紹 Gaussian Mixture Model (GMM) 的文章，並在文章里貼了一段 GMM實現的 Matlab 示例代碼，然后就不斷地有人來問我關於那段代碼的問題，問得最多的就是大家經常發現在跑那段代碼的時候估計出來的 Covariance Matrix 是 singular 的，所以在第 96 行求逆的時候會掛掉。這是今天要介紹的主要話題，我會講得羅嗦一點，把關於那篇文章我被問到的一些其他問題也都提到一下，不過，在步入正題之前，不妨先來小小地閑扯一下。

   我自己以前就非常在痛恨書里看到偽代碼，又不能編譯運行，還搞得像模像樣的像代碼一般，所以我自己在寫 blog 的時候就嘗試給出直接可運行的代碼，但是現在我漸漸理解為什么書的作者喜歡用偽代碼而不是可運行的代碼了（除非是專門講某編程語言的書）。原因很簡單：示例用的代碼和真實項目中的代碼往往是差別很大的，直接給出可運行的示例代碼，讀者就會直接拿去運行（因為包括我在內的大部分人都是很懶的，不是說"懶"是程序員的天性么？），而往往示例代碼為了結構清晰並突出算法的主要脈絡，對很多瑣碎情況都沒有處理，都使用最直接的操作，沒有做任何優化（例如我的那個 GMM 示例代碼里直接用 matlab 的inv 函數來求 Covariance 矩陣的逆），等等，因此直接運行這樣的代碼——特別是在實際項目中使用，是非常不明智的。

   可是那又為什么不直接給出實際項目級別的代碼呢？第一個原因當然是工作量，我想程序員都知道從一個算法到一個健壯、高效的系統之間有多長的路要走，而且很多時候都需要根據不同的項目環境和需要做不同的修改。所以當一個人實際上是在介紹算法的時候讓他去弄一套工業級別的代碼出來作為示例，實在是太費力了。當然更重要的原因是：工業級別的代碼通常經過大量的優化並充斥着大量錯誤處理、以及為了處理其他地方的 bug 或者整個系統混亂的 API 等而做的各種 workaround 等等……也許根本就看不懂，基本完全失去了示例的作用。所以，偽代碼就成為了唯一的選擇。

   羅嗦了半天，現在我們回到正題，那篇文章講的是 GMM 的學習，所以關於 GMM 的問題可以直接參考那篇文章，出問題的地方僅在於 GMM 的每個 component (Gaussian Distribution) 的參數估計上，特別地，在於 Covariance Matrix 的估計上。所以，我們今天主要來探討 Gaussian 分布的協方差矩陣的估計問題。首先來看一下多維 Gaussian 分布的概率密度函數：

   eq: 1 »

   參數估計的過程是給定一堆數據點 ，並假設他們是由某個真實的 Gaussian 分布獨立地 (IID) 采樣出來的，然后根據這堆點來估計該 Gaussian 分布的參數。對於一個多維 Gaussian 分布來說，需要估計的就是均值  和協方差矩陣  兩"個"參數——用引號是因為 和  都不是標量，假設我們的數據是  維的，那么  實際上是  個參數，而作為一個 的矩陣， 則是由  個參數組成。根據協方差矩陣的對稱性，我們可以進一步把  所包含的參數個數降低為  個。

   值得一提的是，我們從參數估計的角度並不能一下子斷定  就是對稱的，需要稍加說明（這實際上是《Pattern Recognition and Machine Learning》書里的習題 2.17）。記 ，假設我們的  不是對稱的，那么我們把他寫成一個對稱矩陣和一個反對稱陣之和（任何一個矩陣都可以這樣表示）：

   將起代入式 (eq: 1) 的指數部分中，注意到對於反對稱部分（書寫方便起見我們暫時另 ）：

   所以指數部分只剩下對稱陣  了。也就是說，如果我們找出一個符合條件的非對稱陣 ，用它的對稱部分  代替也是可以的，所以我們就可以"不失一般性"地假設  是對稱的（因為對稱陣的逆矩陣也是對稱的）。注意我們不用考慮式 (eq: 1) 指數前面的系數（里的那個  的行列式），因為前面的系數是起 normalization 作用的，使得全概率積分等於 1，指數部分變動之后只要重新 normalize 一下即可（可以想像一下，如果  真的是非對稱的，那么 normalization 系數那里應該也不會是  那么簡單的形式吧？有興趣的同學可以自己研究一下）。

   好，現在回到參數估計上。對於給定的數據 ，將其代入式 (eq: 1) 中可以得到其對應的 likelihood ，注意這個並不是  的概率。和離散概率分布不同的是，這里談論單個數據點的概率是沒有多大意義的，因為單點集的測度為零，它的概率也總是為零。所以，這里的 likelihood 出現大於 1 的值也是很正常的——這也是我經常被問到的一個問題。雖然不是概率值，但是我想從 likelihood 或者概率密度的角度去理解，應該也是可以得到直觀的認識的。

   參數估計常用的方法是最大似然估計（Maximum Likelihood Estimation, MLE）。就是將所有給定數據點的似然全部乘起來，得到一個關於參數（這里是  和 ）的函數，然后最大化該函數，所對應的參數值就是最大似然估計的值。通常為了避免數值下溢，會使用單調函數  將相乘變成相加在計算。

   對於最大似然估計，我的理解是尋找最可能生成給定的數據的模型參數。這對於離散型的分布來說是比較好解釋的，因為此時一點的似然實際上就是該點的概率，不過對於像 Gaussian 這樣的針對連續數據的分布來說，就不是那么自然了。不過我們其實有另一種解釋，這是我在Coursera 的 Probabilistic Graphical Model 課上看到的：

   對於一個真實的分布 ，假設我們得到了它的一個估計 ，那么如何來衡量這個估計的好壞呢？比較兩個概率分布的一個常用的辦法是用 KL-Divergence，又叫 Relative Entropy（以及其他若干外號），定義為：

   不過由於我們通常是不知道真實的  的，所以這個值也沒法直接算，於是我們把它簡單地變形一下：

   其中藍色部分是分布  自己的 Entropy，由於我們比較不同的模型估計好壞的時候這個量是不變的，所以我們可以忽略它，只考慮紅色部分。為了使得 KL-Divergence 最小，我們需要將紅色部分（注意不包括前面的負號）最大化。不過由於仍然依賴於真實分布 ，所以紅色部分還是不能直接計算，但是這是一個我們熟悉的形式：某個函數關於真實分布的期望值，而我們手里面又有從真實分布里 IID 采樣出來的  個數據點，所以可以直接使用標准的近似方法：用數據的經驗分布期望來代替原始的期望：

   於是我們很自然地就得到了最大似然估計的目標函數。這也從另一個角度解釋了最大似然估計的合理性：同時對離散型和連續型的數據都適用。

   接下來我們再回到 Gaussian 分布，將 Gaussian 的概率密度函數 (eq: 1) 代入最大似然估計的目標函數，可以得到：

   eq: 2 »

   由於本文主要討論協方差矩陣的估計，我們這里就不管均值矩陣  的估計了，方法上是類似的，下面式子中的  將表示已知的或者同樣通過最大似然估計計算出來的均值。為了找到使得上式最大化的協方差矩陣，我們將它對於  求導：

   另導數等於零，即可得到最優解：

   eq: 3 »

   於是  就是協方差矩陣的最大似然估計。到這里問題就出來了，再回到 Gaussian 分布的概率密度函數 (eq: 1)，如果把  帶進去的話，可以看到是需要對它進行求逆的，這也是我被問得最多的問題：由於  singular 了，導致求逆失敗，於是代碼也就運行失敗了。

   從 (eq: 3) 可以大致看到什么時候  會 singular：由於該式的形式，可以知道  的 rank 最大不會超過 ，所以如果 ，也就是數據的維度大於數據點的個數的話，得到的估計值  肯定是不滿秩的，於是就 singular 了。當然，即使不是這么極端的情況，如果  比較小或者 比較大的時候，仍然會有很大的幾率出現不滿秩的。

   這個問題還可以從更抽象一點的角度來看。我們之前計算過估計 Covariance 矩陣實際上是要估計  這么多個參數，如果  比較大的話，這個數目也會變得非常多，也就意味着模型的復雜度很大，對於很復雜的模型，如果訓練數據的量不夠多，則無法確定一個唯一的最優模型出來，這是機器學習中非常常見的一個問題，通常稱為 overfitting。解決 overfitting 問題的方法有不少，最直接的就是用更多的訓練數據（也就是增大 ），當然，有時候訓練數據只有那么多，於是就只好從反面入手，降低模型復雜度。

   在降低模型復雜度方面最常用的方法大概是正則化了，而正則化最簡單的例子也許是 Ridge Regression 了，通過在目標函數后面添加一項關於參數的 -norm 的項來對模型參數進行限制；此外也有諸如 LASSO 之類的使用 -norm 來實現稀疏性的正則化，這個我們在之前也曾經簡單介紹過。

   在介紹估計 Covariance 的正則化方法之前，我們首先來看一種更直接的方法，直接限制模型的復雜度：特別地，對於我們這里 Gaussian 的例子，比如，我們可以要求協方差矩陣是一個對角陣，這樣一來，對於協方差矩陣，我們需要估計的參數個數就變成了  個對角元素，而不是原來的  個。大大減少參數個數的情況下，overfitting 的可能性也極大地降低了。

   注意 Covariance Matrix 的非對角線上的元素是對應了某兩個維度之間的相關性 (Correlation) 的，為零就表示不相關。而對於 Gaussian 分布來說，不相關和獨立 (Independence) 是等價的，所以說在我們的假設下，向量  的各個維度之間是相互獨立的，換句話說，他們的聯合概率密度可以分解為每個維度各自的概率密度（仍然是 Gaussian 分布）的乘積，這個特性使得我們可以很方便地對協方差矩陣進行估計：只有每個維度上當做一維的 Gaussian 分布各自獨立地進行參數估計就可以了。注意這個時候除非某個維度上的坐標全部是相等的，否則 Covariance 矩陣對角線上的元素都能保證大於零，也就不會出現 singular 的情況了。

   那么直觀上將協方差矩陣限制為對角陣會產生什么樣的模型呢？我們不妨看一個二維的簡單例子。下面的代碼利用了 scikit-learn 里的 GMM 模塊來估計一個單個 component 的 Gaussian，它在選項里已經提供了 Covariance 的限制，這里我們除了原始的 full 和對角陣的 diag 之外，還給出一個 spherical，它假定協方差矩陣是  這樣的形式，也就是說，不僅是對角陣，而且所有對角元素都相等。下面的代碼把三種情況下 fit 出來的模型畫出來，在圖 (fig: 1)中可以看到。

   1. import numpy as np 

   2. import pylab as pl 

   3. from sklearn import mixture 

   4.  

   5. n_samples = 300 

   6. c_types = ['full', 'diag', 'spherical'] 

   7.  

   8. np.random.seed(0) 

   9. C = np.array([[0., -0.7], [3.5, 1.7]]) 

   10. X_train = np.dot(np.random.randn(n_samples, 2), C) 

   11.  

   12. pl.figure(dpi=100,figsize=(3,3)) 

   13. pl.scatter(X_train[:, 0], X_train[:, 1], .8) 

   14. pl.axis('tight') 

   15. pl.savefig('GaussianFit-data.svg') 

   16. pl.close() 

   17.  

   18. for c_type in c_types: 

   19.  clf = mixture.GMM(n_components=1, covariance_type=c_type) 

   20.  clf.fit(X_train) 

   21.  

   22.  x = np.linspace(-15.0, 20.0, num=200) 

   23.  y = np.linspace(-10.0, 10.0, num=200) 

   24.  X, Y = np.meshgrid(x, y) 

   25.  XX = np.c_[X.ravel(), Y.ravel()] 

   26.  Z = np.log(-clf.eval(XX)[0]) 

   27.  Z = Z.reshape(X.shape) 

   28.  

   29.  pl.figure(dpi=100,figsize=(3,3)) 

   30.  CS = pl.contour(X, Y, Z) 

   31.  pl.scatter(X_train[:, 0], X_train[:, 1], .8) 

   32.  

   33.  pl.axis('tight') 

   34.  pl.savefig('GaussianFit-%s.svg' % c_type) 

   35.  pl.close() 

   fig: 1 »

•  

  Training samples.

• 

  Spherical Covariance.

• 

  Diagonal Covariance.

• 

  Full Covariance.

  可以看到，從直觀上來講：spherical Covariance 的時候就是各個維度上的方差都是一樣的，所以形成的等高線在二維上是正圓。而 diagonal Covariance 則可以 capture 更多的特征，因為它允許形成橢圓的等高線，但是它的限制是這個橢圓是不能旋轉的，永遠是正放的。Full Covariance 是 fit 最好的，但是前提是數據量足夠來估計這么多的參數，雖然對於 2 維這樣的低維度通常並不是問題，但是維度高了之后就經常會變得很困難。

  一般來說，使用對角線的協方差矩陣是一種不錯的降低模型復雜度的方式，在 scikit-learn 的 GMM 模塊中這甚至是默認的。不過如果有時候你覺得這樣的限制實在是和你的真實數據想去甚遠的話，也有另外的處理方式。最簡單的處理辦法是在求出了 Covariance 矩陣的估計值之后直接加上一個 ，也就是在對角線上統一加上一個很小的正數 ，通常稱作正則化系數。例如，Matlab 自帶的統計工具箱里的 gmdistribution.fit 函數就有一個可選參數叫做Regularize，就是我們這里說的 。為什么加上  之后就不會 singular 了呢？首先  本身至少肯定是半正定的，所以 ：

  所以這樣之后肯定就是正定的了。此外在 scikit-learn 的 GMM 學習的函數中也有一個類似的參數 min_covar，它的文檔里說的是"Floor on the diagonal of the covariance matrix to prevent overfitting. Defaults to 1e-3."感覺好像是如果原來的 Covariance 矩陣對角線元素小於該參數的時候就將它設置為該參數，不過這樣的做法是不是一定產生正定的協方差矩陣似乎就不像上面那種那樣直觀可以得出結論了。

  不過這樣的做法雖然能夠解決 singular 的問題，但是總讓人有些難以信服：你莫名其妙地在我的協方差矩陣上加上了一些東西，這到底是什么意思呢？最簡單的解釋可以從式 (eq: 1) 那里協方差矩陣你的地方來看，對於矩陣求逆或者解方程的時候出現 singular 的情況，加上一個 也算是數值上的一種標准處理方式，叫做 Tikhonov Regularization，Ridge Regression 本身也可以通過這個角度來解釋。從數值計算的角度來說，這樣處理能夠增加數值穩定性，直觀地來講，穩定性是指  元素的微小數值變化，不會使得求逆（或者解方程之類的）之后的解產生巨大的數值變化，如果  是 singular 的話，通常就不具有這樣的穩定性了。此外，隨着 regularization coefficient  逐漸趨向於零，對應的解也會逐漸趨向於真實的解。

  如果這個解釋還不能令出信服的話，我們還可以從 Bayes 推斷的角度來看這個問題。不過這本身是一個很大的話題，我在這里只能簡略地講一個思路，想了解更多的同學可以參考《Pattern Recognition and Machine Learning》第二章或者其他相關的書籍。

  總的來說，我們這里要通過 Maximum a posteriori (MAP) Estimation 來對協方差矩陣進行估計。做法是首先為  定義一個先驗分布 (prior) ，然后對於數據 ，我們根據 Bayes 公式，有

  eq: 4 »

  其中等式左邊稱作后驗 (posterior)，右邊紅色部分是似然 (likelihood)，藍色部分是用於將概率分布 normalize 到全概率為 1 的，通常並不需要直接通過  去計算，而只要對求出的概率分布進行歸一化就可以了。MAP 的思路是將 posterior 最大的那個  作為最佳估計值，直觀上很好理解，就是給定了數據點之后"最可能的"那個 。而計算上我們是通過 (eq: 4) 來實現的，因為分母是與  無關的，所以在計算的時候（因為只要比較相對大小嘛）忽略掉，這樣一來，就可以看到，其實 MAP 方法和剛才我們用的最大似然 (MLE) 方法唯一的區別就在於前面乘了一個先驗分布。

  雖然從最終的計算公式上來看差別很細微，但是 MAP 卻有很多好處。和本文最相關的當然是先驗的引入，這在很大程度上和我們平時用正則化的目的是一樣的：加入我們自己對模型的先驗知識或者假設或者要求，這樣可以避免在數據不夠的時候產生 overfitting 。此外還有另一個值得一提的好玩特性就是 MAP 可以在線計算：注意到后驗分布本身也是關於  的一個合法分布，所以在我們計算出后驗之后，如果又得到了新的訓練數據，那么我們可以將之前算出來的后驗又作為先驗，代入新的一輪的計算，這樣可以一直不斷地重復下去。:D

  不過，雖然 MAP 框架看上去很美，但是在 general 的情況下計算卻可能會變得很復雜，因此，雖然理論上來說，我們可以把任何先驗知識通過先驗分布的方式加入到模型中來，但是在實際中先驗的選擇往往是根據"哪個更好計算"而不是根據"哪個更合理"的准則來做出的。基於這個准則，有了 Conjugate Prior 的概念，簡單地來說，如果一個 prior 和 likelihood 相乘再歸一化之后得到的 posterior 形式上是和 prior 一樣的話，那么它就是該 likelihood 的一個 conjugate prior。

  選擇 conjugate prior 的好處是顯而易見的，由於先驗和后驗的形式是一樣的，都是某一類型的分布，只是參數的值發生了變化，所以說可以看成是在觀察到數據之后對模型參數進行修正（注意這里的模型是關於  的，而  本身又是關於數據的模型——一個 Gaussian 分布——的參數，不要搞混淆了），並且這種修正有時候可以得到非常直觀的解釋，不過我在這里就不細說了。此外，由於形式沒有發生變化，所以在使用在線計算的時候，每一次觀察到新的數據所使用的計算公式都還是一樣的。

  回到我們的問題，關於 Gaussian 分布的協方差矩陣的 Conjugate Prior 是一個叫做 Inverse Wishart Distribution 的家伙，它的概率密度函數是這個樣子的：

  其中  和參數  都是  的正定矩陣，而  則是 Multivariate Gamma Function。是不是看上去有點恐怖？^_^bbbb我也覺得有點恐怖……接下來讓我們來看 MAP，首先將這個 prior 乘以我們的 likelihood，和 MLE 一樣的，我們在  空間中計算，所以對整個式子取對數，然后丟掉和  無關的常數項（因為是關於  最大化嘛），得到下面的式子：

  可以看到現在已經變得不是那么復雜了，形式上和我們之前的 (eq: 2) 是一樣的，只是多了紅色的一項，接下來對  進行求導並令導數等於零，類似地得到如下方程：

  於是得到最優解為

  這里  是之前的最大似然估計。接下來如果我們將我們 prior 中的參數  選為 ，就可以得到 MAP 估計為：

  就得到了我們剛才的正則化的形式。有一點細微的區別就是這里前面多了一個全局的縮放系數，不過這樣的系數是沒有關系的，可以直接丟掉，因為它總也會在最終歸一化的時候被 Gaussian 分布前面的系數給抵消掉。這樣一來，我們就得到了這種正則化方法的一個看起來比較靠譜的解釋：實際上就是我們先驗地假設協方差矩陣滿足一個均值為 spherical 形狀（所有對角元素都相等的一個對角陣）的 Inverse Wishart 分布（關於 Inverse Wishart 分布的均值的公式可以參見 wikipedia），然后通過數據來對這個先驗進行修正得到后驗概率，並從后驗中取了概率（密度）最大的那個  作為最終的估計。

  這種方式比暴力地直接強制要求協方差矩陣具有 diagonal 或者 spherical 形式要溫柔得多，並且隨着訓練數據的增多，先驗的影響也會逐漸減小，從而趨向於得到真實的參數，而強制結構限制的結果則是不論你最后通過什么手段搞到了多少數據，如果真實模型本身不是 diagonal 或者 spherical 的話，那么你的參數估計的 bias 將永遠都無法被消除，這個我們已經在圖 (fig: 1)中直觀地看到了。當然，在訓練數據本身就不足的情況下，兩種選擇誰好誰好就沒有定論了，不要忘了，強制結構限制的計算復雜度要小很多。:)

  最后，對於不想看長篇大論的人，我做一下簡單的總結：如果碰到估計 Gaussian 的協方差矩陣時結果 singular 的情況，解決方法一般有兩種：

1. 直接限制協方差矩陣的結構，比如要求它是一個對角陣，這種情況下數據的各個維度將會是獨立的，所以可以把每個維度看成一個一維的 Gaussian 分布來單獨進行估計。
2. 在估計出來的協方差矩陣的對角線上統一加上一個很小的正數，使得它變成正定的。這種方法的合理性可以通過正則化或者 MAP 估計來進行解釋。

然后再補充幾點不知道該放在哪里的注意事項：

1. prior、posterior 和 likelihood 比較容易搞混淆，記住 prior 和 posterior 都是關於參數（比如我們這里就是   ）的分布，prior 是我們所假設的參數本來的分布，而 posterior 則是在觀察到訓練數據之后得到的條件分布。而 likelihood 則完全不一樣，一方面它是關於數據   的，另一方面它沒有被歸一化，所以也並不是一個合法的概率分布。
2. 標量值函數關於矩陣變量的求導原理上其實就是和多元標量值函數求導一樣的，不過求起來非常繁瑣，一般不會自己算，網上可以找到一個叫做《Matrix Cookbook》的小冊子，里面有搜集了各種常用的矩陣求導的公式，基本上直接從那里查詢就可以解決大部分問題了。

1. (番外篇): Expectation Maximization

   Expectation Maximization (EM) 是一種以迭代的方式來解決一類特殊最大似然 (Maximum Likelihood) 問題的方法，這類問題通常是無法直接求得最優解，但是如果引入隱含變量，在已知隱含變量的值的情況下，就可以轉化為簡單的情況，直接求得最大似然解。

   我們會看到，上一次說到的 Gaussian Mixture Model 的迭代求解方法可以算是 EM 算法最典型的應用，而最開始說的 K-means 其實也可以看作是 Gaussian Mixture Model 的一個變種（固定所有的  ，並令  即可）。然而 EM 實際上是一種通用的算法（或者說是框架），可以用來解決很多類似的問題，我們最后將以一個中文分詞的例子來說明這一點。

   為了避免問題變得太抽象，我們還是先從上一次的 Gaussian Mixture Model 說起吧。回顧一下我們之前要解決的問題：求以下 Log-likelihood function 的最大值：

   

   但是由於在  函數里面又有加和，沒法直接求。考慮一下 GMM 生成 sample 的過程：先選一個 Component ，然后再從這個 Component 所對應的那個普通的 Gaussian Distribution 里進行 sample 。我們可以很自然地引入一個隱含變量  ，這是一個  維向量，如果第  個 Component 被選中了，那么我們講其第  個元素置為 1 ，其余的全為 0 。那么，再來看看，如果除了之前的 sample  的值之外，我們同時還知道了每個  所對應的隱含變量  的值，情況又會變成怎么樣呢？

   因為我們同時觀察到了  和  ，所以我們現在要最大化的似然函數也變為  。注意到  可以表示為：

   

   而  的概率，當  的第  個元素為 1 的時候，亦即第  個 Component 被選中的時候，這個概率為  ，統一地寫出來就是：

   

   帶入上面個式子，我們得到  的概率是一個大的乘積式（對比之前  是一個和式）。再替換到最開始的那個 Log-likelihood function 中，得到新的同時關於 sample  和隱含變量  的 Log-likelihood:

   

   情況瞬間逆轉，現在  和求和符號換了個位置，直接作用在普通的高斯分布上了，一下子就變成了可以直接求解的問題。不過，事情之所以能發展得如此順利，完全依賴於一個我們偽造的假設：隱含變量的值已知。然而實際上我們並不知道這個值。問題的結症在這里了，如果我們有了這個值，所有的問題都迎刃而解了。回想一下，在類似的地方，我們是如何處理這樣的情況的呢？一個很類似的地方就是（比如，在數據挖掘中）處理缺失數據的情況，一般有幾種辦法：

   1. 用取值范圍類的隨機值代替。
   2. 用平均值代替。
   3. 填 0 或者其他特殊值。

   這里我們采取一種類似於平均值的辦法：取期望。因為這里我們至少有 sample  的值，因此我們可以把這個信息利用起來，針對后驗概率  來取期望。前面說過， 的每一個元素只有 0 和 1 兩種取值，因此按照期望的公式寫出來就是：

   

   中間用貝葉斯定理變了一下形，最后得到的式子正是我們在漫談 GMM 中定義的  。因此，對於上面那個可以直接求解的 Log-likelihood function 中的  ，我們只要用其期望 亦即  代替即可。

   到這里為止，看起來是一個非常完美的方法，不過仔細一看，發現還有一個漏洞：之前的 Log-likelihood function 可以直接求最大值是建立在  是已知情況下，現在雖然我們用  來代替了 ，但是實際上  是一個反過來以非常復雜的關系依賴所要求參數的一個式子，而不是一個"已知的數值"。解決這個問題的辦法就是迭代。

   到此為止，我們就勾勒出了整個 EM 算法的結構，同時也把上一次所講的求解 GMM 的方法又推導了一遍。EM 名字也來源於此：

   1. Expectation   一步對應於求關於后驗概率的期望亦即   ；
   2. Maximization   一步則對應於接下來的正常的最大似然的方法估計相關的參數亦即   、和    。

   如果你還沒有厭煩這些公式的話，讓我們不妨再稍微花一點時間，從偏理論一點的角度來簡略地證明一下 EM 這個迭代的過程每一步都會對結果有所改進，除非已經達到了一個（至少是局部的）最優解。

   現在我們的討論將不局限於 GMM ，並使用一些稍微緊湊一點的符號。用  表示所有的 sample ，同時用  表示所有對應的隱含變量。我們的問題是要通過最大似然的方法估計出  中的參數  。在這里我們假設這個問題很困難，但是要優化  卻很容易。這就是 EM 算法能夠解決的那一類問題。

   現在我們引入一個關於隱含變量的分布 。注意到對於任意的  ，我們都可以對 Log-likelihood Function 作如下分解：

   

   其中  是分布  和  之間的 Kullback-Leibler divergence 。由於 Kullback-Leibler divergence 是非負的，並且只有當兩個分布完全相同的時候才會取到 0 。因此我們可以得到關系  ，亦即  是  的一個下界。

   現在考慮 EM 的迭代過程，記上一次迭代得出的參數為 ，現在我們要選取  以令  最大，由於  並不依賴於  ，因此  的上限（在  固定的時候）是一個定值，它取到這個最大值的條件就是 Kullback-Leibler divergence 為零，亦即  等於后驗概率  。把它帶入到  的表達式中可以得到

   

   其中 const 是常量，而  則正是我們之前所得到的"同時包含了 sample 和隱含變量的 Log-likelihood function 關於后驗概率的期望"，因此這個對應到 EM 中的"E"一步。

   在接下來的"M"一步中，我們講固定住分布  ，再選取合適的  以將  最大化，這次其上界  也依賴於變量  ，並會隨着  的增大而增大（因為我們有前面的不等式成立）。一般情況下  增大的量會比  要多一些，這時候 Kullback-Leibler divergence 在新的參數  下又不為零了，因此我們可以進入下一輪迭代，重新回到"E"一步去求新的  ；另一方面，如果這里 Kullback-Leibler divergence 在新的參數下還是等於 0 ，那么說明我們已經達到了一個（至少是局部的）最優解，迭代過程可以結束了。

   上面的推導中我們可以看到每一次迭代 E 和 M 兩個步驟都是在對解進行改進，因此迭代的過程中得到的 likelihood 會逐漸逼近（至少是局部的）最優值。另外，在 M 一步中除了用最大似然之外，也可以引入先驗使用 Maximum a Posteriori (MAP) 的方法來做。還有一些很困難的問題，甚至在迭代的過程中 M 一步也不能直接求出最大值，這里我們通過把要求放寬──並不一定要求出最大值，只要能夠得到比舊的值更好的結果即可，這樣的做法通常稱作 Generalized EM (GEM)。

   當然，一開始我們就說了，EM 是一個通用的算法，並不只是用來求解 GMM 。下面我們就來舉一個用 EM 來做中文分詞的例子。這個例子來源於論文 Self-supervised Chinese Word Segmentation 。我在上次 MSTC 搜索引擎系列小課堂講文本處理的時候提到過。這里為了把注意力集中到 EM 上，我略去一些細節的東西，簡單地描述一下基本的模型。

   現在我們有一個字符序列 ，並希望得到一個模型 （依賴於參數 ）能用來將其詞的序列  。按照生成模型的方式來考慮，將  看成是由  生成的序列的話，我們自然希望找到那個生成  的可能性最大的  ，亦即通過最大似然的方式來估計參數  。

   然而我們不知道似然函數  應該如何去優化，因此我們引入 latent variable  ，如果我們知道  的話，似然值很容易求得：

   

   其中  的值是從模型  中已知的。但是現在我們不知道  的值，因此我們轉而取其關於后驗概率的期望：

   

   然后將這個期望針對  最大化即完成了 EM 的一次迭代。具體的做法通常是先把一個初始文本（比如所有的  的集合）按照 N-gram 分割（N-gram 在 講 K-medoids 的那篇文章中介紹過）為 ，形成一個最初的辭典，而模型  的參數  實際上就描述了各個 N-gram 的概率  ，初始值可以直接取為頻率值。有了辭典之后對於任意的  ，我們可以根據辭典枚舉出所有可能的分割  ，而每個分割的后驗概率  就是其中的單詞的概率乘積。其他的就是標准的 EM 算法的迭代步驟了。

   實際上在實際產品中我們會使用一些帶了更多啟發式規則的分詞方法（比如 MMSEG），而這個方法更大的用處在於從一堆文本中"學習"出一個詞典來（也就是  ），並且這個過程是全自動的，主要有兩個優點：

   1. 不需要人參與手工分詞、標記等。
   2. 能自動從文本中學習，換句話說，你給它一些某個領域的專業文章，它能從中學習到這個領域的專業詞匯來。

   不管怎樣，以上兩點看起來都非常誘人。不過理論離實際往往還是有不少差距的。我不知道實際分詞系統中有沒有在用這樣的方法來做訓練的。之前我曾經用 Wikipedia （繁體中文）上抓下來的一些數據跑過一下小規模的試驗，結果還可以，但是也並不如想像中的完美。因為當時也並沒有弄明白 EM 是怎么一回事，加上這個過程本身計算負荷還是非常大的，也就沒有深入下去。也許以后有時間可以再嘗試一下。

   總之，EM 是一種應用廣泛的算法，雖然講了點理論也講了例子，但是一沒有貼代碼二沒有貼圖片，似乎實在是不像我的作風，不過精力和篇幅有限，也只能到此為止了。

2. Spectral Clustering

   如果說 K-means 和 GMM 這些聚類的方法是古代流行的算法的話，那么這次要講的 Spectral Clustering 就可以算是現代流行的算法了，中文通常稱為"譜聚類"。由於使用的矩陣的細微差別，譜聚類實際上可以說是一"類"算法。

   Spectral Clustering 和傳統的聚類方法（例如 K-means）比起來有不少優點：

   1. 和  K-medoids  類似，Spectral Clustering 只需要數據之間的相似度矩陣就可以了，而不必像 K-means 那樣要求數據必須是 N 維歐氏空間中的向量。
   2. 由於抓住了主要矛盾，忽略了次要的東西，因此比傳統的聚類算法更加健壯一些，對於不規則的誤差數據不是那么敏感，而且 performance 也要好一些。許多實驗都證明了這一點。事實上，在各種現代聚類算法的比較中，K-means 通常都是作為 baseline 而存在的。  
   3. 計算復雜度比 K-means 要小，特別是在像文本數據或者平凡的圖像數據這樣維度非常高的數據上運行的時候。

   突然冒出這么一個要求比 K-means 要少，結果比 K-means 要好，算得還比 K-means 快的東西，實在是讓人不得不懷疑是不是江湖騙子啊。所以，是騾子是馬，先拉出來溜溜再說。不過，在K-medoids 那篇文章中曾經實際跑過 K-medoids 算法，最后的結果也就是一個 accuracy ，一個數字又不能畫成圖表之類的，看起來實在是沒意思，而且 K-means 跑起來實在是太慢了，所以這里我還是稍微偷懶一下，直接引用一下一篇論文里的結果吧。

   結果來自論文 Document clustering using locality preserving indexing 這篇論文，這篇論文實際上是提的另一種聚類方法（下次如果有機會也會講到），不過在它的實驗中也有 K-means 和 Spectral Clustering 這兩組數據，抽取出來如下所示：

   k	TDT2		Reuters-21578
   	K-means	SC	K-means	SC
   2	0.989	0.998	0.871	0.923
   3	0.974	0.996	0.775	0.816
   4	0.959	0.996	0.732	0.793
   …
   9	0.852	0.984	0.553	0.625
   10	0.835	0.979	0.545	0.615

   其中 TDT2 和 Reuters-21578 分別是兩個被廣泛使用的標准文本數據集，雖然在不同的數據集上得出來的結果並不能直接拿來比較，但是在同一數據集上 K-means 和 SC (Spectral Clustering) 的結果對比就一目了然了。實驗中分別抽取了這兩個數據集中若干個類別（從 2 類到 10 類）的數據進行聚類，得出的 accuracy 分別在上表中列出（我偷懶沒有全部列出來）。可以看到，Spectral Clustering 這里完勝 K-means 。

   

   這么強大的算法，再戴上"譜聚類"這么一個高深莫測的名號，若不是模型無比復雜、包羅宇宙，就肯定是某鎮山之寶、不傳秘籍吧？其實不是這樣，Spectral Clustering 不管從模型上還是從實現上都並不復雜，只需要能求矩陣的特征值和特征向量即可──而這是一個非常基本的運算，任何一個號稱提供線性代數運算支持的庫都理應有這樣的功能。而關於 Spectral Clustering 的秘籍更是滿街都是，隨便從地攤上找來一本，翻開便可以看到 Spectral Clustering 算法的全貌：

   1. 根據數據構造一個 Graph ，Graph 的每一個節點對應一個數據點，將相似的點連接起來，並且邊的權重用於表示數據之間的相似度。把這個 Graph 用鄰接矩陣的形式表示出來，記為    。一個最偷懶的辦法就是：直接用我們前面在 K-medoids 中用的相似度矩陣作為    。
   2. 把    的每一列元素加起來得到    個數，把它們放在對角線上（其他地方都是零），組成一個    的矩陣，記為    。並令    。
   3. 求出    的前    個特征值（在本文中，除非特殊說明，否則"前    個"指按照特征值的大小從小到大的順序）   以及對應的特征向量    。
   4. 把這    個特征（列）向量排列在一起組成一個    的矩陣，將其中每一行看作    維空間中的一個向量，並使用 K-means 算法進行聚類。聚類的結果中每一行所屬的類別就是原來 Graph 中的節點亦即最初的    個數據點分別所屬的類別。

   就是這么幾步，把數據做了一些詭異的變換，然后還在背后偷偷地調用了 K-means 。到此為止，你已經可以拿着它上街去招搖撞騙了。不過，如果你還是覺得不太靠譜的話，不妨再接着往下看，我們再來聊一聊 Spectral Clustering 那幾個"詭異變換"背后的道理何在。

   其實，如果你熟悉 Dimensional Reduction （降維）的話，大概已經看出來了，Spectral Clustering 其實就是通過 Laplacian Eigenmap 的降維方式降維之后再做 K-means 的一個過程──聽起來土多了。不過，為什么要剛好降到  維呢？其實整個模型還可以從另一個角度導出來，所以，讓我們不妨先跑一下題。

   在 Image Processing （我好像之前有聽說我對這個領域深惡痛絕？）里有一個問題就是對圖像進行 Segmentation （區域分割），也就是讓相似的像素組成一個區域，比如，我們一般希望一張照片里面的人（前景）和背景被分割到不同的區域中。在 Image Processing 領域里已經有許多自動或半自動的算法來解決這個問題，並且有不少方法和 Clustering 有密切聯系。比如我們在談 Vector Quantization 的時候就曾經用 K-means 來把顏色相似的像素聚類到一起，不過那還不是真正的 Segmentation ，因為如果僅僅是考慮顏色相似的話，圖片上位置離得很遠的像素也有可能被聚到同一類中，我們通常並不會把這樣一些"游離"的像素構成的東西稱為一個"區域"，但這個問題其實也很好解決：只要在聚類用的 feature 中加入位置信息（例如，原來是使用 R、G、B 三個值來表示一個像素，現在加入 x、y 兩個新的值）即可。

   另一方面，還有一個經常被研究的問題就是 Graph Cut ，簡單地說就是把一個 Graph 的一些邊切斷，讓他被打散成一些獨立聯通的 sub-Graph ，而這些被切斷的邊的權值的總和就被稱為 Cut值。如果用一張圖片中的所有像素來組成一個 Graph ，並把（比如，顏色和位置上）相似的節點連接起來，邊上的權值表示相似程度，那么把圖片分割為幾個區域的問題實際上等價於把 Graph 分割為幾個 sub-Graph 的問題，並且我們可以要求分割所得的 Cut 值最小，亦即：那些被切斷的邊的權值之和最小，直觀上我們可以知道，權重比較大的邊沒有被切斷，表示比較相似的點被保留在了同一個 sub-Graph 中，而彼此之間聯系不大的點則被分割開來。我們可以認為這樣一種分割方式是比較好的。

   實際上，拋開圖像分割的問題不談，在 Graph Cut 相關的一系列問題中，Minimum cut （最小割）本身就是一個被廣泛研究的問題，並且有成熟的算法來求解。只是單純的最小割在這里通常並不是特別適用，很多時候只是簡單地把和其他像素聯系最弱的那一個像素給分割出去了，相反，我們通常更希望分割出來的區域（的大小）要相對均勻一些，而不是一些很大的區塊和一些幾乎是孤立的點。為此，又有許多替代的算法提出來，如 Ratio Cut 、Normalized Cut 等。

   不過，在繼續討論之前，我們還是先來定義一下符號，因為僅憑文字還是很難表述清楚。將 Graph 表示為鄰接矩陣的形式，記為 ，其中  是節點  到節點  的權值，如果兩個節點不是相連的，權值為零。設  和  為 Graph 的兩個子集（沒有交集），那么兩者之間的 cut 可以正式定義為：

   

   先考慮最簡單的情況，如果把一個 Graph 分割為兩個部分的話，那么 Minimum cut 就是要最小化  （其中  表示  的補集）。但是由於這樣經常會出現孤立節點被分割出來的情況，因此又出現了 RatioCut ：

   

   以及 NormalizedCut ：

   

   其中  表示  中的節點數目，而  。兩者都可以算作  的"大小"的一種度量，通過在分母上放置這樣的項，就可以有效地防止孤立點的情況出現，達到相對平均一些的分割。事實上，Jianbo Shi 的這篇 PAMI paper：Normalized Cuts and Image Segmentation 正是把 NormalizedCut 用在圖像分割上了。

   搬出 RatioCut 和 NormalizedCut 是因為它們和這里的 Spectral Clustering 實際上有非常緊密的聯系。看看 RatioCut ，式子雖然簡單，但是要最小化它卻是一個 NP 難問題，不方便求解，為了找到解決辦法，讓我們先來做做變形。

   令  表示 Graph 的所有節點的集合，首先定義一個  維向量 ：

   

   再回憶一下我們最開始定義的矩陣  ，其實它有一個名字，叫做 Graph Laplacian ，不過，我們后面可以看到，其實有好幾個類似的矩陣都叫做這個名字：

   Usually, every author just calls "his" matrix the graph Laplacian.

   其實也可以理解，就好象現在所有的廠家都說自己的技術是"雲計算"一樣。這個  有一個性質就是：

   

   這個是對任意向量  都成立的，很好證明，只要按照定義展開就可以得到了。把我們剛才定義的那個  帶進去，就可以得到

   

   另外，如果令  為各個元素全為 1 的向量的話，直接展開可以很容易得到  和  。由於  是一個常量，因此最小化 RatioCut 就等價於最小化  ，當然，要記得加上附加條件  以及  。

   問題轉化到這個樣子就好求了，因為有一個叫做 Rayleigh quotient 的東西：

   

   他的最大值和最小值分別等於矩陣  的最大的那個特征值和最小的那個特征值，並且極值在 等於對應的特征向量時取到。由於  是常數，因此最小化  實際上也就等價於最小化  ，不過由於  的最小的特征值為零，並且對應的特征向量正好為  （我們這里僅考慮 Graph 是聯通的情況），不滿足  的條件，因此我們取第二個小的特征值，以及對應的特征向量  。

   到這一步，我們看起來好像是很容易地解決了前面那個 NP 難問題，實際上是我們耍了一個把戲：之前的問題之所以 NP 難是因為向量  的元素只能取兩個值  和  中的一個，是一個離散的問題，而我們求的的特征向量  其中的元素可以是任意實數，就是說我們將原來的問題限制放寬了。那如何得到原來的解呢？一個最簡單的辦法就是看  的每個元素是大於零還是小於零，將他們分別對應到離散情況的  和  ，不過我們也可以采取稍微復雜一點的辦法，用  的 K-means 來將  的元素聚為兩類。

   到此為止，已經有 Spectral Clustering 的影子了：求特征值，再對特征向量進行 K-means 聚類。實際上，從兩類的問題推廣到 k 類的問題（數學推導我就不再詳細寫了），我們就得到了和之前的 Spectral Clustering 一模一樣的步驟：求特征值並取前 k 個最小的，將對應的特征向量排列起來，再按行進行 K-means 聚類。分毫不差！

   用類似的辦法，NormalizedCut 也可以等價到 Spectral Clustering 不過這次我就不再講那么多了，感興趣的話（還包括其他一些形式的 Graph Laplacian 以及 Spectral Clustering 和 Random walk 的關系），可以去看這篇 Tutorial ：A Tutorial on Spectral Clustering 。

   為了緩和一下氣氛，我決定貼一下 Spectral Clustering 的一個簡單的 Matlab 實現：

   function idx = spectral_clustering(W, k) 

    D = diag(sum(W)); 

    L = D-W; 

      

    opt = struct('issym', true, 'isreal', true); 

    [V dummy] = eigs(L, D, k, 'SM', opt); 

      

    idx = kmeans(V, k); 

   end

   最后，我們再來看一下本文一開始說的 Spectral Clustering 的幾個優點：

   1. 只需要數據的相似度矩陣就可以了。這個是顯然的，因為 Spectral Clustering 所需要的所有信息都包含在    中。不過一般    並不總是等於最初的相似度矩陣——回憶一下，是我們構造出來的 Graph 的鄰接矩陣表示，通常我們在構造 Graph 的時候為了方便進行聚類，更加強到"局部"的連通性，亦即主要考慮把相似的點連接在一起，比如，我們設置一個閾值，如果兩個點的相似度小於這個閾值，就把他們看作是不連接的。另一種構造 Graph 的方法是將 n 個與節點最相似的點與其連接起來。
   2. 抓住了主要矛盾，忽略了次要的東西，Performance 比傳統的 K-means 要好。實際上 Spectral Clustering 是在用特征向量的元素來表示原來的數據，並在這種"更好的表示形式"上進行 K-means 。實際上這種"更好的表示形式"是用 Laplacian Eigenmap 進行降維的后的結果，如果有機會，下次談 Dimensionality Reduction 的時候會詳細講到。而降維的目的正是"抓住主要矛盾，忽略次要的東西"。
   3. 計算復雜度比 K-means 要小。這個在高維數據上表現尤為明顯。例如文本數據，通常排列起來是維度非常高（比如，幾千或者幾萬）的稀疏矩陣，對稀疏矩陣求特征值和特征向量有很高效的辦法，得到的結果是一些 k 維的向量（通常 k 不會很大），在這些低維的數據上做 K-means 運算量非常小。但是對於原始數據直接做 K-means 的話，雖然最初的數據是稀疏矩陣，但是 K-means 中有一個求 Centroid 的運算，就是求一個平均值：許多稀疏的向量的平均值求出來並不一定還是稀疏向量，事實上，在文本數據里，很多情況下求出來的 Centroid 向量是非常稠密，這時再計算向量之間的距離的時候，運算量就變得非常大，直接導致普通的 K-means 巨慢無比，而 Spectral Clustering 等工序更多的算法則迅速得多的結果。

   說了這么多，好像有些亂，不過也只能到此打住了。最后再多嘴一句，Spectral Clustering 名字來源於 Spectral theory ，也就是用特征分解來分析問題的理論了。

   UPDATE 2011.11.23: 有不少同學問我關於代碼的問題，這里更新兩點主要的問題：

   1. 關於 LE 降維的維度和 Kmeans 聚類的類別數的關系：我上面的代碼里，取成了一樣的，但是這兩者並不是要求一定要一樣的。最初 Spectral Cluster 是分析聚兩類的情況，就降到 1 維，然后用 thresholding 的方法來分成兩類。對於K 類的情況，自然的類比就是降到 K-1 維，這也是和  LDA  保持一致。因為 Laplacian 矩陣的特征向量有一個全一的向量（對應於特征值 0 ），所以可以求 K 個特征向量然后把特征值 0 對應的那個特征向量去掉。但是在實際中並不一定要保持這兩者一致的，也就是說，這個降維的維度可以作為一個參數進行調節的，選擇效果好的參數。
   2. 關於示例代碼：注意除非我在這里注明了是發布某個 software 或者 package 什么的，否則這里貼的代碼主要都是為了示例作用，為了只顯示出算法的主要部分，因此通常會省略掉一些實現細節，可以看成是"可執行的偽代碼"，不推薦直接用這里的代碼去跑實驗之類的（包括其他 post 里的相關代碼）。除非你想自己試驗一下具體實現和改進，否則可以直接找網上現成的專用的 package ，比如 Spectral Clustering 可以考慮這個包。
3. (番外篇): Dimensionality Reduction

   由於總是有各種各樣的雜事，這個系列的文章竟然一下子拖了好幾個月，（實際上其他的日志我也寫得比較少），現在決定還是先把這篇降維的日志寫完。我甚至都以及忘記了在這個系列中之前有沒有講過"特征"（feature）的概念了，這里不妨再稍微提一下。機器學習應用到各個領域里，會遇到許多不同類型的數據要處理：圖像、文本、音頻視頻以及物理、生物、化學等實驗還有其他工業、商業以及軍事上得到的各種數據，如果要為每一種類型的數據都設計獨立的算法，那顯然是非常不現實的事，因此，機器學習算法通常會采用一些標准的數據格式，最常見的一種格式就是每一個數據對應歐幾里德空間里的一個向量。

   如果原始的數據格式不兼容，那么就需要首先進行轉換，這個過程通常叫做"特征提取"（Feature Extraction），而得到的標准數據格式通常叫做 Feature 。例如，一個最簡單的將一個文本 Document 轉化為向量的方法如下：

   1. 選定特征空間，這里采用三維歐氏空間，三個維度（依次）分別由 to 、be 和 the 表示。
   2. 假設待提取的文檔是"To be, or not to be: that is the question:"，首先對其進行一些預處理，例如去掉單詞的時態后綴、濾掉標點符號等，得到"to be or not to be that be the question"。
   3. 統計三個維度所對應的單詞出現的頻率：to 2 次，be 3 次，the 1 次。
   4. 該文檔對應的向量即  [2, 3, 1]  。

   當然，在實際中我們幾乎不會這樣人工設定空間的各個維度所對應的單詞，而通常是從一個數據集中統計出所有出現的詞，再將其中的一些挑選出來作為維度。怎樣挑選呢？最簡單的辦法是根本不做任何挑選，或者簡單地只是把出現頻率太低的單詞（維度）去掉。

   不過，事實上我們通常會做更復雜一些的處理，例如，如果你是在做 sentiment analysis ，那么你通常會更加關注語氣很重的詞，比如 "bad"、"terrible"、"awesome" 等的重要性就比普通的詞要大，此時你可以為每一個維度設一個權重，例如，給 "bad" 設置權重 2 ，那么出現 3 次的話，向量在該維度對應的值就為 2*3 = 6 。當然這樣手工指定權重只在小范圍內可行，如果要為數百萬個維度指定權重，那顯然是不可能的，另一個稍微自動一點的辦法是 tf-idf 。

   tf 就是 Term Frequency ，就和剛才說的單詞出現的次數差不多，而 idf 則是 Inverse Document Frequency ，通常使用如下公式進行計算：

   

   這相當於自動計算出來的每個單詞的權重，其想法很簡單：如果在許多文檔中都出現的詞（例如虛詞、語氣詞等），它所包含的信息量通常會比較小，所以用以上的公式計算出來的權重也會相對較小；另一方面，如果單詞並不是在很多文檔中都出現的話，很有可能就是出現的那些文檔的重要特征，因此權重會相對大一些。

   前面說了一堆，其實主要是想要對 feature 做一些"預"處理，讓它更"好"一些，手工設置維度的權重固然是很人力，其實 tf-idf 也是在假定了原始 feature 是 document 、term 這樣的形式（或者類似的模型）的情況下才適用的（例如，在門禁之類的系統里，feature 可能有聲音、人臉圖像以及指紋等數據，就不能這樣來處理），因此也並不能說是一種通用的方法。

   然而，既然機器學習的算法可以在不考慮原始數據的來源和意義的情況下工作，那么 feature 的處理應該也可以。事實也是如此，包括 feature selection 和 dimensionality reduction 兩個 topic 都有許多很經典的算法。前者通常是選出重要的 feature 的維度（並拋棄不重要的維度），而后者則是更廣泛意義上的把一個高維的向量映射為一個低維向量，亦即"降維"，得到的結果 feature 的值已經不一定是原始的值，也可以把 feature selection 看作是 dimensionality reduction 的一種特殊情況。舉一個例子，tf-idf 實際上不算 feature selection ，因為它（通常）並沒有丟棄低權值的維度，並且處理過后的特征的每個維度都被乘上了一個權值，不再是原來的值了；但是它卻可以被看作一種降維，雖然嚴格意義上來說維度並沒有"降低"。簡單來說降維可以看作一個函數，其輸入是一個 D 維的向量，輸出是一個 M 維的向量。

   按照機器學習的方法的一貫作風，我們需要定義目標函數並進行最優化。不同的目標也就導致了不同的降維算法，這也正是今天要講的話題。

   然而，我們的目的究竟是什么呢？一個比較直接的問題是原始的數據量維度太高，我們無法處理，因此需要降維，此時我們通常希望在最大限度地降低數據的維度的前提下能夠同時保證保留目標的重要的信息，就好比在做有損的數據壓縮時希望壓縮比盡量大但是質量損失不要太多。於是問題又轉化為如何衡量對信息的保留程度了。

   一個最直接的辦法就是衡量 reconstruction error ，即

   

   其中  是  所對應的低維表示再重新構造出來的高維形式，就相當於是壓縮之后解壓出來的結果，不過雖然有許多壓縮方法都是無損的，就是說這個差值會等於零，但是大部分降維的結果都是有損的。不過我們仍然希望把上面的 reconstruction error 最小化。

   另外一種方式是簡單地使用 variance 來衡量所包含信息量，例如，我們要把一些 D 維的向量降為 1 維，那么我們希望這一維的 variance 達到最大化，亦即：

   

   其中  是降維函數。推而廣之，如果是降為 2 維，那么我希望第 2 維去關注第 1 維之外的 variance ，所以要求它在與第一維垂直的情況下也達到 variance 最大化。以此類推。

   然而，當我們把降維函數  限定維線性的時候，兩種途徑會得到同樣的結果，就是被廣泛使用的 Principal Components Analysis(PCA) 。PCA 的降維函數是線性的，可以用一個  維的矩陣  來表示，因此，一個 D 維的向量  經過線性變換  之后得到一個 M 維向量，就是降維的結果。把原始數據按行排列為一個  維的矩陣  ，則  就是降維后的  維的數據矩陣，目標是使其 covariance 矩陣最大。在數據被規則化（即減去其平均值）過的情況下，協方差矩陣 （covariance）  ，當然矩陣不是一個數，不能直接最大化，如果我們采用矩陣的 Trace （亦即其對角線上元素的和）來衡量其大小的話，要對  求最大化，只需要求出  的特征值和特征向量，將 M 個最大的特征值所對應的特征向量按列排列起來組成線性變換矩陣  即可。這也就是 PCA 的求解過程，得到的降維矩陣 可以直接用到新的數據上。如果熟悉 Latent Semantic Analysis (LSA) 的話，大概已經看出 PCA 和 Singular Value Decomposition (SVD) 以及 LSA 之間的關系了。以下這張圖（引自《The Elements of Statistical Learning》）可以直觀地看出來 PCA 做了什么，對於一些原始為二維的數據，PCA 首先找到了 variance 最大的那一個方向：

   

   PCA 作為一種經典的降維方法，被廣泛地應用於機器學習、計算機視覺以及信息檢索等各個領域，其地位類似於聚類中的 K-means ，在現在關於降維等算法的研究中也經常被作為 baseline 出現。然而，PCA 也有一些比較明顯的缺點：一個就是 PCA 降維是線性變換，雖然線性變換計算方便，並且可以很容易地推廣到新的數據上，然而有些時候線性變換也許並不合適，為此有許多擴展被提出來，其中一個就是 Kernel PCA ，用 Kernel Trick 來將 PCA 推廣到非線性的情況。另外，PCA 實際上可以看作是一個具有 Gaussian 先驗和條件概率分布的 latent variable 模型，它假定數據的 mean 和 variance 是重要的特征，並依靠 covariance 最大化來作為優化目標，而事實上這有時候對於解決問題幫助並不大。

   一個典型的問題就是做聚類或分類，回想我們之前談過的 Spectral Clustering ，就是使用 Laplacian eigenmap 降維之后再做 K-means 聚類，如果問題定下來了要對數據進行區分的話，"目的"就變得明朗了一些，也就是為了能夠區分不同類別的數據，再考慮直觀的情況，我們希望如果通過降維把高維數據變換到一個二維平面上的話，可以很容易"看"出來不同類別的數據被映射到了不同的地方。雖然 PCA 極力降低 reconstruction error ，試圖得到可以代表原始數據的 components ，但是卻無法保證這些 components 是有助於區分不同類別的。如果我們有訓練數據的類別標簽，則可以用 Fisher Linear Discriminant Analysis 來處理這個問題。

   同 PCA 一樣，Fisher Linear Discriminant Analysis 也是一個線性映射模型，只不過它的目標函數並不是 Variance 最大化，而是有針對性地使投影之后屬於同一個類別的數據之間的 variance 最小化，並且同時屬於不同類別的數據之間的 variance 最大化。具體的形式和推導可以參見《Pattern Classification》這本書的第三章 Component Analysis and Discriminants 。

   當然，很多時候（比如做聚類）我們並不知道原始數據是屬於哪個類別的，此時 Linear Discriminant Analysis 就沒有辦法了。不過，如果我們假設原始的數據形式就是可區分的的話，則可以通過保持這種可區分度的方式來做降維，MDS 是 PCA 之外的另一種經典的降維方法，它降維的限制就是要保持數據之間的相對距離。實際上 MDS 甚至不要求原始數據是處在一個何種空間中的，只要給出他們之間的相對"距離"，它就可以將其映射到一個低維歐氏空間中，通常是三維或者二維，用於做 visualization 。

   不過我在這里並不打算細講 MDS ，而是想說一下在 Spectral Clustering 中用到的降維方法 Laplacian Eigenmap 。同 MDS 類似，LE 也只需要有原始數據的相似度矩陣，不過這里通常要求這個相似度矩陣  具有局部性質，亦即只考慮局部領域內的相似性，如果點  和  距離太遠的話， 應該為零。有兩種很直接的辦法可以讓普通的相似度矩陣具有這種局部性質：

   1. 通過設置一個閾值，相似度在閾值以下的都直接置為零，這相當於在一個   -領域內考慮局部性。
   2. 對每個點選取 k 個最接近的點作為鄰居，與其他的點的相似性則置為零。這里需要注意的是 LE 要求相似度矩陣具有對稱性，因此，我們通常會在    屬於    的 k 個最接近的鄰居且/或反之的時候，就保留    的值，否則置為零。

   構造好  之后再來考慮降維，從最簡單的情況開始，即降到一維  ，通過最小化如下目標函數來實現：

   

   從直觀上來說，這樣的目標函數的意義在於：如果原來  和  比較接近，那么  會相對比較大，這樣如果映射過后  和  相差比較大的話，就會被權重  放大，因此最小化目標函數就保證了原來相近的點在映射過后也不會彼此相差太遠。

   令  為將  的每一行加起來所得到的對角陣，而  ，被稱作是拉普拉斯矩陣，通過求解如下的特征值問題

   

   易知最小的那個特征值肯定是 0 ，除此之外的最小的特征值所對應的特征向量就是映射后的結果。特征向量是一個 N 維列向量，將其旋轉一下，正好是 N 個原始數據降到一維之后的結果。

   類似地推廣到 M 維的情況，只需要取除去 0 之外的最小的 M 個特征值所對應的特征向量，轉置之后按行排列起來，就是降維后的結果。用 Matlab 代碼寫出來如下所示（使用了 knn 來構造相似度矩陣，並且沒有用 heat kernel ）：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

   % Laplacian Eigenmap ALGORITHM (using K nearest neighbors) % % [Y] = le(X,K,dmax) % % X = data as D x N matrix (D = dimensionality, N = #points) % K = number of neighbors % dmax = max embedding dimensionality % Y = embedding as dmax x N matrix     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     function [Y] = leigs(X,K,d)     [D,N] = size(X); fprintf(1,'LE running on %d points in %d dimensions\n',N,D);         % STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS fprintf(1,'-->Finding %d nearest neighbours.\n',K);     X2 = sum(X.^2,1); distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;     [sorted,index] = sort(distance); neighborhood = index(2:(1+K),:);             % STEP2: Construct similarity matrix W fprintf(1,'-->Constructing similarity matrix.\n');     W = zeros(N, N); for ii=1:N W(ii, neighborhood(:, ii)) = 1; W(neighborhood(:, ii), ii) = 1; end     % STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF L fprintf(1,'-->Computing embedding.\n');     D = diag(sum(W)); L = D-W;     % CALCULATION OF EMBEDDING options.disp = 0; options.isreal = 1; options.issym = 1; [Y,eigenvals] = eigs(L,d+1,0,options); Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0         fprintf(1,'Done.\n');     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   事實上，Laplacian Eigenmap 假設數據分布在一個嵌套在高維空間中的低維流形上， Laplacian Matrix  則是流形的 Laplace Beltrami operator 的一個離散近似。關於流形以及 Laplacian Eigenmap 這個模型的理論知識就不在這里做更多地展開了，下面看一個比較直觀的例子 Swiss Roll 。

   Swiss Roll 是一個像面包圈一樣的結構，可以看作一個嵌套在三維空間中的二維流形，如下圖所示，左邊是 Swiss Roll ，右邊是從 Swiss Roll 中隨機選出來的一些點，為了標明點在流形上的位置，給它們都標上了顏色。

   

   而下圖是 Laplacian Eigenmap 和 PCA 分別對 Swiss Roll 降維的結果，可以看到 LE 成功地把這個流形展開把在流形上屬於不同位置的點映射到了不同的地方，而 PCA 的結果則很糟糕，幾種顏色都混雜在了一起。

   

   另外，還有一種叫做 Locally Linear Embedding 的降維方法，它同 LE 一樣采用了流形假設，並假定平滑流形在局部具有線性性質，一個點可以通過其局部鄰域內的點重構出來。首先它會將下式最小化

   

   以求解出最優的局部線性重構矩陣  ，對於距離較遠的點  和  ， 應當等於零。這之后再把  當作已知量對下式進行最小化：

   

   這里  成了變量，亦即降維后的向量，對這個式子求最小化的意義很明顯了，就是要求如果原來的數據  可以以  矩陣里對應的系數根據其鄰域內的點重構出來的話，那么降維過后的數據也應該保持這個性質。

   經過一些變換之后得到的求解方法和 LE 類似，也是要求解一個特征值問題，實際上，從理論上也可以得出二者之間的聯系（LE 對應於  而 LLE 對應於 ），如果感興趣的話，可以參考Laplacian Eigenmaps for Dimensionality Reduction and Data Representation 這篇論文里的對比。下面是兩種方法在 Swiss Roll 數據上的結果，也是非常相似的：

   

   有一點需要注意的是，LE 和 LLE 都是非線性的方法，PCA 得到的結果是一個投影矩陣，這個結果可以保存下來，在之后對任意向量進行投影，而 LE 和 LLE都是直接得出了數據降維之后的結果，但是對於新的數據，卻沒有得到一個"降維函數"，沒有一個合適的方法得到新的數據的降維結果。所以，在人們努力尋求非線性形式擴展 PCA 的時候，另一些人則提出了 LE 和 LLE 的線性形式，分別叫做 Locality Preserving Projection 和 Neighborhood Preserving Embedding 。

   在 LPP 中，降維函數跟 PCA 中一樣被定為是一個線性變換，用一個降維矩陣  來表示，於是 LE 的目標函數變為

   

   經過類似的推導，最終要求解的特征值問題如下：

   

   得到的按照特征值從小到大排序的特征向量就組成映射矩陣  ，和 LE 不同的是這里不需要去掉第一個特征向量。另一點是在 LE 中的特征值是一個稀疏的特征值問題，在只需要求解最小的幾個特征值的時候可以比較高效地求解，而這里的矩陣在乘以  之后通常就不再稀疏了，計算量會變得比較大，這個問題可以使用 Spectral Regression 的方法來解決，參見 Spectral Regression: A Unified Approach for Sparse Subspace Learning 這篇 paper 。如果采用 Kernel Trick 再把 LPP 非線性化的話，又會回到 LE 。而 LLE 的線性版本 NPE 也是采用了類似的辦法來得到的，就不在這里多講了。

   另外，雖然 LE 是 unsupervised 的，但是如果訓練數據確實有標簽可用，也是可以加以利用的——在構造相似度矩陣的時候，屬於同一類別的相似度要大一些，而不同類別的相似度則會小一些。

   當然，除去聚類或分類之外，降維本身也是一種比較通用的數據分析的方法，不過有許多人批評降維，說得到的結果沒有意義，不用說非線性，就是最簡單的線性降維，除去一些非藏極端的特殊情況的話，通常將原來的分量線性組合一下都不會得到什么有現成的物理意義的量了。然而話也說回來，現在的機器學習幾乎都是更 prefer "黑盒子"式的方法吧，比如決策樹，各個分支對應與變量的話，它的決策過程其實人是可以"看到"或者說"理解"的，但是 SVM 就不那么"直觀"了，如果再加上降維處理，就更加"不透明"了。不過我覺得這沒什么不好的，如果只是靠可以清晰描訴出來的 rule 的話，似乎感覺神秘感不夠，沒法發展出"智能"來啊 ^_^bb 最后，所謂有沒有物理意義，其實物理量不過也都是人為了描述問題方便而定義出來的吧。

4. Hierarchical Clustering

   本系列不小心又拖了好久，其實正兒八經的 blog 也好久沒有寫了，因為比較忙嘛，不過覺得 Hierarchical Clustering 這個話題我能說的東西應該不多，所以還是先寫了吧（我准備這次一個公式都不貼 ）。Hierarchical Clustering 正如它字面上的意思那樣，是層次化的聚類，得出來的結構是一棵樹，如右圖所示。在前面我們介紹過不少聚類方法，但是都是"平坦"型的聚類，然而他們還有一個更大的共同點，或者說是弱點，就是難以確定類別數。實際上，（在某次不太正式的電話面試里）我曾被問及過這個問題，就是聚類的時候如何確定類別數。

   我能想到的方法都是比較 naive 或者比較不靠譜的，比如：

   1. 根據數據的來源使用領域相關的以及一些先驗的知識來進行估計——說了等於沒有說啊……
   2. 降維到二維平面上，然后如果數據形狀比較好的話，也許可以直觀地看出類別的大致數目。
   3. 通過譜分析，找相鄰特征值 gap 較大的地方——這個方法我只了解個大概，而且我覺得"較大"這樣的詞也讓它變得不能自動化了。

   當時對方問"你還有沒有什么問題"的時候我竟然忘記了問他這個問題到底有沒有什么更好的解決辦法，事后真是相當后悔啊。不過后來在實驗室里詢問了一下，得到一些線索，總的來說復雜度是比較高的，待我下次有機會再細說（先自己研究研究）。

   不過言歸正傳，這里要說的 Hierarchical Clustering 從某種意義上來說也算是解決了這個問題，因為在做 Clustering 的時候並不需要知道類別數，而得到的結果是一棵樹，事后可以在任意的地方橫切一刀，得到指定數目的 cluster ，按需取即可。

   聽上去很誘人，不過其實 Hierarchical Clustering 的想法很簡單，主要分為兩大類：agglomerative（自底向上）和 divisive（自頂向下）。首先說前者，自底向上，一開始，每個數據點各自為一個類別，然后每一次迭代選取距離最近的兩個類別，把他們合並，直到最后只剩下一個類別為止，至此一棵樹構造完成。

   看起來很簡單吧？  其實確實也是比較簡單的，不過還是有兩個問題需要先說清除才行：

5. 如何計算兩個點的距離？這個通常是 problem dependent 的，一般情況下可以直接用一些比較通用的距離就可以了，比如歐氏距離等。
6. 如何計算兩個類別之間的距離？一開始所有的類別都是一個點，計算距離只是計算兩個點之間的距離，但是經過后續合並之后，一個類別里就不止一個點了，那距離又要怎樣算呢？到這里又有三個變種：
   1. Single Linkage：又叫做 nearest-neighbor ，就是取兩個集合中距離最近的兩個點的距離作為這兩個集合的距離，容易造成一種叫做 Chaining 的效果，兩個 cluster 明明從"大局"上離得比較遠，但是由於其中個別的點距離比較近就被合並了，並且這樣合並之后 Chaining 效應會進一步擴大，最后會得到比較松散的 cluster 。
   2. Complete Linkage：這個則完全是 Single Linkage 的反面極端，取兩個集合中距離最遠的兩個點的距離作為兩個集合的距離。其效果也是剛好相反的，限制非常大，兩個 cluster 即使已經很接近了，但是只要有不配合的點存在，就頑固到底，老死不相合並，也是不太好的辦法。
   3. Group Average：這種方法看起來相對有道理一些，也就是把兩個集合中的點兩兩的距離全部放在一起求一個平均值，相對也能得到合適一點的結果。

   總的來說，一般都不太用 Single Linkage 或者 Complete Linkage 這兩種過於極端的方法。整個 agglomerative hierarchical clustering 的算法就是這個樣子，描述起來還是相當簡單的，不過計算起來復雜度還是比較高的，要找出距離最近的兩個點，需要一個雙重循環，而且 Group Average 計算距離的時候也是一個雙重循環。

   另外，需要提一下的是本文一開始的那個樹狀結構圖，它有一個專門的稱呼，叫做 Dendrogram，其實就是一種二叉樹，畫的時候讓子樹的高度和它兩個后代合並時相互之間的距離大小成比例，就可以得到一個相對直觀的結構概覽。不妨再用最開始生成的那個三個 Gaussian Distribution 的數據集來舉一個例子，我采用 Group Average 的方式來計算距離，agglomerative clustering 的代碼很簡單，沒有做什么優化，就是直接的雙重循環：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

   def do_clustering(nodes): # make a copy, do not touch the original list nodes = nodes[:] while len(nodes) > 1: print "Clustering [%d]..." % len(nodes) min_distance = float('inf') min_pair = (-1, -1) for i in range(len(nodes)): for j in range(i+1, len(nodes)): distance = nodes[i].distance(nodes[j]) if distance < min_distance: min_distance = distance min_pair = (i, j) i, j = min_pair node1 = nodes[i] node2 = nodes[j] del nodes[j] # note should del j first (j > i) del nodes[i] nodes.append(node1.merge(node2, min_distance))     return nodes[0]

   數據點又一千多個，畫出來的 dendrogram 非常大，為了讓結果看起來更直觀一點，我把每個葉節點用它本身的 label 來染色，並且向上合並的時候按照權重混合一下顏色，最后把圖縮放一下得到這樣的一個結果（點擊查看原圖）：

   或者可以把所有葉子節點全部拉伸一下看，在右邊對齊，似乎起來更加直觀一點：

   

   

   從這個圖上可以很直觀地看出來聚類的結果，形成一個層次，而且也在總體上把上個大類分開來了。由於這里我把圖橫過來畫了，所以在需要具體的 flat cluster 划分的時候，直觀地從圖上可以看成豎着划一條線，打斷之后得到一片"森林"，再把每個子樹里的所有元素變成一個"扁平"的集合即可。完整的 Python 代碼如下：

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106

   from scipy.linalg import norm from PIL import Image, ImageDraw     def make_list(obj): if isinstance(obj, list): return obj return [obj]     class Node(object): def __init__(self, fea, gnd, left=None, right=None, children_dist=1): self.__fea = make_list(fea) self.__gnd = make_list(gnd) self.left = left self.right = right self.children_dist = children_dist     self.depth = self.__calc_depth() self.height = self.__calc_height()     def to_dendrogram(self, filename): height_factor = 3 depth_factor = 20 total_height = int(self.height*height_factor) total_depth = int(self.depth*depth_factor) + depth_factor im = Image.new('RGBA', (total_depth, total_height)) draw = ImageDraw.Draw(im) self.draw_dendrogram(draw, depth_factor, total_height/2, depth_factor, height_factor, total_depth) im.save(filename)         def draw_dendrogram(self,draw,x,y,depth_factor,height_factor,total_depth): if self.is_terminal(): color_self = ((255,0,0), (0,255,0), (0,0,255))[int(self.__gnd[0])] draw.line((x, y, total_depth, y), fill=color_self) return color_self else: y1 = int(y-self.right.height*height_factor/2) y2 = int(y+self.left.height*height_factor/2) xc = int(x + self.children_dist*depth_factor) color_left = self.left.draw_dendrogram(draw, xc, y1, depth_factor, height_factor, total_depth) color_right = self.right.draw_dendrogram(draw, xc, y2, depth_factor, height_factor, total_depth)     left_depth = self.left.depth right_depth = self.right.depth sum_depth = left_depth + right_depth if sum_depth == 0: sum_depth = 1 left_depth = 0.5 right_depth = 0.5 color_self = tuple([int((a*left_depth+b*right_depth)/sum_depth) for a, b in zip(color_left, color_right)]) draw.line((xc, y1, xc, y2), fill=color_self) draw.line((x, y, xc, y), fill=color_self) return color_self         # use Group Average to calculate distance def distance(self, other): return sum([norm(x1-x2) for x1 in self.__fea for x2 in other.__fea]) \ / (len(self.__fea)*len(other.__fea))     def is_terminal(self): return self.left is None and self.right is None     def __calc_depth(self): if self.is_terminal(): return 0 return max(self.left.depth, self.right.depth) + self.children_dist     def __calc_height(self): if self.is_terminal(): return 1 return self.left.height + self.right.height     def merge(self, other, distance): return Node(self.__fea + other.__fea, self.__gnd + other.__gnd, self, other, distance)         def do_clustering(nodes): # make a copy, do not touch the original list nodes = nodes[:] while len(nodes) > 1: print "Clustering [%d]..." % len(nodes) min_distance = float('inf') min_pair = (-1, -1) for i in range(len(nodes)): for j in range(i+1, len(nodes)): distance = nodes[i].distance(nodes[j]) if distance < min_distance: min_distance = distance min_pair = (i, j) i, j = min_pair node1 = nodes[i] node2 = nodes[j] del nodes[j] # note should del j first (j > i) del nodes[i] nodes.append(node1.merge(node2, min_distance))     return nodes[0]

   agglomerative clustering 差不多就這樣了，再來看 divisive clustering ，也就是自頂向下的層次聚類，這種方法並沒有 agglomerative clustering 這樣受關注，大概因為把一個節點分割為兩個並不如把兩個節點結合為一個那么簡單吧，通常在需要做 hierarchical clustering 但總體的 cluster 數目又不太多的時候可以考慮這種方法，這時可以分割到符合條件為止，而不必一直分割到每個數據點一個 cluster 。

   總的來說，divisive clustering 的每一次分割需要關注兩個方面：一是選哪一個 cluster 來分割；二是如何分割。關於 cluster 的選取，通常采用一些衡量松散程度的度量值來比較，例如 cluster 中距離最遠的兩個數據點之間的距離，或者 cluster 中所有節點相互距離的平均值等，直接選取最"松散"的一個 cluster 來進行分割。而分割的方法也有多種，比如，直接采用普通的 flat clustering 算法（例如 k-means）來進行二類聚類，不過這樣的方法計算量變得很大，而且像 k-means 這樣的和初值選取關系很大的算法，會導致結果不穩定。另一種比較常用的分割方法如下：

7. 待分割的 cluster 記為 G ，在 G 中取出一個到其他點的平均距離最遠的點 x ，構成新 cluster H；
8. 在 G 中選取這樣的點 x' ，x' 到 G 中其他點的平均距離減去 x' 到 H 中所有點的平均距離這個差值最大，將其歸入 H 中；
9. 重復上一個步驟，直到差值為負。

   到此為止，我的 hierarchical clustering 介紹就結束了。總的來說，在我個人看來，hierarchical clustering 算法似乎都是描述起來很簡單，計算起來很困難（計算量很大）。並且，不管是 agglomerative 還是 divisive 實際上都是貪心算法了，也並不能保證能得到全局最優的。而得到的結果，雖然說可以從直觀上來得到一個比較形象的大局觀，但是似乎實際用處並不如眾多 flat clustering 算法那么廣泛。

10. Deciding the Number of Clusterings

    如何確定聚類的類別個數這個問題經常有人問我，也是一直以來讓我自己也比較困惑的問題。不過說到底其實也沒什么困惑的，因為這個問題本身就是一個比較 ill posed 的問題呀：給一堆離散的點，要你給出它們屬於幾個 cluster，這個基本上是沒有唯一解或者說是沒有合適的標准來衡量的。比如如果簡單地用每個類別里的點到類中心的距離之和來衡量的話，一下子就會進入到"所有的點都獨立成一類"這樣的尷尬境界中。

    但是 ill posed 也並不是一個很好的理由，因為我們其實大部分時候都是在處理 ill posed 的問題嘛，比如 Computer Vision 整個一個領域基本上就沒有啥問題是 well posed 的…… =.=bb，比如下面盜用一下 Bill Freeman 的 Slides 中的一張講 deblur 的問題的圖：

    Blurry image Sharp image Blur kernel

    		=
    		=
    		=

    從 Stata Center 的模糊照片找出清晰照片和模糊核的這個過程（特別對於計算機來說）就是非常 ambigous 的。（順便這個 Slides 本身也是很有意思的，推薦看一下。）

    所以么，還是讓我們先拋開各種借口，回到問題本身。當然這次的標題取得有點宏大，其實寫這篇日志的主要目的不過是想吐槽一下上一次 6.438 課的一道作業題而已……=.=bb

    總而言之呢，像 Kmeans 之類的大部分經典的聚類算法，都是需要事先指定一個參數K作為類別數目的。但是很多時候這個K值並不是那么好確定的，更麻煩的是，即使你使用很暴力的方法把某個范圍內的整數K全部都試一遍，通常也沒法知道哪個K才是最好的。

    所以呢，我們自然會問：那有沒有算法不需要指定類別數呢？當然有！最簡單的就是層次聚類，它會將你的數據點用一個樹形結構給連起來，可是我想要的是聚類結果啊！這也好辦，只要在合適的樹深度的地方把樹切開，變成一個一個的子樹，每個子樹就是一個cluster。不過問題就來了，究竟什么深度才是合適的深度呢？事實上不同的深度會產生不同的類別數目，這基本上把原來的選擇K的問題轉化成了選擇樹深度的問題。

    沒關系，我們還有其他貨，比如著名的 Mean Shift 算法，也是不需要指定類別數目就可以聚類的，而且聚類的結果也不是一棵樹那種奇奇怪怪的東西。具體 Mean Shift 算法是什么樣的今天就不在這里講了，不過它其實也有一個參數 Window Size 需要選擇——對，正如大家所料，這個 Window Size 其實也是會影響最終聚類的類別數目。換句話說，原來選擇K的問題現在被轉化為了選擇 Window Size 的問題。

    

    圖 2    Factor graph demonstration for Affinity Propagation of 3 points.

    然后我們進入貝葉斯推斷的領域，之前跟別人聊天的時候就聽說過有個叫做 Dirichlet Process 的東西，特別神奇，能夠自動地根據數據 adaptive 地確定合適的類別數目。我感覺這背后可能像其他算法一樣會隱藏着其他的參數需要調節，不過鑒於我對 DP 完全不懂，就不在這里說胡話了，感興趣的同學可以去研究一下。

    然后貝葉斯推斷聚類里的另一個成員就是 Affinity Propagation，其實是在講了 Loopy Belief Propagation 算法之后的一道作業題，給了一個 factor graph，讓把消息傳遞的公式推出來然后實現出來。看起來好像蠻容易的樣子，實際上折騰了好幾天，因為實現的時候碰到各種 tricky 的問題。

    AP 算法把聚類問題看成一個 MAP 推斷問題，假設我們有n個數據點，這里我們要求類別中心實際上是這些數據點中的一個子集，聚類的結果可以用這個binary variable 來表示，其中表示點被歸到點所代表的那個類別中。

    給定個隨機變量之后，接下來是要通過 local factor 來定義他們的 (unnormalized) joint distribution。首先我們要求已知一個 Affinity Matrix S，其中表示點和點之間的相似程度（這里  並不一定要求是對稱的），一個簡單的做法就是令

    然后整個 joint distribution 主要由兩類 factor 構成，在圖 (fig: 2) 中給出了 3 個數據點時候的 factor graph 的例子。其中橙色的 factor 定義為

    直觀地來講，該 factor 的第一部分是說每個點只能且必須被 assign 到一個 cluster；第二部分是講被assign 到的時候 factor 的值為 ，也就是說會偏向於 assign 到相似度較大的那個 cluster 去。

    而綠色的那些 factor 則定義為：

    這樣的 indicator factor 要求如果有任何點被 assign 到的話，那么必須要被 assign 到它自己。換句話說，每個 cluster 的中心必須是屬於它自己那個類的，這也是非常合理的要求。

    好了，模型的定義就是這么簡單，這樣我們會得到一個 loopy factor graph，接下來只要在這個 graph 上跑一下 Max-Product Algorithm 就可以得到最優的 ，從而得到聚類結果了，並且都不用指定類別數目什么的，因為合適的類別數目可以從數據中自動地推斷出來！

    簡直太神奇了，我都不敢相信這是真的！結果，果然童話里都是騙人的……因為隱藏在背后的還有許多細節需要處理。首先有一些 trick 來降低算法的復雜程度，比如說，由於所有的變量都是 binary 的，因此在消息傳遞的過程中我們可以只傳遞等於 1 和等於 0 兩種狀態時的消息差；然后由於數值下溢的原因，一般都會對所有的消息取log ，這樣 Max-Product 算法會變成 Max-Sum（或者 Min-Sum）；然后有些消息是不做任何操作就直接 literally 傳到接下來的節點的，這樣的消息可以從中間步驟中去掉，最后化簡之后的結果會只留下兩種消息，在經典的 Affinity Propagation 算法中分別被成為 responsibility 消息和 availability 消息，可以有另外的 intuitive 的解釋，具體可以參見 (Frey & Dueck, 2007)。

    

    圖 3    Clustering result of Affinity Propagation.

    不過這些還是不太夠，因為 Loopy BP 的收斂性其實是沒有什么保證的，隨便實現出來很可能就不收斂。比如說，如果直接進行 parallel 地 message updating 的話，很可能就由於 update 的幅度過大導致收斂困難，因此需要用 dampen 的方法將 fixed-point 的式子

    改為

    然后用這個來做 updating。當然也可以用 0.5 以外的其他 dampen factor。不過這樣還是不夠的，反正我折騰了好久最后發現必須把 parallel updating 改成 in-place updating 才能收斂，這樣一來更新的幅度就更小了一些，不過 dampen 仍然是需要的。另外就是一同討論的另外的同學還發現實現中消息更新的順序也是極端重要……把順序變一下就立馬從完全不 work 變成結果超好了。-_-!!於是突然想起來 CV 課上老師拿 Convolutional Network 開玩笑的時候說，這個算法效果超級好，但是就是太復雜了，在發明這個算法的那個 lab 所在地的方圓 50 米之外完全沒人能成功訓練 Convolutional Network……

    

    圖 4    Clustering result of Affinity Propagation, with large diagonal values for the affinity matrix.

    此外收斂性的判斷也是有各種方法，比如消息之差的絕對值、或者是連續多少輪迭代產生的聚類中心都沒有發生變化等等。即使在收斂之后，如何確定 assignment 也並不是顯而易見的。因為 Max-Product 算法得到的是一些 Max-Marginal，如果全局最優的 configuration 不是 unique 的話，general 地是不能直接局部地通過各自的 Max-Marginal 來確定全局 configuration 的，於是可能需要實現 back-tracking，總之就是動態規划的東西。不過也有簡單的 heuristic 方法可以用，比如直接通過的Max-Marginal 來確定點到底是不是一個 cluster center，確定了 center 之后則再可以通過 message 或者直接根據 Affinity 矩陣來確定最終的類別 Assignment。圖 (fig: 3) 展示了一個二維情況下 100 個點通過 AP 聚類的結果。

    不過，你有沒有發現一點不對勁的地方？一切似乎太完美了：什么參數都不用給，就能自動確定類別數？果然還是有問題的。事實上，如果你直接就去跑這個算法，肯定得不到圖上的結果，而是會得到一個n個類別的聚類結果：所有的點都成為了中心並且被歸類到了自己。這樣的結果也是可以理解的，因為我們計算 Affinity 的方式（距離的相反數）導致每個點和自己的 Affinity 是最大的……所以呢，為了解決這個問題，我們需要調整一個參數，就是S矩陣的對角線上的值，這個值反應了每個點自己成為類別中心的偏好程度，對的，如大家所料，這個就是背后的隱藏參數，直接影響類別數目。

    

    圖 5    Clustering result of Affinity Propagation, with small diagonal values for the affinity matrix.

    如圖 (fig: 4) 和 (fig: 5) 分別是把對角線上的元素值設置得比較大和比較小的情況，可以看到前者產生了更多的類別，因為大家都比較傾向於讓自己成為類別中心；而后者則只產生了兩個類，因為大家都很不情願……（我發現這種畫一些線連到類別中心的 visualization 方法似乎會給人產生很強的暗示感覺，於是各種聚類結果都看起來好像很"對"的樣子）

    所以呀，到最后好像還是竹籃打水一場空的樣子，怎么繞怎么繞好像也繞不開類別數目的這個參數。不過好像也並不是完全沒有任何進展的，因為雖然各種算法都或多或少地需要指定一些參數，但是這些參數從不同的角度來闡釋對應的問題，比如說在 AP 中，一般可以直接將  的對角線元素設置為其他所以元素的中位數，通常就能得到比較合理的結果；有時候對於特定的實際問題來說，從某一個角度來進行參數選擇可能會有比較直觀的 heuristic 可以用呢。所以說在是否需要指定參數這個問題上，還是不用太鑽牛角尖了啊。

    最后我還想提一下，最近看到的 lab 的兩篇 NIPS 文章 (Canas, Poggio & Rosasco, 2012),(Canas & Rosasco, 2012)，里面的視角也是很有意思的。

    我們知道，比如說，KMeans 算法的目標函數是這樣定義的：

    其中是一個由K個點組成的集合，而一個點到的距離被定義為到S中所有點的距離中最小的那一個值。這樣可以看成是用的這K個點去近似原本的n個點所帶來的誤差。然而這樣的目標函數對於選擇合適的K並沒有什么指導意義，因為隨着K的增大我們得到的最優的會使得目標函數越來越小。

    但是我們可以將這個目標函數推廣一下，類似於 Supervised Learning Theory 中那樣引入概率空間。具體地來說，我們假設數據點是從一個未知的概率分布中采樣而得到的，並且目標函數也從對n個數據點的近似推廣到整個生成流形的近似：

    當然由於是未知的，所以是沒法計算的，但是我們可以從 data sample 得到的來對它進行近似，並得到 bound 之類的。直觀地來講，此時對數值K的選擇將會影響到模型空間的復雜度，於是會出現一個trade-off，於是從這個角度下去探討"最優的K值"就變成一個很合理的問題了。感興趣的同學可以詳細參考 paper 以及里面的參考文獻，都在 arXiv 上可以下載到的。

    References

    1. Canas, G. D., & Rosasco, L. (2012). Learning Probability Measures with respect to Optimal Transport Metrics. In   NIPS.
    2. Canas, G. D., Poggio, T., & Rosasco, L. (2012). Learning Manifolds with K-Means and K-Flats. In   NIPS.
    3. Frey, B. J., & Dueck, D. (2007). Clustering by passing messages between data points.Science,   315, 972–977.