【機器學習】K-means聚類演算法初探
資料聚類是對於靜態資料分析的一門技術，在許多領域內都被廣泛地應用，包括機器學習、資料探勘、模式識別、影象分析、資訊檢索以及生物資訊等。聚類是把相似的物件通過靜態分類的方法分成不同的組別或者更多的子集，這樣讓在同一個子集中的成員物件都有相似的一些屬性，常見的包括在座標系中更加短的空間距離等。
我們拿2維特徵的例項作為例子。我們以這個2維特徵向量作為座標，在一個2維空間中用點標註出這些例項，如圖1所示，這裡是隨機生成的100個例項。圖1中帶顏色的方框同樣是例項，在之後介紹的K-means演算法中，這些方框是初始被隨機選擇出來的聚類中心點。
圖1 在二維空間中100個具有兩個特徵的例項
我們希望聚類演算法能夠將特徵相近的例項聚整合一個集合，最後形成多個由特徵相近的例項聚整合的聚類。如圖2
圖2 聚成10個類的效果
通常，我們將聚類分析作為一種無監督機器學習演算法來看待。與有監督學習不同，如常見的分類問題，我們有標註好分類類別的訓練資料，通過這些資料訓練出一個模型來對新的資料的分類進行預測。在聚類分析中，我們的資料是沒有標註類別的，我們通過資料中例項的特徵值相似程度將相似的例項劃分到一類中。
在接下來的內容中，我們將介紹一種簡單的聚類演算法----K-means，第二部分它的原理、過程以及可能存在的問題和解決方案。第三部分介紹
1 K-Means演算法介紹
1.1 演算法初探
K-means是一種基於距離的迭代式演算法[1]。它將n個觀察例項分類到k個聚類中，以使得每個觀察例項距離它所在的聚類的中心點比其他的聚類中心點的距離更小。
其中，距離的計算方式可以是歐式距離（2-norm distance），或者是曼哈頓距離（Manhattan distance,1-norm distance）或者其他。這裡我們使用歐式距離。
要將每個觀測例項都歸類到距離它最近的聚類中心點，我們需要找到這些聚類中心點的具體位置。然而，要確定聚類中心點的位置，我們必須知道該聚類中包含了哪些觀測例項。這似乎是一個“先有蛋還是先有雞”的問題。從理論上來說，這是一個NP-Hard
但是，我們可以通過啟發式的演算法來近似地解決問題。找到一個全域性最優的方案是NP-hard問題，但是，降低問題的難度，如果能夠找到多個區域性最優方案，然後通過一種評價其聚類結果優劣的方式，把其中最優的解決方案作為結果，我們就可以得到一個不錯的聚類結果[2]。
這就是為什麼我們說K-means演算法是一個迭代式的演算法。演算法[2]的過程如下：
1）所有的觀測例項中隨機抽取出k個觀測點，作為聚類中心點，然後遍歷其餘的觀測點找到距離各自最近的聚類中心點，將其加入到該聚類中。這樣，我們就有了一個初始的聚類結果，這是一次迭代的過程。
2）我們每個聚類中心都至少有一個觀測例項，這樣，我們可以求出每個聚類的中心點（means），作為新的聚類中心，然後再遍歷所有的觀測點，找到距離其最近的中心點，加入到該聚類中。然後繼續執行2）。
3）如此往復2），直到前後兩次迭代得到的聚類中心點一模一樣。
這樣，演算法就穩定了，這樣得到的k個聚類中心，和距離它們最近的觀測點構成k個聚類，就是我們要的結果。
實驗證明，演算法試可以收斂的[2]。
計算聚類的中心點有三種方法如下：
1）Minkowski Distance 公式 —— λ 可以隨意取值，可以是負數，也可以是正數，或是無窮大。
公式（1）
2）Euclidean Distance 公式 —— 也就是第一個公式 λ=2 的情況
公式（2）
3）CityBlock Distance 公式 —— 也就是第一個公式 λ=1 的情況
公式（3）
這三個公式的求中心點有一些不一樣的地方，我們看下圖（對於第一個 λ 在 0-1之間）。
（1）Minkowski Distance （2）Euclidean Distance （3） CityBlock Distance
上面這幾個圖的大意是他們是怎麼個逼近中心的，第一個圖以星形的方式，第二個圖以同心圓的方式，第三個圖以菱形的方式。我們採用第二種方式，在實際實現中，我們並沒有開根號。
我們演算法十分的簡單。下面對100個隨機生成的2維特徵的觀測例項在2維空間中執行這個K-means演算法的過程用圖表示出來。整個過程用了10次迭代後收斂。假設k=10。
100個觀測點的初始分佈如圖1所示。方塊為隨機挑出的10個聚類中心點。
圖3顯示了將各個觀測值劃分到距離最近的聚類後的結果。
圖3
然後根據每個聚類的點的情況，計算出新的中心點，得到新的聚類中心，再以此往復，直到達到收斂的條件。
在第9次迭代和第10次迭代中，聚類中心沒有任何變化，所以演算法收斂，結束。從圖5中可以看出，聚類中心在這10次迭代中位置變化的過程。黑色方塊為初始位置，白色為收斂後位置，其餘同一顏色的為中間過程中的位置。
圖5 聚類中心位置變化的過程
那麼，如何評價一個聚類結果呢？我們計算所有觀測點距離它對應的聚類中心的距離的平方和即可，我們稱這個評價函式為evaluate(C)。它越小，說明聚類越好。
1.2 K-means的問題及解決方案
K-means演算法非常簡單，然而卻也有許多問題。
1）首先，演算法只能找到區域性最優的聚類，而不是全域性最優的聚類。而且演算法的結果非常依賴於初始隨機選擇的聚類中心的位置。我們通過多次執行演算法，使用不同的隨機生成的聚類中心點執行演算法，然後對各自結果C通過evaluate(C)函式進行評估，選擇多次結果中evaluate(C)值最小的那一個。
2）關於初始k值選擇的問題。首先的想法是，從一個起始值開始，到一個最大值，每一個值執行k-means演算法聚類，通過一個評價函式計算出最好的一次聚類結果，這個k就是最優的k。我們首先想到了上面用到的evaluate(C)。然而，k越大，聚類中心越多，顯然每個觀測點距離其中心的距離的平方和會越小，這在實踐中也得到了驗證。第四節中的實驗結果分析中將詳細討論這個問題。
3）關於效能問題。原始的演算法，每一次迭代都要計算每一個觀測點與所有聚類中心的距離。有沒有方法能夠提高效率呢？是有的，可以使用k-d tree或者ball tree這種資料結構來提高演算法的效率。特定條件下，對於一定區域內的觀測點，無需遍歷每一個觀測點，就可以把這個區域內所有的點放到距離最近的一個聚類中去。這將在第三節中詳細地介紹。
2 實現
2.1 k-d tree 與 ball tree
k-d tree[5]或者ball tree[4]這個資料結構在K近鄰演算法中被用到來提高演算法的效率，在K-means中能否使用這些資料結構呢？當然能，而且，它們用在這裡更加高效。在這裡，我們使用ball tree來優化演算法的效能。我們首先介紹k-d tree。
1）k-d tree[5]
把n維特徵的觀測例項放到n維空間中，k-d tree每次通過某種演算法選擇一個特徵(座標軸)，以它的某一個值作為分界做超平面，把當前所有觀測點分為兩部分，然後對每一個部分使用同樣的方法，直到達到某個條件為止。
上面的表述中，有幾個地方下面將會詳細說明：（1）選擇特徵（座標軸）的方法 （2）以該特徵的哪一個為界 （3）達到什麼條件演算法結束。
(1)選擇特徵的方法
計算當前觀測點集合中每個特徵的方差，選擇方差最大的一個特徵，然後畫一個垂直於這個特徵的超平面將所有觀測點分為兩個集合。
（2)以該特徵的哪一個值為界 即垂直選擇座標軸的超平面的具體位置。
第一種是以各個點的方差的中值（median）為界。這樣會使建好的樹非常地平衡，會均勻地分開一個集合。這樣做的問題是，如果點的分佈非常不好地偏斜的，選擇中值會造成連續相同方向的分割，形成細長的超矩形(hyperrectangles)。
替代的方法是計算這些點該座標軸的平均值，選擇距離這個平均值最近的點作為超平面與這個座標軸的交點。這樣這個樹不會完美地平衡，但區域會傾向於正方地被劃分，連續的分割更有可能在不同方向上發生。
（3）達到什麼條件演算法結束
實際中，不用指導葉子結點只包含兩個點時才結束演算法。你可以設定一個預先設定的最小值，當這個最小值達到時結束演算法。
圖 6 一個k-d tree劃分二維空間
圖6中，星號標註的是目標點，我們在k-d tree中找到這個點所處的區域後，依次計算此區域包含的點的距離，找出最近的一個點（黑色點），如果在其他region中還包含更近的點則一定在以這兩個點為半徑的圓中。假設這個圓如圖中所示包含其他區域。先看這個區域兄弟結點對應區域，與圓不重疊；再看其雙親結點的兄弟結點對應區域。從它的子結點對應區域中尋找（圖中確實與這個雙親結點的兄弟結點的子結點對應區域重疊了）。在其中找是否有更近的結點。
k-d tree的優勢是可以遞增更新。新的觀測點可以不斷地加入進來。找到新觀測點應該在的區域，如果它是空的，就把它新增進去，否則，沿著最長的邊分割這個區域來保持接近正方形的性質。這樣會破壞樹的平衡性，同時讓區域不利於找最近鄰。我們可以當樹的深度到達一定值時重建這棵樹。
然而，k-d tree也有問題。矩形並不是用到這裡最好的方式。偏斜的資料集會造成我們想要保持樹的平衡與保持區域的正方形特性的衝突。另外，矩形甚至是正方形並不是用在這裡最完美的形狀，由於它的角。如果圖6中的圓再大一些，即黑點距離目標點點再遠一些，圓就會與左上角的矩形相交，需要多檢查一個區域的點，而且那個區域是當前區域雙親結點的兄弟結點的子結點。
為了解決上面的問題，我們引入了ball tree。
2）ball tree[4]
解決上面問題的方案就是使用超球面而不是超矩形劃分區域。使用球面可能會造成球面間的重疊，但卻沒有關係。ball tree就是一個k維超球面來覆蓋這些觀測點，把它們放到樹裡面。圖7（a)顯示了一個2維平面包含16個觀測例項的圖,圖7（b）是其對應的ball tree，其中結點中的數字表示包含的觀測點數。
圖 7 ball tree對二維平面的劃分和ball tree
不同層次的圓被用不同的風格畫出。樹中的每個結點對應一個圓，結點的數字表示該區域保含的觀測點數，但不一定就是圖中該區域囊括的點數，因為有重疊的情況，並且一個觀測點只能屬於一個區域。實際的ball tree的結點儲存圓心和半徑。葉子結點儲存它包含的觀測點。
使用ball tree時，先自上而下找到包含target的葉子結點，從此結點中找到離它最近的觀測點。這個距離就是最近鄰的距離的上界。檢查它的兄弟結點中是否包含比這個上界更小的觀測點。方法是：如果目標點距離兄弟結點的圓心的距離大於這個圓的圓心加上前面的上界的值，則這個兄弟結點不可能包含所要的觀測點。（如圖8）否則，檢查這個兄弟結點是否包含符合條件的觀測點。
圖 8 點與超圓
那麼，ball tree的分割演算法是什麼呢？
選擇一個距離當前圓心最遠的觀測點i1，和距離i1最遠的觀測點 i2，將圓中所有離這兩個點最近的觀測點都賦給這兩個簇的中心，然後計算每一個簇的中心點和包含所有其所屬觀測點的最小半徑。對包含n個觀測點的超圓進行分割，只需要線性的時間。
與k-d tree一樣，如果結點包含的觀測點到達了預先設定的最小值，這個頂點就可以不再分割了。
2.2 在K-means演算法中使用ball tree
在k-means演算法中使用k-d tree和ball tree效率更高，因為在此所有的觀測點都是一起處理，而在k近鄰中，測試的觀測點單獨處理。
首先，根據所有的觀測點建立一個包含它們的k-d tree 或者是ball tree。
在k-means演算法的每一次迭代中，設定的聚類中心集合為C={Ci，i= 1,...,k},每個簇中心有一個對應其歸屬它的觀測點集合，初始為空。從根結點開始遍歷樹，直到找到葉子結點為止，判斷其中的觀測點距離哪一個中心近，然後賦給那個中心。有可能在瀏覽一個內部結點時，它所包含的點完全包含在一個聚類的區域內。這時只需要直接把其下所有的觀測點加入到該聚類中即可。
我們可以在每個簇的中心點儲存一個向量，儲存聚類中點各個座標之和和點的個數。在計算其中心點時只需一除即可。
在多維空間中，如何判斷一個內部結點包含的觀測點全部落在一個聚類中呢？ 我們採用下面描述的方法。Stackoverflow上的討論見這裡：點選這裡
如果一個空間容器比如這裡的超球面與某一簇中心的最大距離小於其與其它所有簇中心的最小距離，即：
公式（4）
其中，表示聚類的中心點，如果上式成立，則可推出：
公式（5）
即此容器中的任意觀測點距離這個的距離都小於距離其他聚類中心點的距離。此時，可以斷定這個容器中的所有觀測點都在這個所在的聚類中。
在我們這個環境即ball tree中，任意點與超球面的最大距離為其與圓心距離加圓的半徑，最小距離為其與圓心距離減去圓的半徑。在實際的執行中，這個策略可以明顯地提高程式執行的效能。經大量的實驗得出，在對ball tree中的結點訪問時，能夠避免繼續向下探索的次數佔所有訪問結點次數平均達到17.8764%。每次避免訪問其子樹時，就大量減少了訪問總結點的時間。
3 實驗結果分析
我們從k=2開始，一直呼叫K-means演算法（每次呼叫10次，選擇evaluate最小那一次聚類的結果返回）直到k=14.我們記錄下每個k聚類後評估函式evaluate(Ck)的值，然後畫出其變化的曲線。理論上，隨著k的增大，evaluate(Ck)的值一定是非遞增的。因為中心越多，各個觀測點距離中心的平均距離越近。所以，我們通過另一種方式來確定k。我們在一個要聚類資料集上我們按照上面的方法執行多次，其evaluate關於k的趨勢圖如圖9所示：
圖9 多次執行上述過程的evaluate關於k的變化圖
通過我們的多次試驗（上面只給出了3次的evaluate變化過程），多次結果的共同特徵是，每次k=9~k=10間，evaluate函式的計算值都有極小的變化，而k<9時變化又極大，在k=9之後明顯下降變緩。所以，我們k=9或者k=10作為初始的聚類個數，即k-means中的k比較好。根據MDL原則，這裡選擇k = 9。
上面的例子中以實際的過程提供了一種選擇k的思路，就是根據evaluate(C)值下降的“拐點”和MDL原則來選擇k。另外，還有許多更科學的方法選擇k，見這裡：點選這裡