?一. 前言

?

?

?

本來想寫關(guān)于聚類系列算法的介紹，但是聚類系列的其它幾個(gè)算法原理比較簡單，網(wǎng)上有大量的教程可以查閱。這里主要是介紹一下譜聚類算法，做一個(gè)學(xué)習(xí)筆記，同時(shí)也希望對想要了解該算法的朋友有一個(gè)幫助。關(guān)于聚類的其他系列算法，這里推薦一個(gè)寫的很不錯(cuò)的博客。

譜聚類在最近幾年變得受歡迎起來，主要原因就是它實(shí)現(xiàn)簡單，聚類效果經(jīng)常優(yōu)于傳統(tǒng)的聚類算法（如K-Means算法）。剛開始學(xué)習(xí)譜聚類的時(shí)候，給人的感覺就是這個(gè)算法看上去很難，但是當(dāng)真正的深入了解這個(gè)算法的時(shí)候，其實(shí)它的原理并不難，但是理解該算法還是需要一定的數(shù)學(xué)基礎(chǔ)的。如果掌握了譜聚類算法，會(huì)對矩陣分析，圖論和降維中的主成分分析等有更加深入的理解。

本文首先會(huì)簡單介紹一下譜聚類的簡單過程，然后再一步步的分解這些過程中的細(xì)節(jié)以及再這些過程中譜聚類是如何實(shí)現(xiàn)的，接著總結(jié)一下譜聚類的幾個(gè)算法，再接著介紹譜聚類算法是如何用圖論知識(shí)來描述，最后對本文做一個(gè)總結(jié)。下面我們就來介紹一下譜聚類的基本原理。

二. 譜聚類的基本原理介紹

2.1 譜和譜聚類

2.1.1 譜

方陣作為線性算子，它的所有特征值的全體統(tǒng)稱為方陣的譜。方陣的譜半徑為最大的特征值。矩陣A的譜半徑是矩陣

的最大特征值。

2.1.2 譜聚類

譜聚類是一種基于圖論的聚類方法，通過對樣本數(shù)據(jù)的拉普拉斯矩陣的特征向量進(jìn)行聚類，從而達(dá)到對樣本數(shù)據(jù)聚類的母的。譜聚類可以理解為將高維空間的數(shù)據(jù)映射到低維，然后在低維空間用其它聚類算法（如KMeans）進(jìn)行聚類。

2.2 譜聚類算法簡單描述

輸入：n個(gè)樣本點(diǎn)

和聚類簇的數(shù)目k；

輸出：聚類簇

（1）使用下面公式計(jì)算

的相似度矩陣W；

W為

組成的相似度矩陣。

（2）使用下面公式計(jì)算度矩陣D；

，即相似度矩陣W的每一行元素之和

D為

組成的

對角矩陣。

（3）計(jì)算拉普拉斯矩陣

；

（4）計(jì)算L的特征值，將特征值從小到大排序，取前k個(gè)特征值，并計(jì)算前k個(gè)特征值的特征向量

；

（5）將上面的k個(gè)列向量組成矩陣

，

；

（6）令

是

的第

行的向量，其中

；

（7）使用k-means算法將新樣本點(diǎn)

聚類成簇

；

（8）輸出簇

，其中，

.

上面就是未標(biāo)準(zhǔn)化的譜聚類算法的描述。也就是先根據(jù)樣本點(diǎn)計(jì)算相似度矩陣，然后計(jì)算度矩陣和拉普拉斯矩陣，接著計(jì)算拉普拉斯矩陣前k個(gè)特征值對應(yīng)的特征向量，最后將這k個(gè)特征值對應(yīng)的特征向量組成

的矩陣U，U的每一行成為一個(gè)新生成的樣本點(diǎn)，對這些新生成的樣本點(diǎn)進(jìn)行k-means聚類，聚成k類，最后輸出聚類的結(jié)果。這就是譜聚類算法的基本思想。相比較PCA降維中取前k大的特征值對應(yīng)的特征向量，這里取得是前k小的特征值對應(yīng)的特征向量。但是上述的譜聚類算法并不是最優(yōu)的，接下來我們一步一步的分解上面的步驟，總結(jié)一下在此基礎(chǔ)上進(jìn)行優(yōu)化的譜聚類的版本。

