譜聚類（spectral clustering）是廣泛使用的聚類演算法，比起傳統的K-Means演算法，譜聚類對資料分佈的適應性更強，聚類效果也很優秀，同時聚類的計算量也小很多，更加難能可貴的是實現起來也不復雜。在處理實際的聚類問題時，個人認為譜聚類是應該首先考慮的幾種演算法之一。下面我們就對譜聚類的演算法原理做一個總結。

1. 譜聚類概述

　　　　譜聚類是從圖論中演化出來的演算法，後來在聚類中得到了廣泛的應用。它的主要思想是把所有的資料看做空間中的點，這些點之間可以用邊連線起來。距離較遠的兩個點之間的邊權重值較低，而距離較近的兩個點之間的邊權重值較高，通過對所有資料點組成的圖進行切圖，讓切圖後不同的子圖間邊權重和儘可能的低，而子圖內的邊權重和儘可能的高，從而達到聚類的目的。

　　　　乍一看，這個演算法原理的確簡單，但是要完全理解這個演算法的話，需要對圖論中的無向圖，線性代數和矩陣分析都有一定的瞭解。下面我們就從這些需要的基礎知識開始，一步步學習譜聚類。

2. 譜聚類基礎之一：無向權重圖

　　　　由於譜聚類是基於圖論的，因此我們首先溫習下圖的概念。對於一個圖$G$，我們一般用點的集合$V$和邊的集合$E$來描述。即為$G(V,E)$。其中$V$即為我們資料集裡面所有的點$(v_1, v_2,...v_n)$。對於$V$中的任意兩個點，可以有邊連線，也可以沒有邊連線。我們定義權重$w_{ij}$為點$v_i$和點$v_j$之間的權重。由於我們是無向圖，所以$w_{ij} = w_{ji}$。

　　　　對於有邊連線的兩個點$v_i$和$v_j$，$w_{ij} > 0$,對於沒有邊連線的兩個點$v_i$和$v_j$，$w_{ij} = 0$。對於圖中的任意一個點$v_i$，它的度$d_i$定義為和它相連的所有邊的權重之和，即$$d_i = \sum\limits_{j=1}^{n}w_{ij}$$

　　　　利用每個點度的定義，我們可以得到一個nxn的度矩陣$D$,它是一個對角矩陣，只有主對角線有值，對應第i行的第i個點的度數，定義如下：

$$\mathbf{D} =
\left( \begin{array}{ccc}
d_1 & \ldots & \ldots \\
\ldots & d_2 & \ldots \\  
\vdots & \vdots & \ddots \\  
\ldots & \ldots & d_n \end{array} \right)$$

　　　　利用所有點之間的權重值，我們可以得到圖的鄰接矩陣$W$，它也是一個nxn的矩陣，第i行的第j個值對應我們的權重$w_{ij}$。

　　　　除此之外，對於點集$V$的的一個子集$A \subset V$，我們定義：$$|A|: = 子集A中點的個數$$ $$vol(A): = \sum\limits_{i \in A}d_i$$

3. 譜聚類基礎之二：相似矩陣

　　　　在上一節我們講到了鄰接矩陣$W$，它是由任意兩點之間的權重值$w_{ij}$組成的矩陣。通常我們可以自己輸入權重，但是在譜聚類中，我們只有資料點的定義，並沒有直接給出這個鄰接矩陣，那麼怎麼得到這個鄰接矩陣呢？

　　　　基本思想是，距離較遠的兩個點之間的邊權重值較低，而距離較近的兩個點之間的邊權重值較高，不過這僅僅是定性，我們需要定量的權重值。一般來說，我們可以通過樣本點距離度量的相似矩陣$S$來獲得鄰接矩陣$W$。

　　　　構建鄰接矩陣$W$的方法有三類。$\epsilon$-鄰近法，K鄰近法和全連線法。

　　　　對於$\epsilon$-鄰近法，它設定了一個距離閾值$\epsilon$，然後用歐式距離$s_{ij}$度量任意兩點$x_i$和$x_j$的距離。即相似矩陣的$s_{ij} = ||x_i-x_j||_2^2$,  然後根據$s_{ij}$和$\epsilon$的大小關係，來定義鄰接矩陣$W$如下：

$$W_{ij}=
\begin{cases}
0& {s_{ij} > \epsilon}\\
\epsilon& {{s_{ij} \leq \epsilon}}
\end{cases}$$

