K-Means是最為經典的無監督聚類（Unsupervised Clustering）演算法，其主要目的是將n個樣本點劃分為k個簇，使得相似的樣本盡量被分到同一個聚簇。K-Means衡量相似度的計算方法為歐氏距離（Euclid Distance）。

本文將會介紹以下幾個部分的內容：

  1. K-Means迭代求解
  2. K-Means缺點和優化
  3. Speed up K-Means with Random Approximation
  4. 實驗部分

1. K-Means迭代求解

數據集 D 有 n 個樣本點 {x_1, x_2, cdots, x_n} ，假設現在要將這些樣本點聚類為 k 個簇， k 個簇中心為 {mu_1, mu_2, cdots, mu_k} 。定義指示變數 gamma_{ij} in {0,1} ，如果第 i 個樣本被聚類到第 k 個簇，那麼 gamma_{ij} = 1 ，否則 gamma_{ij} = 0 。在K-Means中每個樣本只能屬於一個簇，這是Hard-Clustering，也稱為Non-fuzzy Clustering，與之對應的是Soft-Clustering或者Fuzzy-Clustering，指的是一個樣本可以被歸為多個簇。在K-Means裡面有 sum_{j}gamma_{ij} = 1 ，那麼K-Means的優化目標為：

J(gamma, mu) = sum_{i = 1}^nsum_{j=1}^kgamma_{ij}Vert x_i - mu_jVert_2^2

目標需要優化的變數有 gamma_{ij} 和 mu_j ，由於 gamma_{ij} 是01離散變數，優化起來是NP-Hard問題。所以採用啟發式的求解，即Lloyds algorithm。求解過程分為Assignment和Update兩步。

Assignment Step:

gamma_{ij} = egin{cases} 1,& j = argmin_{j^1} ||mathbf x_i - mu_{j^1}||^2 \ 0, & ext{otherwise} end{cases}

Update Step:

mu_j = frac{sum_{i=1}^n gamma_{ij}mathbf x_i}{sum_{i=1}^n gamma_{ij}}

在Assignment Step，把樣本歸為最近的簇中心所代表的簇，很明顯，這一步更新會減小目標函數的值；在Update Step中，更新聚簇中心，用這個簇的所有樣本的中心點當作新的聚簇中心，不難證明這一步是最小化簇內樣本到聚簇中心歐式距離平方和的最優解。所以，按照上述兩步，只要有 gamma_{ij} 和 mu_j被更新了，那麼目標函數值肯定會下降。由於將 n 個樣本點分配到 k 個聚簇中心最多隻有 k^n 個劃分方案，每一種劃分方案對應一個目標函數值。採用上述兩步演算法，每步更新都會使得目標函數值減少，故上述演算法流程肯定會收斂。但是由於是採用的啟發式演算法，不一定能保證收斂到全局最優解。

事實上，K-Means的兩步更新可以看作是EM演算法的兩步，即Expectation和Maximization兩步，這裡不再詳細介紹。

2. K-Means缺點和優化

K-Means有以下幾個缺點：

  1. 需要確定K的值。K的取值需要事先確定，然而在無監督聚類任務上，由於並不知道數據集究竟有多少類別，所以很難確定K的取值。舉例來說，利用K-Means做 彩色圖像分割（Segmentation）任務，如果圖像只有一個物體和背景，那麼很容易確定將K設置為2，並且較容易得到比較好的分割效果，但是對於一副有更多物體的圖像，K的設定就很難了，因為我們並不總能人工去判斷這幅圖像裡面有幾個類別。對於文檔聚類更是如此，比如一批新聞 數據集，可以按照主題聚類，但是有多少主題呢，這個不能事先獲得。
  2. 對異常點敏感。K-Means很容易受到異常點（outliers）的影響，由於K-Means在Update Step取得是簇內樣本均值，那麼就會很容易受到異常點的影響，比如某個簇內樣本在某個維度上的值特別大，這就使得聚簇中心偏向於異常點，從而導致不太好的聚類效果。
  3. 凸形聚類。K-Means由於採用歐氏距離來衡量樣本之間相似度，所以得到的聚簇都是凸的，就不能解決「S型」數據分佈的聚類，這就使得K-Means的應用範圍受限，難以發現數據集中一些非凸的性質。
  4. 聚簇中心初始化，收斂到局部最優，未考慮密度分佈等等。

