上一次我們瞭解了一個最基本的 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) 和中位數 (

然而並不是所有的資料都能滿足這樣的要求，對於數值型別的特徵，比如身高，可以很自然地用這樣的方式來處理，但是類別 (categorical) 型別的特徵就不行了。舉一個簡單的例子，如果我現在要對犬進行聚類，並且希望直接在所有犬組成的空間中進行，k-means 就無能為力了，因為歐氏距離

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

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

除此之外，由於中心點是在已有的資料點裡面選取的，因此相對於 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 。