　　　　從上式可見，兩點間的權重要不就是$\epsilon$,要不就是0，沒有其他的資訊了。距離遠近度量很不精確，因此在實際應用中，我們很少使用$\epsilon$-鄰近法。

　　　　第二種定義鄰接矩陣$W$的方法是K鄰近法，利用KNN演算法遍歷所有的樣本點，取每個樣本最近的k個點作為近鄰，只有和樣本距離最近的k個點之間的$w_{ij} > 0$。但是這種方法會造成重構之後的鄰接矩陣W非對稱，我們後面的演算法需要對稱鄰接矩陣。為了解決這種問題，一般採取下面兩種方法之一：

　　　　第一種K鄰近法是隻要一個點在另一個點的K近鄰中，則保留$S_{ij}$

$$W_{ij}=W_{ji}=
\begin{cases}
0& {x_i \notin KNN(x_j) \;and \;x_j \notin KNN(x_i)}\\
exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& {x_i \in KNN(x_j)\; or\; x_j \in KNN(x_i})
\end{cases}$$

　　　　第二種K鄰近法是必須兩個點互為K近鄰中，才能保留$S_{ij}$

$$W_{ij}=W_{ji}=
\begin{cases}
0& {x_i \notin KNN(x_j) \;or\;x_j \notin KNN(x_i)}\\
exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& {x_i \in KNN(x_j)\; and \; x_j \in KNN(x_i})
\end{cases}$$

　　　　第三種定義鄰接矩陣$W$的方法是全連線法，相比前兩種方法，第三種方法所有的點之間的權重值都大於0，因此稱之為全連線法。可以選擇不同的核函式來定義邊權重，常用的有多項式核函式，高斯核函式和Sigmoid核函式。最常用的是高斯核函式RBF，此時相似矩陣和鄰接矩陣相同：$$W_{ij}=S_{ij}=exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})$$

　　　　在實際的應用中，使用第三種全連線法來建立鄰接矩陣是最普遍的，而在全連線法中使用高斯徑向核RBF是最普遍的。

4. 譜聚類基礎之三：拉普拉斯矩陣

　　　　單獨把拉普拉斯矩陣(Graph Laplacians)拿出來介紹是因為後面的演算法和這個矩陣的性質息息相關。它的定義很簡單，拉普拉斯矩陣$L= D-W$。$D$即為我們第二節講的度矩陣，它是一個對角矩陣。而$W$即為我們第二節講的鄰接矩陣，它可以由我們第三節的方法構建出。

　　　　拉普拉斯矩陣有一些很好的性質如下：

　　　　1）拉普拉斯矩陣是對稱矩陣，這可以由$D$和$W$都是對稱矩陣而得。

　　　　2）由於拉普拉斯矩陣是對稱矩陣，則它的所有的特徵值都是實數。

　　　　3）對於任意的向量$f$,我們有$$f^TLf = \frac{1}{2}\sum\limits_{i,j=1}^{n}w_{ij}(f_i-f_j)^2$$

　　　　　　這個利用拉普拉斯矩陣的定義很容易得到如下：

$$f^TLf = f^TDf - f^TWf = \sum\limits_{i=1}^{n}d_if_i^2 - \sum\limits_{i,j=1}^{n}w_{ij}f_if_j$$ $$=\frac{1}{2}( \sum\limits_{i=1}^{n}d_if_i^2 - 2 \sum\limits_{i,j=1}^{n}w_{ij}f_if_j + \sum\limits_{j=1}^{n}d_jf_j^2) = \frac{1}{2}\sum\limits_{i,j=1}^{n}w_{ij}(f_i-f_j)^2 $$

　　　　4) 拉普拉斯矩陣是半正定的，且對應的n個實數特徵值都大於等於0，即$0 =\lambda_1 \leq \lambda_2 \leq... \leq \lambda_n$， 且最小的特徵值為0，這個由性質3很容易得出。

5. 譜聚類基礎之四：無向圖切圖

　　　　對於無向圖$G$的切圖，我們的目標是將圖$G(V,E)$切成相互沒有連線的k個子圖，每個子圖點的集合為：$A_1,A_2,..A_k$，它們滿足$A_i \cap A_j = \emptyset$,且$A_1 \cup A_2 \cup ... \cup A_k = V$.

　　　　對於任意兩個子圖點的集合$A, B \subset V$, $A \cap B =  \emptyset$, 我們定義A和B之間的切圖權重為：$$W(A, B) = \sum\limits_{i \in A, j \in B}w_{ij}$$