2.3 譜聚類算法中的重要屬性

2.3.1 相似度矩陣介紹

相似度矩陣就是樣本點(diǎn)中的任意兩個(gè)點(diǎn)之間的距離度量，在聚類算法中可以表示為距離近的點(diǎn)它們之間的相似度比較高，而距離較遠(yuǎn)的點(diǎn)它們的相似度比較低，甚至可以忽略。這里用三種方式表示相似度矩陣：一是

-近鄰法（

-neighborhood graph），二是k近鄰法（k-nearest nerghbor graph），三是全連接法（fully connected graph）。下面我們來介紹這三種方法。

（1）

-neighborhood graph：

，表示樣本點(diǎn)中任意兩點(diǎn)之間的歐式距離

用此方法構(gòu)造的相似度矩陣表示如下：

該相似度矩陣由于距離近的點(diǎn)的距離表示為

，距離遠(yuǎn)的點(diǎn)距離表示為0，矩陣種沒有攜帶關(guān)于數(shù)據(jù)集的太多的信息，所以該方法一般很少使用，在sklearn中也沒有使用該方法。

（2）k-nearest nerghbor graph：

由于每個(gè)樣本點(diǎn)的k個(gè)近鄰可能不是完全相同的，所以用此方法構(gòu)造的相似度矩陣并不是對稱的。因此，這里使用兩種方式表示對稱的knn相似度矩陣，第一種方式是如果

在

的k個(gè)領(lǐng)域中或者

在

的k個(gè)領(lǐng)域中，則

為

與

之間的距離，否則為

；第二種方式是如果

在

的k個(gè)領(lǐng)域中并且

在

的k個(gè)領(lǐng)域中，則

為

與

之間的距離，否則為

。很顯然第二種方式比第一種方式生成的相似度矩陣要稀疏。這兩種方式用公式表達(dá)如下：

第一種方式：

第二種方式：

（3）fully connected graph:

該方法就是在算法描述中的高斯相似度方法，公式如下：

該方法也是最常用的方法，在sklearn中默認(rèn)的也是該方法，表示任意兩個(gè)樣本點(diǎn)都有相似度，但是距離較遠(yuǎn)的樣本點(diǎn)之間相似度較低，甚至可以忽略。這里面的參數(shù)

控制著樣本點(diǎn)的鄰域?qū)挾?，?/p>

越大表示樣本點(diǎn)與距離較遠(yuǎn)的樣本點(diǎn)的相似度越大，反之亦然。

2.3.2 拉普拉斯矩陣介紹

對于譜聚類來說最重要的工具就是拉普拉斯矩陣了，下面我們來介紹拉普拉斯矩陣的三種表示方法。

（1）未標(biāo)準(zhǔn)化的拉普拉斯矩陣：

未標(biāo)準(zhǔn)化的拉普拉斯矩陣定義如下：

其中W是上節(jié)所說的相似度矩陣，D是度矩陣，在算法描述中有介紹。很顯然，W與D都是對稱矩陣。

未標(biāo)準(zhǔn)化的拉普拉斯矩陣L滿足下面幾個(gè)性質(zhì)：

（a）對任意一個(gè)向量

都有：

證明如下：

（b）L是對稱的和半正定的，證明如下：

因?yàn)?/p>

，所以

，所以為半正定矩陣。由于W和D都是對稱矩陣，所以L為對稱矩陣。

（c）L最小的特征值為0，且特征值0所對應(yīng)的特征向量為全1向量，證明如下：

令

表示

的全1向量，則

由D和W的定義可以得出上式。

（d）L有n個(gè)非負(fù)的實(shí)數(shù)特征值：

（2）標(biāo)準(zhǔn)化拉普拉斯矩陣

標(biāo)準(zhǔn)化拉普拉斯矩陣有兩種表示方法，一是基于隨機(jī)游走（Random Walk）的標(biāo)準(zhǔn)化拉普拉斯矩陣

和對稱標(biāo)準(zhǔn)化拉普拉斯矩陣

，定義如下：

標(biāo)準(zhǔn)化的拉普拉斯矩陣滿足如下性質(zhì)：

（a）對任意一個(gè)向量

都有：

（b）當(dāng)且僅當(dāng)

是

的特征值，對應(yīng)的特征向量為

時(shí)，則

是

