K-Means聚類演算法（二）：演算法實現及其優化

01-26

（最近在車間幹活的時候把手砸傷了，所以打字還是有點不便，大家原諒我更新的慢，加上贊比較少，心情比較低落TAT）

首先介紹一下題圖，這個是個小蘿莉的照片（如下），每一個點都可以視作是一個三維向量（點？）（RGB三通道圖片），那麼，使用K-Means演算法對這些點進行聚類，我們就很容易得到幾個中心點和幾類，把同一類的數據點（像素點）用中心點表示就可以得到壓縮後的圖片，以上分別是把3通道×每通道8bit=24bit的像素點壓縮為1bit，2bit，3bit和4bit的效果，可以看到，雖然信息損失了6倍，但是對畫質的影響沒有想像的大。

忘了說了，我把圖像壓縮的代碼傳上去了，但是圖片沒辦法傳送，因為本人的GitHub抽風了（準確的說是網路抽風）

DataScienceNote/KMeans_Img.m at master · TsingJyujing/DataScienceNote · GitHub

請準備好loli.jpg文件，並且配合k_means.m一起食用~

k_means.m代碼放在了GitHub上，大家自己下載著玩:P

DataScienceNote/k_means.m at master · TsingJyujing/DataScienceNote · GitHub

其次修正一下上次的說法：並不只是歐氏空間可以用這個方法，只要是對於提供了以下幾種運算元且滿足以下性質的空間都能用：

首先有個計算任意兩點距離（或者叫不相似度）的算符記作 $<vec{a},vec{b}>$

且有交換性：

$<vec{a},vec{b}>=<vec{b},vec{a}>$

其次有求任意一堆點的平均值的演算法： $vec{mu}=Mean(vec{{a}_{1}},..., vec{{a}_{n}})$

求出以後使得：

$sum_{i=1}^{n}{(vec{mu}-vec{{a}_{i}})}^{2}$

就可以了。

如果這一段有錯誤，求指出哈~

這一篇依舊不會涉及數學原理，我在糾結中，因為要從頭說清楚EM演算法是一件很費工夫的事情：首先是含有隱變數(Latent Variables)的模型去估計實際的數據，然後從概率論出發，到極大似然估計確定參數和E（Expectation）步和M（Maximization）步。以上是我清楚的，還有我到現在有雲里霧裡的收斂性證明。

並不是說費工夫就不寫了，逃了，而是用最少的公式和步驟說清楚和K-Means有關的部分。

下一章會比較高能，說是高能，其實也就是求個偏導數，極大似然估計啥的而已……而已……

Part Ⅰ 編程實戰的時候那些坑：

實戰的時候是會遇到很多坑的，首先說第一個：

一、如何初始化中心點

事實上，初始化的不好，是會震蕩的（當然，有了後面的解決震蕩的辦法要好得多）！

初始化的時候，不妨先估計一下數據的中心在哪裡，數據的尺度（叫Scale合適么）有多大，然後在生成數據的時候以中心點為中心，以尺度為因子乘以服從標準正態分布的數據。

大致就是這樣的一段代碼：

group_center = mean(data_set);group_range = range(data_set);centers = (randn(K,dim).*repmat(group_range,K,1)./3+repmat(group_center,K,1));

至於為什麼要除以3，是因為這裡 ${sigma}=1$ ，一般數據不會越過 $±3{sigma}$ 。

這裡根據每一個維度分別做了Scale，使得其能充填在數據中心點附近，然而又不太遠。

二、如何計算沒有「爭取」到任意一點的中心點？

有的時候某些中心點太強勢，直接把點搶沒了：主要出現在點少類多的情況下。

這個時候，就需要忽略這些中心點。但是中心點不調整也不好，有可能就在那個地方註定孤獨一生了。

這個時候我們把這個點初始化到中心點附近的地方：

centers(i,:) = (randn(1,dim).*group_range./3+group_center);

讓他重新來過。

這一點在有權重參與計算的時候特別重要麼，因為沒有爭取到點，所以這個時候權重會是0，就出現除零錯誤了。

三、我怎麼知道有沒有收斂呢？