　　　　那麼對於我們k個子圖點的集合：$A_1,A_2,..A_k$，我們定義切圖cut為：$$cut(A_1,A_2,...A_k) = \frac{1}{2}\sum\limits_{i=1}^{k}W(A_i, \overline{A}_i )$$

 　　　　其中$\overline{A}_i $為$A_i$的補集，意為除$A_i$子集外其他V的子集的並集。

　　　　那麼如何切圖可以讓子圖內的點權重和高，子圖間的點權重和低呢？一個自然的想法就是最小化$cut(A_1,A_2,...A_k)$, 但是可以發現，這種極小化的切圖存在問題，如下圖：

　　　　我們選擇一個權重最小的邊緣的點，比如C和H之間進行cut，這樣可以最小化$cut(A_1,A_2,...A_k)$, 但是卻不是最優的切圖，如何避免這種切圖，並且找到類似圖中"Best Cut"這樣的最優切圖呢？我們下一節就來看看譜聚類使用的切圖方法。

6. 譜聚類之切圖聚類

　　　　為了避免最小切圖導致的切圖效果不佳，我們需要對每個子圖的規模做出限定，一般來說，有兩種切圖方式，第一種是RatioCut，第二種是Ncut。下面我們分別加以介紹。

6.1 RatioCut切圖

　　　　RatioCut切圖為了避免第五節的最小切圖，對每個切圖，不光考慮最小化$cut(A_1,A_2,...A_k)$，它還同時考慮最大化每個子圖點的個數，即：$$RatioCut(A_1,A_2,...A_k) = \frac{1}{2}\sum\limits_{i=1}^{k}\frac{W(A_i, \overline{A}_i )}{|A_i|}$$

　　　　那麼怎麼最小化這個RatioCut函式呢？牛人們發現，RatioCut函式可以通過如下方式表示。

　　　　我們引入指示向量$h_j \in \{h_1, h_2,..h_k\}\; j =1,2,...k$,對於任意一個向量$h_j$, 它是一個n維向量（n為樣本數），我們定義$h_{ij}$為：

$$h_{ij}=
\begin{cases}
0& { v_i \notin A_j}\\
\frac{1}{\sqrt{|A_j|}}& { v_i \in A_j}
\end{cases}$$

　　　　那麼我們對於$h_i^TLh_i$,有：

$$ \begin{align} h_i^TLh_i & = \frac{1}{2}\sum\limits_{m=1}\sum\limits_{n=1}w_{mn}(h_{im}-h_{in})^2 \\& =\frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}(\frac{1}{\sqrt{|A_i|}} - 0)^2 +  \sum\limits_{m \notin A_i, n \in A_i}w_{mn}(0 - \frac{1}{\sqrt{|A_i|}} )^2\\& = \frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}\frac{1}{|A_i|} +  \sum\limits_{m \notin A_i, n \in A_i}w_{mn}\frac{1}{|A_i|}\\& = \frac{1}{2}(cut(A_i, \overline{A}_i) \frac{1}{|A_i|} + cut(\overline{A}_i, A_i) \frac{1}{|A_i|}) \\& =  \frac{cut(A_i, \overline{A}_i)}{|A_i|} \end{align}$$

　　　　上述第（1）式用了上面第四節的拉普拉斯矩陣的性質3. 第二式用到了指示向量的定義。可以看出，對於某一個子圖i，它的RatioCut對應於$h_i^TLh_i$,那麼我們的k個子圖呢？對應的RatioCut函式表示式為：$$RatioCut(A_1,A_2,...A_k) = \sum\limits_{i=1}^{k}h_i^TLh_i = \sum\limits_{i=1}^{k}(H^TLH)_{ii} = tr(H^TLH)$$

　　　　其中$tr(H^TLH)$為矩陣的跡。也就是說，我們的RatioCut切圖，實際上就是最小化我們的$tr(H^TLH)$。注意到$H^TH=I$,則我們的切圖優化目標為:$$\underbrace{arg\;min}_H\; tr(H^TLH) \;\; s.t.\;H^TH=I $$

　　　　注意到我們H矩陣裡面的每一個指示向量都是n維的，向量中每個變數的取值為0或者$\frac{1}{\sqrt{|A_j|}}$，就有$2^n$種取值，有k個子圖的話就有k個指示向量，共有$k2^n$種H，因此找到滿足上面優化目標的H是一個NP難的問題。那麼是不是就沒有辦法了呢？

