作者：桂。
時間：2017-04-13 19:14:48
鏈接：http://www.cnblogs.com/xingshansi/p/6702174.html
聲明：本文大部分內容來自：劉建平Pinard博客的內容。
前言
之前對非負矩陣分解（NMF）簡單梳理了一下，總覺得NMF與聚類非常相似，像是譜聚類的思想。在此將譜聚類的知識梳理一下，內容無法轉載，不然直接轉載劉建平Pinard的博文了，常用的譜聚類有RatioCut和Ncut算法，全文主要梳理RatioCut算法：
1）背景知識；
2）理論推導；
3）應用實例
內容為自己的學習記錄，其中參考他人的部分，最后一並給出鏈接。
一、背景知識
關於圖的基本概念，以及常用到的拉普拉斯矩陣，之前已經有博文介紹過。直接從圖的分割說起：
A-鄰接矩陣
鄰接矩陣的構造方法，常用的有KNN、全連接等方法，這里僅以全連接中的高斯核為例：
$W_{ij}=S_{ij}=exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})$
B-無向圖切圖
對於無向圖$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$的補集。
那么如何切圖可以讓子圖內的點權重和高，子圖間的點權重和低呢？一個自然的想法就是最小化$cut(A_1,A_2,...A_k)$，但是可以發現，這種極小化的切圖存在問題，如下圖：
找到類似圖中"Best Cut"這樣的最優切圖呢？一個自然的想法就是，類似為了防止過擬合而添加正則項一樣，可以添加新的限定，這就是譜聚類的思想。
二、理論推導（RatioCut）
定義$|A_i|$: = 子集$A_i$中點的個數。現在對每個切圖，不光考慮最小化$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 =\{h_1, h_2,..h_k\}\; j =1,2,...k$，對於任意一個向量$h_j$它是一個n維向量（n為樣本數），我們定義$h_{ji}$為：
$h_{ji}= \begin{cases} 0& { v_i \notin A_j}\\ \frac{1}{\sqrt{|A_j|}}& { v_i \in A_j} \end{cases}$
借助拉普拉斯矩陣特性，我們對於$h_i^TLh_i$有：
可以看出，對於某一個子圖i，它的RatioCut對應於$h_i^TLh_i$，那么我們的k個子圖呢？對應的RatioCut函數表達式為：
注意到$H^TH=I$，優化函數轉化為：
因為每一個h的取值有兩種可能，因此該准則函數需要k*2n種H，這是一個NP難問題。
如果對條件適當放松呢？比如這樣：
h不再看作只有兩種取值的離散變量，而是具有連續取值的變量。
這樣一來，上面的優化函數就可以對h利用拉格朗日乘子法進行求解。這種求解方法是瑞利熵求解的一類，關於瑞利熵前文有介紹。因為這里放寬了h的限定，使得h從離散量變為連續量，如何與之前的對應呢？最簡單的辦法就是看求解的h離h原始的兩個取值，哪個更近，對應的就算做哪一類。離哪個更近？沒錯，這正是Kmeans的思想，故后處理也可以用調Kmeans來完成。Kmeans之前，通常將求解的h每一列分別歸一化。
至此完成了RatioCut的步驟。
三、代碼實現
首先根據上文的理論分析，給出RatioCut的算法步驟：
步驟一：求解拉普拉斯矩陣L
步驟二：對L進行特征值分解，並取K個最小特征值對應的特征向量（K為類別數目）
步驟三：將求解的K個特征向量（並分別歸一化），構成新的矩陣，對該矩陣進行Kmeans處理
Kmeans得到的類別標簽，就是原數據的類別標簽，至此完成RatioCut聚類。
給出對應代碼：
sigma2 = 0.002;
%%Step1: Calculate Laplace matrix
for i = 1:N
for j =1:N
W(i,j) = exp(-sqrt(sum((X(i,:)-X(j,:)).^2))/2/sigma2);
end
end
W = W-diag(diag(W));% adjacency matrix
D = diag(sum(W)); %degree matrix
L = D-W;%laplace matrix
%%Step2:Eigenvalues decomposition
K = 3;
[Qini,V] = eig(L);
%%Step3:New matrix Q
[~,pos] = sort(diag(V),'ascend');
Q = Qini(:,pos(1:K));
Q = Q./repmat(sqrt(diag(Q'*Q)'),N,1);
[idx,ctrs] = kmeans(Q,K);
測試一下，按數據為3類進行譜聚類，可以看出來還是有效的，譜聚類中高斯權重涉及到$\sigma$如何取值，不過這里就不做進一步討論了。
參考：