這個在書中已經說了，使用代價函數：

$widetilde{J}=sum_{i=1}^{C}{sum_{j=1}^{N}{{r}_{ij} imes{ u({x}_{j},{mu}_{i})}}}$

其中：

$u({x}_{j},{mu}_{i})={||{x}_{j}-{mu}_{i}||}^{2}$

${r}_{ij}$ 的解釋：如果第j個數據點屬於第i類，那就記作1，否則記作0的一個 $N imes C$ 大小的矩陣。代價函數的差分值小於一定數值的時候（N次越不過最小值點）即可認為收斂了。

（以下一段算是題外廢話）

這裡需要補充另外兩種評價聚類演算法的方法：F-Measure和信息熵

但是使用這兩個代價函數以後就不能用原來提出的迭代演算法了，因為原來的EM演算法是針對上面的 $widetilde{J}$ 那個代價函數的。

關於這兩個代價函數和如何推導會在下一篇講清楚。

四、代價函數不收斂，怎麼破？

其實有一些論文和偏理論的書講了代價函數的收斂性和初值敏感性的證明，然而不管是本文還是下一篇文章我都不準備討論這些證明，主要原因是：我！不！太！懂！……（逃

（如果有希望知道這些但是不知道哪裡找資料的大犇，我推薦你《統計學習方法》中關於EM演算法的收斂性討論的部分，只不過該書沒有說K-Means演算法，直說了硬幣遊戲模型和Gaussian混合模型，真正找EM演算法在K-Means上的推導在Pattern Recognition and Machine Learning上，不過沒有求偏導和導出求中心點的詳細過程，直說了句we can easily solve，翻譯過來就是「易證」、「易知」、「易求」或者「我們把這個簡單有趣的證明留給讀者」之類的鳥話）

首先說一下什麼時候容易發生震蕩：在數據點個數比較少而且比較稀疏的時候容易發生這種事情，發生的原因大約有兩種常見的：

1、陷入某個環里，然後開始震蕩，你就會看見某個中心點繞了一圈一圈的……低頻振蕩。

2、兩個點交換來交換去，每次交換不改變J的值就收斂了，如果交換以後不幸影響了其它的點，就出現了高頻振蕩（多一句嘴，為什麼稱之為高頻呢？因為達到了奈奎斯特採樣頻率最大值啊）。

這個時候給出一種簡單的解決方案：阻尼。

簡而言之，就是更新自己位置的時候考慮一下原來的位置，一般阻尼比（在0~1之間取值）決定收斂速度，收斂的慢了也就不容易震蕩，也就越容易陷入局部極小值，也就是說，不震蕩的情況下我們應該把阻尼比儘可能取小一點

$vec{{C}}^{upd}=vec{{C}}^{new} imes(1-xi)+vec{{C}}^{old} imesxi$

其中：

$vec{{C}}^{upd}$ 是最後中心點的取值， $vec{{C}}^{new}$ 是當前集合的中心點， $vec{{C}}^{old}$ 是原來的中心點坐標。

代碼中的體現就是：

dumper = 0.1;...centers = (1-dumper).*centers + dumper.*lst_cnt;

五、如何給某些點設置權重

首先明確一個問題：為什麼要設置權重？

設置權重一般發生在兩種情況下：

1，數據集中有大量的重複點，而且數據量比較大計算非常的燒CPU。

2，數據本身有輕重緩急之分，例如只買一兩條褲子還扭扭捏捏挑挑選選的我這樣的小屌絲客戶和喬布斯那種為了自己喜歡的黑色套頭衫讓停產的生產線再運行的情懷客戶的權重肯定是不一樣的（請一口氣讀完）。

可以看看第一篇文章的標題圖（封面）就是帶權聚類的一個結果：

那麼問題來了：如何設置權重？

設置權重其實不困難，首先要明確代價函數是什麼：是點距離的和啊！那麼更新策略是什麼？是找中心點啊！那麼現在我們就需要引入加權平均值的演算法：

${Mean}_{weight}({vec{{v}_{i}}},vec{w})=sum_{i=1}^{n}{{w}_{i} imesvec{{v}_{i}}}$