　　　　注意觀察$tr(H^TLH)$中每一個優化子目標$h_i^TLh_i$,其中$h$是單位正交基， L為對稱矩陣，此時$h_i^TLh_i$的最大值為L的最大特徵值，最小值是L的最小特徵值。如果你對主成分分析PCA很熟悉的話，這裡很好理解。在PCA中，我們的目標是找到協方差矩陣（對應此處的拉普拉斯矩陣L）的最大的特徵值，而在我們的譜聚類中，我們的目標是找到目標的最小的特徵值，得到對應的特徵向量，此時對應二分切圖效果最佳。也就是說，我們這裡要用到維度規約的思想來近似去解決這個NP難的問題。

　　　　對於$h_i^TLh_i$，我們的目標是找到最小的L的特徵值，而對於$tr(H^TLH) = \sum\limits_{i=1}^{k}h_i^TLh_i$，則我們的目標就是找到k個最小的特徵值，一般來說，k遠遠小於n，也就是說，此時我們進行了維度規約，將維度從n降到了k，從而近似可以解決這個NP難的問題。

　　　　通過找到L的最小的k個特徵值，可以得到對應的k個特徵向量，這k個特徵向量組成一個nxk維度的矩陣，即為我們的H。一般需要對H矩陣按行做標準化，即$$h_{ij}^{*}= \frac{h_{ij}}{(\sum\limits_{t=1}^kh_{it}^{2})^{1/2}}$$

　　　　由於我們在使用維度規約的時候損失了少量資訊，導致得到的優化後的指示向量h對應的H現在不能完全指示各樣本的歸屬，因此一般在得到nxk維度的矩陣H後還需要對每一行進行一次傳統的聚類，比如使用K-Means聚類.

6.2 Ncut切圖

　　　　Ncut切圖和RatioCut切圖很類似，但是把Ratiocut的分母$|Ai|$換成$vol(A_i)$. 由於子圖樣本的個數多並不一定權重就大，我們切圖時基於權重也更合我們的目標，因此一般來說Ncut切圖優於RatioCut切圖。$$NCut(A_1,A_2,...A_k) = \frac{1}{2}\sum\limits_{i=1}^{k}\frac{W(A_i, \overline{A}_i )}{vol(A_i)}$$

　　　　，對應的，Ncut切圖對指示向量$h$做了改進。注意到RatioCut切圖的指示向量使用的是$\frac{1}{\sqrt{|A_j|}}$標示樣本歸屬，而Ncut切圖使用了子圖權重$\frac{1}{\sqrt{vol(A_j)}}$來標示指示向量h，定義如下:

$$h_{ij}=
\begin{cases}
0& { v_i \notin A_j}\\
\frac{1}{\sqrt{vol(A_j)}}& { v_i \in A_j}
\end{cases}$$

　　　　那麼我們對於$h_i^TLh_i$,有：

$$ \begin{align} h_i^TLh_i & = \frac{1}{2}\sum\limits_{m=1}\sum\limits_{n=1}w_{mn}(h_{im}-h_{in})^2 \\& =\frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}(\frac{1}{\sqrt{vol(A_j)}} - 0)^2 +  \sum\limits_{m \notin A_i, n \in A_i}w_{mn}(0 - \frac{1}{\sqrt{vol(A_j)}} )^2\\& = \frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}\frac{1}{vol(A_j)} +  \sum\limits_{m \notin A_i, n \in A_i}w_{mn}\frac{1}{vol(A_j)}\\& = \frac{1}{2}(cut(A_i, \overline{A}_i) \frac{1}{vol(A_j)} + cut(\overline{A}_i, A_i) \frac{1}{vol(A_j)}) \\& =  \frac{cut(A_i, \overline{A}_i)}{vol(A_j)} \end{align}$$

　　　　推導方式和RatioCut完全一致。也就是說，我們的優化目標仍然是$$NCut(A_1,A_2,...A_k) = \sum\limits_{i=1}^{k}h_i^TLh_i = \sum\limits_{i=1}^{k}(H^TLH)_{ii} = tr(H^TLH)$$

　　　　但是此時我們的$H^TH \neq I$,而是$H^TDH = I$。推導如下：$$ h_i^TDh_i = \sum\limits_{j=1}^{n}h_{ij}^2d_j =\frac{1}{vol(A_i)}\sum\limits_{v_j \in A_i}w_{v_j} = \frac{1}{vol(A_i)}vol(A_i) =1$$