對於以上給出的很多問題，有很多相應的解決方案，從而得到了很多K-Means衍生出來的演算法以及一些其餘的聚類演算法，下面大致列舉一下：

K-Medians。不同於K-Means使用歐式距離計算目標函數，K-Medians使用曼哈頓距離（Manhattan Distance），對應的則是 l1 範數。在目標求解的Update過程中，用簇內樣本每個維度的中位數來當作聚簇中心在該維度的值（K-Means使用的是平均值），K-Medians的Median就是這個意思；在Assignment過程中，使用曼哈頓距離將每個樣本劃歸為最近的聚簇。K-Medians可以很好地解決異常點問題，比如Figure 1所示，ABCD是一個簇內的四個點，D是異常點，現在執行Update Step更新聚簇中心，使用K-Means得到紅色的聚簇中心L2，使用K-Medians得到藍色的聚簇中心L1。不難發現，在異常點D的影響下，K-Medians更魯棒，能較好的代表ABC三點。

Figure 1 : K-Medians解決outlier示意圖

K-Medoids。在K-Means中，聚簇中心不一定是樣本點，而是樣本點的中心。但是在K-Medoids中，聚簇中心就是樣本點。比如在初始化過程中，K-Means是隨機初始化k個聚簇中心，K-Medoids是在樣本集中選擇k個作為聚簇中心。即便K-Means初始過程中也可以選擇樣本點，但是在後續Step中，K-Medoids一直保持聚簇中心為樣本點，K-Means則不然。另外，K-Medoids使用的是和K-Medians一樣的曼哈頓距離，更新聚簇中心時，在簇內選擇一個使得簇內曼哈頓距離之和最小的樣本點作為聚簇中心。選擇樣本點作為聚簇中心的一個好處是聚簇中心可能是稀疏的，比如在文本聚類中，文本向量大多數是稀疏向量，係數向量計算相似度或距離速度更快。

K-Means++。初始聚簇中心的選擇是K-Means的一個重要難點，K-Means++使用了一種方法來尋找到一組比較好的初始聚簇中心。其主要思想是：首先，選擇一個樣本點作為第一個聚簇中心，然後再選擇下一個聚簇中心時要盡量離第一個已經選擇的簇中心遠一些，選擇的方法是根據一個概率分佈，這個概率分佈和樣本到聚簇中心的距離有關，比如 p sim k d^2(x, mu_1) ，即已選擇的聚簇中心越遠概率越大；對於第三個聚簇中心的選擇，需要加權考慮到第一、第二個已選擇聚簇中心的距離。重複上述過程直到選出k個聚簇中心。K-Means++可以很好地解決K-Means的初始化問題。

其他聚類方法。為瞭解決K-Means的K初始值選取的問題，可以採用層次聚類（Hierarchical Clustering）的方法；為瞭解決K-Means未考慮密度分佈以及非凸聚類的缺陷，可以考慮使用密度聚類DBSCAN演算法。此外還有很多聚類演算法，比如混合高斯聚類，這裡不再列舉。

3. Speed up K-Means with Random Approximation

其實，本文之前主要是想介紹這一部分內容的，但是考慮到這一部分比較枯燥，比較偏理論，所以在前面加了兩部分的科普內容，從簡入繁吧。

在K-Means的Assignment過程中，將樣本劃分為距離最近的一個聚簇，這個過程可以使用下面的式子來表示：

Vert x_i - mu_j Vert_2^2 = Vert x_i Vert_2^2 + Vert mu_j Vert_2^2 - 2x_i^Tmu_j