特征值，對應(yīng)的特征向量為u；

（c）當(dāng)且僅當(dāng)

時(shí)，

是

的特征值，對應(yīng)的特征向量為u；

（d）0是

的特征值，對應(yīng)的特征向量為

，

為

的全1向量；0也是

的特征值，對應(yīng)的特征向量為

；

（e）

和

是半正定矩陣并且有非負(fù)實(shí)數(shù)特征值：

.

關(guān)于各個(gè)版本的譜聚類算法的不同之處，就是在于相似度矩陣的計(jì)算方式不同和拉普拉斯矩陣的表示方法不同，其它步驟基本相同。下面就來介紹關(guān)于譜聚類的兩個(gè)比較流行的標(biāo)準(zhǔn)化算法。

2.4 標(biāo)準(zhǔn)化譜聚類算法介紹

2.4.1 隨機(jī)游走拉普拉斯矩陣的譜聚類算法描述

輸入：n個(gè)樣本點(diǎn)

和聚類簇的數(shù)目k；

輸出：聚類簇

（1）計(jì)算

的相似度矩陣W；

（2）計(jì)算度矩陣D；

（3）計(jì)算拉普拉斯矩陣

；

（4）計(jì)算

的特征值，將特征值從小到大排序，取前k個(gè)特征值，并計(jì)算前k個(gè)特征值的特征向量

；

（5）將上面的k個(gè)列向量組成矩陣

，

；

（6）令

是

的第

行的向量，其中

；

（7）使用k-means算法將新樣本點(diǎn)

聚類成簇

；

（8）輸出簇

，其中，

.

2.4.2 對稱拉普拉斯矩陣的譜聚類算法描述

輸入：n個(gè)樣本點(diǎn)

和聚類簇的數(shù)目k；

輸出：聚類簇

（1）計(jì)算

的相似度矩陣W；

（2）計(jì)算度矩陣D；

（3）計(jì)算拉普拉斯矩陣

；

（4）計(jì)算

的特征值，將特征值從小到大排序，取前k個(gè)特征值，并計(jì)算前k個(gè)特征值的特征向量

；

（5）將上面的k個(gè)列向量組成矩陣

，

；

（6）令

是

的第

行的向量，其中

；

（7）對于

，將

依次單位化，使得

；

（8）使用k-means算法將新樣本點(diǎn)

聚類成簇

；

（9）輸出簇

，其中，

.

上面兩個(gè)標(biāo)準(zhǔn)化拉普拉斯算法加上未標(biāo)準(zhǔn)化拉普拉斯算法這三個(gè)算法中，主要用到的技巧是將原始樣本點(diǎn)

轉(zhuǎn)化為新的樣本點(diǎn)

，然后再對新樣本點(diǎn)使用其它的聚類算法進(jìn)行聚類，在這里最后一步用到的聚類算法不一定非要是KMeans算法，也可以是其它的聚類算法，具體根據(jù)實(shí)際情況而定。在sklearn中默認(rèn)是使用KMeans算法，但是由于KMeans聚類對初始聚類中心的選擇比較敏感，從而導(dǎo)致KMeans算法不穩(wěn)定，進(jìn)而導(dǎo)致譜聚類算法不穩(wěn)定，所以在sklearn中有另外一個(gè)可選項(xiàng)是'discretize'，該算法對初始聚類中心的選擇不敏感。

三. 用切圖的觀點(diǎn)來解釋譜聚類

聚類算法給我們的直觀上的感覺就是根據(jù)樣本點(diǎn)的相似性將他們劃分成不同的組，使得在相同組內(nèi)的數(shù)據(jù)點(diǎn)是相似的，不同組之間的數(shù)據(jù)點(diǎn)是不相似的。對于給定的樣本點(diǎn)計(jì)算相似度，形成相似度圖，因而譜聚類的問題可以被重新描述如下：我們想要找到圖形的一個(gè)分區(qū)，使得不同分區(qū)之間的邊具有非常低的權(quán)重（這意味著不同分區(qū)中的點(diǎn)彼此不相似）并且分區(qū)內(nèi)的邊具有高權(quán)重（這意味著其中的點(diǎn)彼此相似）。在這個(gè)小節(jié)我們將討論如何推導(dǎo)譜聚類為近似的圖分區(qū)的問題。