每次移到這個地方就可以了，其中 $vec{w}$ 滿足 ${L}^{1}$ 範數是1，這是文藝說法。

普通說法就是 $sum{{w}_{i}}=1$

如果不滿足也懶得歸一化可以：

${Mean}_{weight}({vec{{v}_{i}}},vec{w})=frac{1}{sum{{w}_{i}}}sum_{i=1}^{n}{{w}_{i} imesvec{{v}_{i}}}$

這就是我程序中使用的公式。

具體在代碼中的表現就是：

centers(i,:) = sum(data_set(cls_idx,:).*repmat(weight(cls_idx),1,dim))... ./sum(weight(cls_idx));

Part Ⅱ Matlab編程技巧：

一、如何愉快的求取 ${L}^{2}$ 範數/距離？

首先，我們一起確認一下 ${L}^{2}$ 範數的公式：

對於任意一個向量： $vec{v}=({v}_{1},...,{v}_{d})$

其 ${L}^{2}$ 範數是：

${left| x ight|} ^{2}= sum_{i=1}^{D}{{{v}_{i}}^{2}}$

隨後，我們確認一下求範數的幾個步驟：

對於已知的一點pnt和行向量組data_set，首先確認行向量組的向量個數：

[N,~]=size(data_set )

隨後用repmat把pnt複製N遍，變成和data_set同樣大小的一個矩陣。

然後那麼一減啊，再求個平方，就可以開始求和了（就是那個Σ），但是是橫向拍扁求和，因為我們是行向量啊。

所以給sum函數第二個參數：2

最後開平方即可：

dist_list = sqrt(sum((data_set - repmat(pnt,N,1)).^2,2))

二、如何愉快的計算某一個點屬於哪一類：

我們剛才求出了很多的列（dist_list ），這些就是到每個點的距離的列表，這裡一共有k個點，所以做成一個k*N的矩陣：

for i = 1:K %求到每一個中心的距離向量，且不用張量求更加節約空間 dist_vec(:,i) = sqrt(sum((data_set - repmat(centers(i,:),N,1)).^2,2)); end

這個時候，到第一個點距離最短那麼就是1，第二個就是2，如果還用For循環一個個判斷你就out了，min函數的第二個返回值是min值所在的位置：

[~,cls_vec] = min(dist_vec,[],2);

只要你樂意，其實還可以生成{r}_{n,k}矩陣（PRML書里用的，偏理論）

rnk = zeros(k,N);rnk = rnk((0N-1)).*k+cls_vec)";

其實沒必要。

三、用什麼做返回值

這是個仁者見仁丹、智者見智障的問題，Matlab中用剛才計算的cls_vec作為返回值，其實還可以設計一些更加實用的返回值（看你怎麼用了）

not_empty_cls = unique(cls_vec);k_means_result = cell(length(not_empty_cls),1);for i = 1:length(not_empty_cls) cls = not_empty_cls(i); sub_res = cell(3,1); cls_idx = find(cls_vec==cls); sub_res{1} = centers(cls,:); sub_res{3} = cls_idx; sub_res{2} = data_set(cls_idx,:); k_means_result{i} = sub_res;end

首先看一下有那幾類不空的，好、假設有N類。

然後生成一個 ${N} imes {1}$ 的Cell向量，每一個Cell裡面放3個Cell：

1，中心點，2，屬於這一類的原始數據，3，剛才那些數據在輸入數據里的Index。

這樣的好處是：信息比較全面，調用方便。

缺點是不夠直觀么，層次複雜。

你只要開心可以兩種都返回，沒人攔著你（就好像我代碼里示範的那樣）
推薦閱讀：

※對比傳統K-Means等聚類演算法，LDA主題模型在文本聚類上有何優缺點？
※海量數據的聚類通常如何做？
※哪種聚類演算法可以不需要指定聚類的個數，而且可以生成聚類的規則？
※如何根據每個策略的 daily return 對不同策略進行最為有效的分類？

TAG:机器学习 | 聚类算法 | 图片压缩 |