原理部分主要來自大牛zouxy09和trnadomeet兩個人的博客；后面的代碼詳細講解為自己精心編寫

一、概述

?????????非監督學習的一般流程是：先從一組無標簽數據中學習特征，然后用學習到的特征提取函數去提取有標簽數據特征，然后再進行分類器的訓練和分類。之前說到，一般的非監督學習算法都存在很多hyper-parameters需要調整。而，最近我們發現對于上面同樣的非監督學習流程中，用K-means聚類算法來實現特征學習，也可以達到非常好的效果，有時候還能達到state-of-the-art的效果。亮瞎了凡人之俗眼。

?????????托“bag of features?”的福，K-means其實在特征學習領域也已經略有名氣。今天我們就不要花時間迷失在其往日的光芒中了。在這里，我們只關注，如果要K-means算法在一個特征學習系統中發揮良好的性能需要考慮哪些因素。這里的特征學習系統和其他的Deep Learning算法一樣：直接從原始的輸入（像素灰度值）中學習并構建多層的分級的特征。另外，我們還分析了K-means算法與江湖中其他知名的特征學習算法的千絲萬縷的聯系（天下武功出少林，哈哈）。

?????????經典的K-means聚類算法通過最小化數據點和最近鄰中心的距離來尋找各個類中心。江湖中還有個別名，叫“矢量量化vector quantization”（這個在我的博客上也有提到）。我們可以把K-means當成是在構建一個字典D?R^nxk，通過最小化重構誤差，一個數據樣本x⁽ⁱ⁾?Rⁿ可以通過這個字典映射為一個k維的碼矢量。所以K-means實際上就是尋找D的一個過程：

?????? 這里，s⁽ⁱ⁾就是一個與輸入x⁽ⁱ⁾對應的碼矢量。D^(j)是字典D的第j列。K-means畢生的目標就是尋找滿足上面這些條件的一個字典D和每個樣本x⁽ⁱ⁾對應的碼矢量s⁽ⁱ⁾。我們一起來分析下這些條件。首先，給定字典D和碼矢量s⁽ⁱ⁾，我們需要能很好的重構原始的輸入x⁽ⁱ⁾。數學的表達是最小化x⁽ⁱ⁾和它的重構D s⁽ⁱ⁾。這個目標函數的優化需要滿足兩個約束。首先，|| s⁽ⁱ⁾||₀<=1，意味著每個碼矢量s⁽ⁱ⁾被約束為最多只有一個非零元素。所以我們尋找一個x⁽ⁱ⁾對應的新的表達，這個新的表達不僅需要更好的保留x⁽ⁱ⁾的信息，還需要盡可能的簡單。第二個約束要求字典的每列都是單位長度，防止字典中的元素或者特征變得任意大或者任意小。否則，我們就可以隨意的放縮D^(j)和對應的碼矢量，這樣一點用都木有。

?????????這個算法從精神層面與其他學習有效編碼的算法很相似，例如sparse coding：

???????? Sparse coding也是優化同樣類型的重構。但對于編碼復雜度的約束是通過在代價函數中增加一個懲罰項λ|| s⁽ⁱ⁾||₁，以限制s⁽ⁱ⁾是稀疏的。這個約束和K-means的差不多，但它允許多于一個非零值。在保證s⁽ⁱ⁾簡單的基礎上，可以更準確的描述x⁽ⁱ⁾。

?????????雖然Sparse coding比K-means性能要好，但是Sparse coding需要對每個s⁽ⁱ⁾重復的求解一個凸優化問題，當拓展到大規模數據的時候，這個優化問題是非常昂貴的。但對于K-means來說，對s⁽ⁱ⁾的優化求解就簡單很多了：

?????????這個求解是很快的，而且給定s求解D也很容易，所以我們可以通過交互的優化D和s（可以了解下我博客上的EM算法的博文）來快速的訓練一個非常大的字典。另外，K-means還有一個讓人青睞的地方是，它只有一個參數需要調整，也就是要聚類的中心的個數k。

二、數據、數據預處理

2.1數據：?

下載的CIFAR-10數據庫；其中的每個data_batch都是10000x3072大小的，即有1w個樣本圖片，每個圖片都是32*32且rgb三通道的，這里的每一行表示一個樣本。

提取40W個訓練patches：