3.1 最小切（mincut）

對于無向圖G，我們的目標(biāo)是將圖

切成互相沒有連接的子圖，每個(gè)子圖點(diǎn)的集合為

，其中

且

。

對于任意兩個(gè)子圖點(diǎn)的集合

，

，我們定義A和B之間的權(quán)重切圖為：

對于k個(gè)子圖集合

，我們定義切圖cut為：

其中，

為A的補(bǔ)集。

我們可以想到的切圖方法就是最小化上式，也就是使各個(gè)組之間的權(quán)重盡可能的小，但是在許多情況下mincut只是簡單的將圖中的一個(gè)定點(diǎn)與其余的頂點(diǎn)分開，在并不是我們想要的結(jié)果，合理的切分結(jié)果應(yīng)該是組內(nèi)的樣本點(diǎn)盡可能的多。所以mincut在實(shí)際中并不常用。下面我們介紹另外兩個(gè)切圖方式RatioCut和Ncut。

3.2 RatioCut切圖

在RatioCut切圖中，不僅要考慮使不同組之間的權(quán)重最小化，也考慮了使每個(gè)組中的樣本點(diǎn)盡量多。

定義

為子集A中的頂點(diǎn)個(gè)數(shù)，RatioCut的表達(dá)式如下：

將圖中頂點(diǎn)V分成k個(gè)分區(qū)

，我們定義指示向量為：

，其中：

我們令

為包含k個(gè)指示向量作為列向量的矩陣，注意到H的列向量彼此正交，即

，然后我們計(jì)算下式：

，證明如下：

結(jié)合上面可以得到：

其中Tr(A)表示矩陣A的跡，也就是矩陣A的對角線之和。

所以我們此時(shí)的目標(biāo)函數(shù)成為：

，約束條件為：

注意到我們H矩陣?yán)锩娴拿恳粋€(gè)指示向量都是n維的，向量中每個(gè)變量的取值為0或者

，就有

種取值，有k個(gè)子圖的話就有k個(gè)指示向量，共有

種H，因此找到滿足上面優(yōu)化目標(biāo)的H是一個(gè)NP難的問題。那么是不是就沒有辦法了呢？

注意觀察

中每一個(gè)優(yōu)化子目標(biāo)

，其中

是單位正交基，

為對稱矩陣。在PCA中，我們的目標(biāo)是找到協(xié)方差矩陣（對應(yīng)此處的拉普拉斯矩陣L）的最大的特征值，而在我們的譜聚類中，我們的目標(biāo)是找到目標(biāo)的最小的特征值，得到對應(yīng)的特征向量，此時(shí)對應(yīng)二分切圖效果最佳。也就是說，我們這里要用到維度規(guī)約的思想來近似去解決這個(gè)NP難的問題。

對于

，我們的目標(biāo)是找到最小的L的特征值，而對于

，則我們的目標(biāo)就是找到k個(gè)最小的特征值，一般來說，k遠(yuǎn)遠(yuǎn)小于n，也就是說，此時(shí)我們進(jìn)行了維度規(guī)約，將維度從n降到了k，從而近似可以解決這個(gè)NP難的問題。

通過找到L的最小的k個(gè)特征值，可以得到對應(yīng)的k個(gè)特征向量，這k個(gè)特征向量組成一個(gè)nxk維度的矩陣，即為我們的H。一般需要對H矩陣按行做標(biāo)準(zhǔn)化，即

由于我們在使用維度規(guī)約的時(shí)候損失了少量信息，導(dǎo)致得到的優(yōu)化后的指示向量h對應(yīng)的H現(xiàn)在不能完全指示各樣本的歸屬，因此一般在得到nxk維度的矩陣H后還需要對每一行進(jìn)行一次傳統(tǒng)的聚類，比如使用K-Means聚類。

3.3 Ncut切圖

Ncut在最小化損失函數(shù)之外，還考慮了子圖之間的權(quán)重大小。Ncut切圖與Ratiocut類似，只是把RatioCut分母中的

替換成了

，其中：

由于子圖樣本的個(gè)數(shù)多并不一定權(quán)重就大，我們切圖時(shí)基于權(quán)重也更合我們的目標(biāo)，因此一般來說Ncut切圖優(yōu)于RatioCut切圖。所以Ncut的目標(biāo)函數(shù)如下：

