漫談 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 二、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.