patches = zeros(numPatches, rfSize*rfSize*3);%400000*108

for i=1:numPatches

? if (mod(i,10000) ==0) fprintf('Extracting patch: %d / %d\n', i, numPatches);end

? r= random('unid', CIFAR_DIM(1) - rfSize +1);%從1到27中采用均勻分布提取數

? c= random('unid', CIFAR_DIM(2) - rfSize +1);

? %使用mod(i-1,size(trainX,1))是因為對每個圖片樣本，提取出numPatches/size(trainX,1)=8個patch；每個size(trainX,1)=50000為一輪，提取8輪，總計40W個patches

?patch = reshape(trainX(mod(i-1,size(trainX,1))+1, :), CIFAR_DIM);%32*32*3

?patch = patch(r:r+rfSize-1,c:c+rfSize-1,:);%6*6*3

?patches(i,:) = patch(:)';%patches的每一行代表一個小樣本

end

相當于在原始圖片上隨機抽取8個6*6的rgb 圖像patches。

為了建立一個“完備complete”的字典（至少有256個類中心的字典），我們需要保證有足夠的用以訓練的patches，這樣每個聚類才會包含一定合理數量的輸入樣本。實際上，訓練一個k-means字典比其他的算法（例如sparse coding）需要的訓練樣本個數要多，因為在k-means中，一個樣本數據只會貢獻于一個聚類的中心，換句話說一個樣本只能屬于一個類，與其他類就毫無瓜葛了，該樣本的信息全部傾注給了它的歸宿（對應類的中心）。（以上分析來自zouxy09博客）。所以本文作者用40W個隨機patches，來訓練108維的聚類特征。

2.2、預處理?Pre-processing

?????????給定多張圖片構成的一個矩陣（其中每張圖片看成是一個向量，多張圖片就可以看做是一個矩陣了）。要對這個矩陣進行whitening操作，而在這之前是需要均值化的。在以前的實驗中，有時候是對每一張圖片內部做均值，也就是說均值是針對每張圖片的所有維度，而有的時候是針對矩陣中圖片的每一維做均值操作，那么是不是有矛盾呢？其實并不矛盾，主要是這兩種均值化的目的不同。如果是算該均值的協方差矩陣，或者將一些訓練樣本輸入到分類器訓練前，則應該對每一維采取均值化操作（因為協方差均值是描述每個維度之間的關系）。如果是為了增強每張圖片亮度的對比度，比如說在進行whitening操作前，則需要對圖片的內部進行均值化（此時一般還會執行除以該圖像內部的標準差操作）。

在訓練之前，我們需要對所有的訓練樣本patches先進行亮度和對比度的歸一化。具體做法是，對每個樣本，我們減去灰度的均值和除以標準差。另外，在除以標準差的時候，為了避免分母為0和壓制噪聲，我們給標準差增加一個小的常數。對于[0, 255]范圍的灰度圖，給方差加10一般是ok的：

樣本歸一化代碼：

patches = bsxfun(@rdivide, bsxfun(@minus,patches, mean(patches,2)), sqrt(var(patches,[],2)+10));

?

?????????對訓練數據進行歸一化后，我們就可以在上面運行k-means算法以獲得相應的聚類中心了（字典D的每一列），可視化在圖a中，可以看到，k-means趨向于學習到低頻的類邊緣的聚類中心。但不幸的是，這樣的特征不會有比較好的識別效果，因為鄰域像素的相關性（圖像的低頻變化）會很強。這樣K-means就會產生很多高度相關的聚類中心。而不是分散開聚類中心以更均勻的展開訓練數據。所以我們需要先用白化來去除數據的相關性，以驅使K-means在正交方向上分配更多的聚類中心。

??????圖a是由k-means從沒有經過白化處理的自然圖像中學習到的聚類中心。圖b展示了沒有和有白化的效果。左邊是沒有經過白化的，因為數據存在相關性，所以聚類中心會跑偏。右邊是經過白化后的，可以看到聚類中心是更加正交的，這樣學習到的特征才會不一樣。圖c展示的是從經過白化后的圖像patches學習到的聚類中心。

?????????實現whitening白化一個比較簡單的方法是ZCA白化（可以參考UFLDL）。我們先對數據點x的協方差矩陣進行特征值分解cov(x)=VDV^T。然后白化后的點可以表示為：