，證明同上

在Ncut中使用子圖權(quán)重

來表示指示向量h，定義如下：

我們的優(yōu)化目標(biāo)函數(shù)是：

但是此時(shí)我們的

，而是

。推導(dǎo)如下：

也就是說，此時(shí)我們的優(yōu)化目標(biāo)最終為：

約束條件：

令

，則問題變成如下形式：

約束條件為：

可以發(fā)現(xiàn)這個(gè)式子和RatioCut基本一致，只是中間的L變成了

。這樣我們就可以繼續(xù)按照RatioCut的思想，求出

的最小的前k個(gè)特征值，然后求出對應(yīng)的特征向量，并標(biāo)準(zhǔn)化，得到最后的特征矩陣

,最后對

進(jìn)行一次傳統(tǒng)的聚類（比如K-Means）即可。

一般來說，

相當(dāng)于對拉普拉斯矩陣

做了一次標(biāo)準(zhǔn)化，即

，所以Ncut會(huì)產(chǎn)生標(biāo)準(zhǔn)化的譜聚類，而RatioCut會(huì)產(chǎn)生未標(biāo)準(zhǔn)化的譜聚類。

四. 譜聚類算法的優(yōu)缺點(diǎn)

4.1 優(yōu)點(diǎn)

（1）當(dāng)聚類的類別個(gè)數(shù)較小的時(shí)候，譜聚類的效果會(huì)很好，但是當(dāng)聚類的類別個(gè)數(shù)較大的時(shí)候，則不建議使用譜聚類；

（2）譜聚類算法使用了降維的技術(shù)，所以更加適用于高維數(shù)據(jù)的聚類；

（3）譜聚類只需要數(shù)據(jù)之間的相似度矩陣，因此對于處理稀疏數(shù)據(jù)的聚類很有效。這點(diǎn)傳統(tǒng)聚類算法（比如K-Means）很難做到

（4）譜聚類算法建立在譜圖理論基礎(chǔ)上，與傳統(tǒng)的聚類算法相比，它具有能在任意形狀的樣本空間上聚類且收斂于全局最優(yōu)解

4.2 缺點(diǎn)

（1）譜聚類對相似度圖的改變和聚類參數(shù)的選擇非常的敏感；

（2）譜聚類適用于均衡分類問題，即各簇之間點(diǎn)的個(gè)數(shù)相差不大，對于簇之間點(diǎn)個(gè)數(shù)相差懸殊的聚類問題，譜聚類則不適用；

close all clear all clc; %讀入一副圖像 % ?global slide_window interpixel_distance qutity_level k tic %讀圖 I=imread('oil3-01.jpg'); % I=rgb2gray(I); figure,imshow(I); I1=double(I); % ?I1=rgb2gray(I1); L=harrl(I1); LL=harrll(L); matrix=LL; II=mat2gray(matrix); figure,imshow(II); II=matrix(:); m=size(II,1); %聚類數(shù) k=2; Delta=0.0012; d=zeros(m,m); W=zeros(m,m);%尺度系數(shù) for i=1:m ? ?for j=1:m ? ? ? ?d(i,j)=abs(II(i)-II(j))^2; ? ? ? ?W(i,j)=exp(-d(i,j)/2*Delta^2); ? ?end end ? ? ?% W=affiliate(ent,Delta); Y=cut(W); Y=Y(:,2); % matrix=mat2gray(Y); % MATRIX=reshape(matrix,[60,60]); % figure,imshow(MATRIX); % imwrite(MATRIX,'featureIMAGE'); Y=Y./sum(Y); %Y=Y'; % plot(m,n) %T=0.5*(max(Y)+min(Y)); %while true % ? g=Y>=T; % ?Tn=0.5*(mean(Y(g)+mean(Y(~g)))); % done=abs(T-Tn)<0.1; % ?T=Tn; % end % ?figure,imshow(image); label_matrix1=k_meanclusterings(Y,k); label_matrix2=reshape(label_matrix1,64,64); label_map=label2rgb(label_matrix2); figure,imshow(label_map); imwrite(label_map,'oil3-02.jpg'); toc % im_matric=label_matrix; % ?im_matric1=reshape(im_matric,[m,n]); %label_map=mat2gray(im_matric1); %figure,imshow(label_map);