1. 譜聚類

      給你部落格園上若干個部落格，讓你將它們分成K類，你會怎樣做？想必有很多方法，本文要介紹的是其中的一種——譜聚類。
      聚類的直觀解釋是根據樣本間相似度，將它們分成不同組。譜聚類的思想是將樣本看作頂點，樣本間的相似度看作帶權的邊，從而將聚類問題轉為圖分割問題：找到一種圖分割的方法使得連線不同組的邊的權重儘可能低（這意味著組間相似度要儘可能低），組內的邊的權重儘可能高（這意味著組內相似度要儘可能高）。將上面的例子代入就是將每一個部落格當作圖上的一個頂點，然後根據相似度將這些頂點連起來，最後進行分割。分割後還連在一起的頂點就是同一類了。更具體的例子如下圖所示：

在上圖中，一共有6個頂點（部落格），頂點之間的連線表示兩個頂點的相似度，現在要將這圖分成兩半（兩個類），要怎樣分割（去掉哪邊條）？根據譜聚類的思想，應該去掉的邊是用虛線表示的那條。最後，剩下的兩半就分別對應兩個類了。
      根據這個思想，可以得到unnormalized譜聚類和normalized譜聚類，由於前者比後者簡單，所以本文介紹unnormalized譜聚類的幾個步驟（假設要分K個類）：
(a)建立similarity graph，並用 W 表示similarity graph的帶權鄰接矩陣
(b)計算unnormalized graph Laplacian matrix L(L = D - W, 其中D是degree matrix)
(c)計算L的前K個最小的特徵向量
(d)把這k個特徵向量排列在一起組成一個N*k的矩陣，將其中每一行看作k維空間中的一個向量，並使用 K-means 演算法進行聚類

2. 演算法原理解析

     這一節主要從大體上解釋unnormalized譜聚類的四個步驟是怎麼來的，不涉及具體的公式推導。
(a)譜聚類的思想就是要轉化為圖分割問題。因此，第一步就是將原問題轉化為圖。轉為圖有兩個問題要解決：一是兩個頂點的邊要怎樣定義；二是要保留哪些邊。
      對於第一個問題，如果兩個點在一定程度上相似，就在兩個點之間新增一條邊。相似的程度由邊的權重表示（上圖中邊上面的數值就是權重了）。因此，只要是計算相似度的公式都可用，不過常用的是Gaussian similarity function
      要保留部分邊的原因有：邊太多了不好處理；權重太低的邊是多餘的。常用的保留邊的方法是建立k-nearest neighbor graph。在這種圖中，每個頂點只與K個相似度最高的點連邊。

(b)unnormalized graph Laplacian matrix(以下用L表示)有很多很好的性質，也正是這個原因，才要在第二步中計算這麼一個矩陣。最重要的性質是下面這一組性質：

這一組性質將在之後的公式推導中起到決定性作用。

(c)將原問題轉化為圖後，接下來的工作就是決定怎樣分割了。圖分割問題實際上就是最小割問題(mincut problem)。最小割問題可定義為最小化以下目標函式：

其中k表示分成k個組，Ai表示第i個組，表示第Ai的補集，W(A,B)表示第A組與第B組之間的所有邊的權重之和。
      這個式子的直觀意義：如果要分成K個組，那麼其代價就是進行分割時去掉的邊的權重的總和。可惜的是直接最小化這式子通常會導致不好的分割。以分成2類為例，這個式子通常會將圖分成這樣的兩類：一個點為一類，剩下的所有點為另一類。顯然，這樣的分割是很不好的。因為我們期望著每個類都有合理的大小。所以，要對這個式子進行改進，改進後的公式（稱為RatioCut)如下：

其中|A|表示A組中包含的頂點數目。
      在RatioCut中，如果某一組包含的頂點數越少，那麼它的值就越大。在一個最小化問題中，這相當於是懲罰，也就是不鼓勵將組分得太小。現在只要將最小化RatioCut解出來，分割就完成了。不幸的是，這是個NP難問題。想要在多項式時間內解出來，就要對這個問題作一個轉化了。在轉化的過程中，就用到上面提到的L的那一組性質，經過若干推導，最後可以得到這樣的一個問題：

其中H是一個矩陣，它的元素的定義(Eq.(5))如下：

      如果H矩陣的元素不為0，則說明第i個點屬於第j個類。也就是說，只要得到H矩陣，就能知道要怎樣分割了。可惜的是，這個問題仍然是NP難問題。但是，如果我們讓H矩陣的元素能夠取任意值，這個問題就變成多項式時間內可解的了，此時問題變為：

      根據Rayleigh-Ritz theorem，這個問題的解是L的前k個最小的特徵向量組成的矩陣H，其中特徵向量是按列來排，即H的每一列，均為一個特徵向量。

(d)在第三步中，我們為了鬆馳NP難問題，讓H矩陣取任意值，因此，解出來的H矩陣不再具有原來的性質——元素值能指出哪個點屬於哪一類。儘管如此，對於k-means來說，將H矩陣的每一行當作一個點進行聚類還是挺輕鬆的。因此，用k-means對H矩陣進行聚類作為譜聚類的最終結果。

3. 譜聚類的實現

     以下是unnormalized譜聚類的MATLAB版實現（部落格園的程式碼格式選擇中居然沒有Matlab的。。。這裡選個C++的）：

function [ C, L, D, Q, V ] = SpectralClustering(W, k)
% spectral clustering algorithm
% input: adjacency matrix W; number of cluster k 
% return: cluster indicator vectors as columns in C; unnormalized Laplacian L; degree matrix D;
%         eigenvectors matrix Q; eigenvalues matrix V

% calculate degree matrix
degs = sum(W, 2);
D = sparse(1:size(W, 1), 1:size(W, 2), degs);

% compute unnormalized Laplacian
L = D - W;

% compute the eigenvectors corresponding to the k smallest eigenvalues
% diagonal matrix V is NcutL's k smallest magnitude eigenvalues 
% matrix Q whose columns are the corresponding eigenvectors.
[Q, V] = eigs(L, k, 'SA');

% use the k-means algorithm to cluster V row-wise
% C will be a n-by-1 matrix containing the cluster number for each data point
C = kmeans(Q, k);

% convert C to a n-by-k matrix containing the k indicator vectors as columns
C = sparse(1:size(D, 1), C, 1);

end

4. 相關資料