??????? ?_zca是一個很小的常數。對于對比度歸一化后的數據，對16x16的patch，可以設?_zca=0.01，對8x8的patch，可以設?_zca=0.1?。需要注意的一點是，這個常數不能太小，如果太小會加強數據的高頻噪聲，會使特征學習更加困難。另外，因為數據的旋轉對K-means沒有影響，所以可以使用其他的白化變換方法，例如PCA白化（與ZCA不同只在于其旋轉了一個角度）。

在白化后的數據中運行k-means可以得到清晰的邊緣特征。這些特征和sparse coding啊，ICA啊等等方法學到的初級特征差不多。如圖c所示。另外，對于新的數據，可能需要調整歸一化的參數?和白化的參數?_zca，以滿足好的效果（例如，圖像patches具有高對比度，低噪聲和少低頻波動）。

數據白化代碼：

C = cov(patches);%計算patches的協方差矩陣

? M= mean(patches);

?[V,D] = eig(C);

? P= V * diag(sqrt(1./(diag(D) +0.1))) * V';%P是ZCA Whitening矩陣

? %對數據矩陣白化前，應保證每一維的均值為0

?patches = bsxfun(@minus, patches, M) * P;%注意patches的行列表示的意義不同時，白化矩陣的位置也是不同的。

?

三、用kmeans算法獲取聚類中心

詳見：k-means聚類原理代碼分析

四、特征提取

?????????我們用K-means聚類算法來從輸入數據中學習K個聚類中心c^(k)。當學習到這K個聚類中心后，我們可以有兩種特征映射f的方法。第一種是標準的1-of-K，屬于硬分配編碼：

??????這個f_k(x)表示樣本x的特征向量f的第k個元素，也就是特征分量。為什么叫硬分配呢，因為一個樣本只能屬于某類。也就是說每個樣本x對應的特征向量里面只有一個1，其他的全是0。K個類有K個中心，樣本x與哪個中心歐式距離最近，其對應的位就是1，其他位全是0，屬于稀疏編碼的極端情況，高稀疏度啊。

第二種是采用非線性映射。屬于軟編碼。

??????這里z_k=||x-c^(k)||₂，u(z)是向量z的均值。也就是先計算出對應樣本與k個類中心點的平均距離d，然后如果那些樣本與類別中心點的距離大于d的話都設置為0，小于d的則用d與該距離之間的差來表示。這樣基本能夠保證一半以上的特征都變成0了，也是具有稀疏性的，且考慮了更多那些距類別中心距離比較近的值。為什么叫軟分配呢。軟表示這個樣本以一定的概率屬于這個類，它還以一定的概率屬于其他的類，有點模糊聚類的感覺。而硬則表示這個樣本與哪個類中心距離最近，這個樣本就只屬于這個類。與其他類沒有任何關系。

特征提取代碼：

1 ，原始圖像數據拆分變換：

im2col函數：

該函數是將一個大矩陣按照小矩陣取出來，并把取出的小矩陣展成列向量。比如說B =im2col(A,[m n],block_type):就是把A按照m*n的小矩陣塊取出，取出后按照列的方式重新排列成向量,然后多個列向量組成一個矩陣。而參數block_type表示的是取出小矩形框的方式，有兩種值可以取，分別為’distinct’和’sliding’。Distinct方式是指在取出的各小矩形在原矩陣中是沒有重疊的，元素不足的補0。而sliding是每次移動一個元素，即各小矩形之間有元素重疊，但此時沒有補0元素的說法。如果該參數不給出，則默認的為’sliding’模式。

代碼：

patches = [ im2col(reshape(X(i,1:1024),CIFAR_DIM(1:2)), [rfSize rfSize]) ;?

? ? ? ? ? ? ? ? ? ?im2col(reshape(X(i,1025:2048),CIFAR_DIM(1:2)), [rfSize rfSize]) ;??????????????

? ? ? ? ? ? ? ? ? ?im2col(reshape(X(i,2049:end),CIFAR_DIM(1:2)), [rfSize rfSize]) ]';

最后，通過類似“卷積拆分”變換將原始圖像變為一個729*108的矩陣，每行為原始圖像提取的6*6rgb? patch的向量表達；卷積塊總共移動729次。