每個數據點樣本都要和聚簇中心做一次上述操作，所以一次迭代過程中計算複雜度為 O(nkd) ， d 是 指 樣本 維度。觀察上面的式子， Vert x_i Vert_2^2 的計算是可以預先完成的， Vert mu_j Vert_2^2 可以在每次迭代時預先計算好，而不必對每個樣本都計算一次。所以更新的關鍵操作就是 x_i^Tmu_j 的計算，即計算數據樣本矩陣 A 和 聚簇 中心 mu_j 的內積 A^Tmu_j 。所以我們可以想辦法在樣本矩陣上做點操作來減小計算的複雜度。

這裡說明一下，K-Means相比於其餘聚類演算法速度已經很快了，如果再加上好的初始聚簇中心，那麼收斂速度更快。但是，K-Means的應用場景經常是文本聚類、圖像聚類，數據多，維度高，所以計算量還是比較大的。

所以，如何對樣本矩陣 A 進行 變換才能使得既可以保證聚類結果不會變得很差又能提高收斂速度呢？這就是ICML2018的一篇論文《K-means clustering using random matrix sparsification》以及NIPS2010的一篇論文《Random Projection for k-means Clustering》所做的事情了。

首先介紹文章《Random Projection for k-means Clustering》提出的使用Random Projection將樣本矩陣 Ain R^{n imes d}變為矩陣 A^{}in R^{n imes d^{}}，其中 d^{} ll d 。即將樣本進行降維（Dimension Reduction）到低維空間進行聚類，可能讀者會好奇，這不就是PCA嗎？當然，PCA可以完成這麼一件事情，並且效果比較好，事實上可以這麼做。但是，Random Projection與PCA不一樣，其操作簡單，只要構建一個投影矩陣即可，而PCA降維還要做SVD，計算開銷比較大。下面Figure 2給出了演算法流程，圖來自論文截圖：

Figure 2 : Random Projection for K-Means

演算法首先選擇 t geq ck/epsilon^2 作為投影后的樣本矩陣的維度，構造一個投影矩陣 Rin R^{d imes t} ，每個元素 R_{ij} in {-frac{1}{sqrt{t}}, frac{1}{sqrt{t}}} ，概率分別為1/2。然後即可將樣本矩陣投影到維度比較小的空間中去，在新的樣本矩陣上做K-Means得到的結果會以大概率逼近在原來樣本矩陣上聚類得到的結果。至於理論上的bound這裡就不仔細介紹了，可以參考論文。

上面給出的Random Projection的方法 簡單，很容易實現，並且有很好的理論支持。下面介紹一下ICML2018的《K-means clustering using random matrix sparsification》，同樣是對樣本矩陣進行變換，這篇文章不是進行降維，而是進行稀疏化。 因為稀疏化的矩陣可以加快矩陣乘法的計算速度，所以也不失為一種有效的加速方法。文章題目裡面的Random Sparsification即是在矩陣中隨機選取一部分元素，將其變為零，達到稀疏化的目的。

文章中給出了兩種Random Sparsification的方法：第一種是uniform的，即每個元素變為零元素的概率相同，參考Figure 3；第二種是non-uniform的，為每個元素根據值的大小與矩陣中絕對值最大的元素的比值進行確定變為零的概率，簡單說就是元素絕對值越大，變為零的概率越大，參考Figure 4。對矩陣進行稀疏化之後，然後進行K-Means聚類。

Figure 3 : uniform sparsification

Figure 4 : non-uniform sparsification

在第一篇文章中，通過Random Projection得到一個投影后的矩陣widetilde{A}，這個 widetilde{A}是有 名字 的，稱之為Projection-Cost Preserving Sketches，它的定義是Figure 5（圖片來自論文《Dimensionality Reduction for k-Means Clustering and Low Rank Approximation》截圖）：

Figure 5 : Projection-Cost Preserving Sketches

同樣地，第二篇論文通過Random Sparsification尋求一個 widetilde{A} 來代替 A 。兩篇論文得到的矩陣都能實現以大概率逼近在原來樣本矩陣 A 上得到的聚類結果。

4. 實驗部分

由於兩篇論文的代碼很容易實現，所以我在Mnist測試了一下，速度的確提升了很多，並且聚類性能也較好。實驗結果和代碼這裡就不展示了，放在我的Github上，地址為：

lxcnju/approximated_kmeans?

github.com圖標
推薦閱讀：
相關文章