　　　　也就是說，此時我們的優化目標最終為：$$\underbrace{arg\;min}_H\; tr(H^TLH) \;\; s.t.\;H^TDH=I $$

　　　　此時我們的H中的指示向量$h$並不是標準正交基，所以在RatioCut裡面的降維思想不能直接用。怎麼辦呢？其實只需要將指示向量矩陣H做一個小小的轉化即可。

　　　　我們令$H = D^{-1/2}F$, 則：$H^TLH = F^TD^{-1/2}LD^{-1/2}F$，$H^TDH=F^TF = I$,也就是說優化目標變成了: $$\underbrace{arg\;min}_F\; tr(F^TD^{-1/2}LD^{-1/2}F) \;\; s.t.\;F^TF=I $$

　　　　可以發現這個式子和RatioCut基本一致，只是中間的L變成了$D^{-1/2}LD^{-1/2}$。這樣我們就可以繼續按照RatioCut的思想，求出$D^{-1/2}LD^{-1/2}$的最小的前k個特徵值，然後求出對應的特徵向量，並標準化，得到最後的特徵矩陣$F$,最後對$F$進行一次傳統的聚類（比如K-Means）即可。

　　　　一般來說， $D^{-1/2}LD^{-1/2}$相當於對拉普拉斯矩陣$L$做了一次標準化，即$\frac{L_{ij}}{\sqrt{d_i*d_j}}$

7. 譜聚類演算法流程

　　　　鋪墊了這麼久，終於可以總結下譜聚類的基本流程了。一般來說，譜聚類主要的注意點為相似矩陣的生成方式（參見第二節），切圖的方式（參見第六節）以及最後的聚類方法（參見第六節）。

　　　　最常用的相似矩陣的生成方式是基於高斯核距離的全連線方式，最常用的切圖方式是Ncut。而到最後常用的聚類方法為K-Means。下面以Ncut總結譜聚類演算法流程。

　　　　輸入：樣本集D=$(x_1,x_2,...,x_n)$，相似矩陣的生成方式, 降維後的維度$k_1$, 聚類方法，聚類後的維度$k_2$

　　　　輸出： 簇劃分$C(c_1,c_2,...c_{k_2})$.　

　　　　1) 根據輸入的相似矩陣的生成方式構建樣本的相似矩陣S

　　　　2）根據相似矩陣S構建鄰接矩陣W，構建度矩陣D

　　　　3）計算出拉普拉斯矩陣L

　　　　4）構建標準化後的拉普拉斯矩陣$D^{-1/2}LD^{-1/2}$

　　　　5）計算$D^{-1/2}LD^{-1/2}$最小的$k_1$個特徵值所各自對應的特徵向量$f$

　　　　6) 將各自對應的特徵向量$f$組成的矩陣按行標準化，最終組成$n \times k_1$維的特徵矩陣F

　　　　7）對F中的每一行作為一個$k_1$維的樣本，共n個樣本，用輸入的聚類方法進行聚類，聚類維數為$k_2$。

　　　　8）得到簇劃分$C(c_1,c_2,...c_{k_2})$.　　　　　　　　　

8. 譜聚類演算法總結

　　　　譜聚類演算法是一個使用起來簡單，但是講清楚卻不是那麼容易的演算法，它需要你有一定的數學基礎。如果你掌握了譜聚類，相信你會對矩陣分析，圖論有更深入的理解。同時對降維裡的主成分分析也會加深理解。

　　　　下面總結下譜聚類演算法的優缺點。

　　　　譜聚類演算法的主要優點有：

　　　　1）譜聚類只需要資料之間的相似度矩陣，因此對於處理稀疏資料的聚類很有效。這點傳統聚類演算法比如K-Means很難做到

　　　　2）由於使用了降維，因此在處理高維資料聚類時的複雜度比傳統聚類演算法好。

　　　　譜聚類演算法的主要缺點有：

　　　　1）如果最終聚類的維度非常高，則由於降維的幅度不夠，譜聚類的執行速度和最後的聚類效果均不好。

　　　　2) 聚類效果依賴於相似矩陣，不同的相似矩陣得到的最終聚類效果可能很不同。

（歡迎轉載，轉載請註明出處。歡迎溝通交流： [email protected]）