2 ，數據的標準化，白化處理

標準化

patches =

bsxfun(@rdivide,bsxfun(@minus, patches, mean(patches,2)), sqrt(var(patches,[],2)+10));

白化

patches = bsxfun(@minus, patches, M) * P;

?

3 ，特征提取

此處采用是采用非線性映射，屬于軟編碼。

??????這里z_k=||x-c^(k)||₂，u(z)是向量z的均值。也就是先計算出對應樣本與k個類中心點的平均距離d，然后如果那些樣本與類別中心點的距離大于d的話都設置為0，小于d的則用d與該距離之間的差來表示。這樣基本能夠保證一半以上的特征都變成0了，也是具有稀疏性的，且考慮了更多那些距類別中心距離比較近的值。為什么叫軟分配呢。軟表示這個樣本以一定的概率屬于這個類，它還以一定的概率屬于其他的類，有點模糊聚類的感覺。而硬則表示這個樣本與哪個類中心距離最近，這個樣本就只屬于這個類。與其他類沒有任何關系。

計算樣本與k個中心的距離

xx = sum(patches.^2, 2);%729列，的列向量

cc = sum(centroids.^2, 2)';%1600列，的列向量

xc = patches * centroids';%729*1600的矩陣，

z = sqrt( bsxfun(@plus, cc, bsxfun(@minus,xx, 2*xc)) ); % distances = xx^2+cc^2-2*xx*cc;

z為729*1600的矩陣，z矩陣的每行為patches中一行（patches每行為一個rgb卷積塊6*6*3=108）與1600個聚類中心的距離。因為一個原始圖片有729個卷積塊，每個卷積塊都與1600個卷積中心計算距離，所以有729*1600個距離。

?

軟編碼：

其實就是也就是先計算出對應樣本與k個類中心點的平均距離d，然后如果那些樣本與類別中心點的距離大于d的話都設置為0，小于d的則用d與該距離之間的差來表示。

這樣基本能夠保證一半以上的特征都變成0了，也是具有稀疏性的，且考慮了更多那些距類別中心距離比較近的值。

mu = mean(z, 2); % 每個卷積塊與1600個中心的平均距離

patches = max(bsxfun(@minus, mu, z), 0);%patches中每一行保存的是：小樣本與這1600個類別中心距離的平均值減掉與每個類別中心的距離，限定最小距離為0；

a=[ -1 2 3]; max(a,0)=[0 2 3];%a中的元素與0比較；若小于0，取值為0，大于0取值本身

把patches從729*1600，變形為27*27*1600；上面的軟編碼實現了類似用6*6的卷積塊，卷積圖像的操作，1600個中心，相當于1600個特征提取卷積塊。如圖

Patches的每列為：

根據卷積運算圖像原則，一個卷積塊在原圖像上依次滑動，也就是一個特征卷積提取整個圖像；按照這個原理，patches的每列應該是每個一個聚類中心點（“卷積塊”），與原圖像每個卷積模塊求距離后，再進行軟編碼。不應該是上面原圖像的每個卷積塊，與1600個中心點求距離在進行軟計算，這里有些不明白？

?

4pool操作

prows = CIFAR_DIM(1)-rfSize+1;

pcols = CIFAR_DIM(2)-rfSize+1;

patches = reshape(patches, prows, pcols,numCentroids);

pool代碼：

halfr = round(prows/2);

halfc = round(pcols/2);

q1 = sum(sum(patches(1:halfr, 1:halfc, :), 1),2);%區域內求和，是個列向量，1600*1

q2 = sum(sum(patches(halfr+1:end,1:halfc, :), 1),2);

q3 = sum(sum(patches(1:halfr, halfc+1:end, :),1),2);

q4 = sum(sum(patches(halfr+1:end, halfc+1:end, :),1),2);

% concatenate into feature vector

XC(i,:) = [q1(:);q2(:);q3(:);q4(:)]';%類似于pooling操作

這么pool是為了數據降維，但是這是什么pooling方式，為什么用這種方式有些不明白。總之最后向量化后，就提取完特征了。

參考文獻：

Deep Learning論文筆記之（三）單層非監督學習網絡分析

Deep learning：二十五(Kmeans單層網絡識別性能)

總結

以上是生活随笔為你收集整理的经典的K-means聚类算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。