引入
聚類演算法一般可以分為兩類：
- Compactness。代表的演算法有 K-means，GMM 等。但這類演算法只能處理凸集，為了處理非凸的樣本集，必須引⼊核技巧。
- Connectivity。這類以 spectral clustering 為代表。
舉個例子，將下述資料採用聚類演算法進行聚類，可以採用 GMM 或 K-Means 的方法：
然而對於下述資料卻並不能使用上述兩種演算法：
此時可以考慮採用譜聚類（spectral clustering）的方法。
譜聚類演算法(Spectral Clustering)
主要思想是把所有的資料看做空間中的點，這些點之間可以用邊連線起來。距離較遠的兩個點之間的邊權重值較低，而距離較近的兩個點之間的邊權重值較高，通過對所有資料點組成的圖進行切圖，讓切圖後不同的子圖間邊權重和儘可能的低，而子圖內的邊權重和儘可能的高，從而達到聚類的目的。
1 基礎知識
1.1 無向帶權圖
對於一個圖 $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}$
利用節點的權重值，可以得到圖的鄰接矩陣 $W$ ，是一個 $n \times n$ 的矩陣，第 $i$ 行的第 $j$ 個值對應的權重為 $w_{ij}$。
除此之外，對於點集 $V$ 的的一個子集 $A \subset V$ ，我們定義：
-
- $|A|: = 子集A中點的個數$
- $ vol(A): = \sum\limits_{i \in A}d_i$
利用每個節點的度，可以得到一個 $n \times n$ 的度矩陣 $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)$
1.2 相似矩陣
一般情況下，鄰接矩陣 $W$，通常可以自己輸入權重，但在譜聚類中，我們只有資料點的定義，並沒有直接給出這個鄰接矩陣，那麼怎麼得到這個鄰接矩陣呢？
基本思想：距離較遠的兩個點之間的邊權重值較低，而距離較近的兩個點之間的邊權重值較高。一般來說，可以通過樣本點距離度量的相似矩陣 $S$ 來獲得鄰接矩陣 $W$。
構建鄰接矩陣 $W$ 的方法有三類：
- $\epsilon$ - 鄰近法
- $K$ 鄰近法
- 全連線法
1.2.1 $\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$- 鄰近法。
1.2.2 $K$ 鄰近法
基本思想：利用 KNN 演算法遍歷所有的樣本點，取每個樣本最近的 $k$ 個點作為近鄰，只有和樣本距離最近的 $k$ 個點之間的 $w_{ij} > 0$ 。
但是這種方法會造成重構之後的鄰接矩陣 $W$ 非對稱，且後面的演算法需要對稱鄰接矩陣。為解決這種問題，一般採取下面兩種方法之一：
Method1：第一種 $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}$
Method2：第二種 $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}$
1.2.3 全連線法
相比前兩種方法，第三種方法所有的點之間的權重值都大於 $0$，因此稱之為全連線法。
基本思想：選擇不同的核函式來定義邊權重，常用的有多項式核函式，高斯核函式和 Sigmoid 核函式。
最常用的是高斯核函式RBF，此時相似矩陣和鄰接矩陣相同：
$w_{ij}=s_{ij}=exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})$
在實際的應用中，使用第三種全連線法來建立鄰接矩陣是最普遍的，而在全連線法中使用高斯徑向核 RBF 是最普遍的。
2 拉普拉斯矩陣
拉普拉斯矩陣 $L= D-W$。$D$ 即為度矩陣，它是一個對角矩陣。而 $W$ 即為 鄰接矩陣。（參考《圖神經網路基礎二：譜圖理論》）
舉例：
普拉斯矩陣的性質如下：
1、拉普拉斯矩陣是 對稱矩陣;
2、$L$ 的行和為零 ；
3、$\mathbf{L} $ 有一個特徵值為零 ；
4、$\mathrm{L} $ 是半正定矩陣 ；
5、對於任意向量 $f$，有：
$f^{T} L f=\frac{1}{2} \sum \limits _{i=1}^{N} \sum \limits _{j=1}^{N} w_{i j}\left(f_{i}-f_{j}\right)^{2}$
這一性質利用拉普拉斯矩陣的性質很容易可以得到：
$\begin{array}{l}f^{T} L f&=f^{T} D f-f^{T} W f \\&=\sum \limits _{i=1}^{N} d_{i} f_{i}^{2}-\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i} f_{j} \\&=\frac{1}{2}\left(\sum\limits_{i=1}^{N} d_{i} f_{i}^{2}-2 \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i} f_{j}+\sum\limits_{j=1}^{N} d_{j} f_{j}^{2}\right) \\&=\frac{1}{2}\left(\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i}^{2}-2 \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i} f_{j}+\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{j}^{2}\right) \\&=\frac{1}{2} \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j}\left(f_{i}-f_{j}\right)^{2}\end{array}$
對於上面的性質5，如果 $f$ 為網路中訊號的值的向量，那麼 $f^{T} L f$ 稱為圖訊號的總變差 (Total Variation)，可以刻畫圖訊號整體的平滑度
3 無向圖切圖
對於無向圖 $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$ 。然而這種方法存在問題，也就是最小化的 $cut$ 對應的切圖不一定就是符合要求的最優的切圖，如下圖：
在上面的例子中，我們選擇一個權重最小的邊緣的點，比如 $C$ 和 $H$ 之間進行 $cut$，這樣可以最小化 $cut(A_1,A_2,...A_k)$，但是卻不是最優的切圖，如何避免這種切圖，並且找到類似圖中 "Best Cut" 這樣的最優切圖呢？
接下介紹譜聚類使用的切圖方法。
4 譜聚類之切圖聚類
為避免最小切圖導致的切圖效果不佳，需要對每個子圖的規模做出限定，一般來說，有兩種切圖方式，
- 第一種是 RatioCut
- 第二種是 Ncut
4.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 函式，引入指示向量 $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{array}{l} 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{array}$
可以看出，對於某一個子圖 $i$，它的 RatioCut 對應於 $h_i^TLh_i$，那麼 $k$個子圖呢？對應的 RatioCut 函式表示式為：
$\begin{array}{l}RatioCut \left(A_{1}, A_{2}, \ldots A_{k}\right)&=\sum\limits _{i=1}^{k} h_{i}^{T} L h_{i}\\&=\sum\limits_{i=1}^{k}\left(H^{T} L H\right)_{i i} \\&=\operatorname{tr}\left(H^{T} L H\right)\end{array}$
上式中 $\operatorname{tr}\left(H^{T} L H\right) $ 為矩陣 $H^{T} L H$ 的跡， $H=\left(\begin{array}{llll}h_{1} & h_{2} & \cdots & h_{k}\end{array}\right) $ ，需要注意這裡 的 $H$ 滿足 $H^{T} H=I $ ，並且 $H$ 的元素只能取 $0$ 或者 $\frac{1}{\left|A_{i}\right|} $ 。
所以我們需要優化以下目標函式：
$\begin{array}{c} \underset{H}{\operatorname{argmin}} \quad \operatorname{tr}\left(H^{T} L H\right) \\\text { s.t. } H^{T} H=I\end{array}$
注意到 $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 ，我們的目標是找到協方差矩陣（對應此處的拉普拉斯矩陣 $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$ 現在不能完全指示各樣本的歸屬，因此一般在得到 $n\times k$ 維度的矩陣 $H$ 後還需要對每一行進行一次傳統的聚類，比如使用 K-Means 聚類。
4.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{array}{l} 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_i)}} - 0)^2 + \sum\limits_{m \notin A_i, n \in A_i}w_{mn}(0 - \frac{1}{\sqrt{vol(A_i)}} )^2\\& = \frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}\frac{1}{vol(A_i)} + \sum\limits_{m \notin A_i, n \in A_i}w_{mn}\frac{1}{vol(A_i)}\\& = \frac{1}{2}(cut(A_i, \overline{A}_i) \frac{1}{vol(A_i)} + cut(\overline{A}_i, A_i) \frac{1}{vol(A_i)}) \\& = \frac{cut(A_i, \overline{A}_i)}{vol(A_i)} \end{array}$
推導方式和 RatioCut 完全一致。也就是說，我們的優化目標仍然是
$\begin{array}{l}\operatorname{NCut}\left(A_{1}, A_{2}, \cdots, A_{k}\right)&=\sum\limits _{i=1}^{k} h_{i}^{T} L h_{i} \\&=\sum\limits_{i=1}^{k}\left(H^{T} L H\right)_{i i} \\&=\operatorname{tr}\left(H^{T} L H\right)\end{array}$
但是此時我們的 $H^TH \neq I$，而是 $H^TDH = I$。推導如下：
$\begin{array}{l}H^{T} D H&=\left(\begin{array}{c}h_{1}^{T} \\h_{2}^{T} \\\vdots \\h_{k}^{T}\end{array}\right)\left(\begin{array}{ccccc}d_{1} & & & \\& d_{2} & & \\& & \ddots & \\& & & d_{N}\end{array}\right)\left(\begin{array}{cccc}h_{1} & h_{2} & \cdots & h_{k}\end{array}\right) \\&=\left(\begin{array}{cccc}h_{11} d_{1} & h_{12} d_{2} & \cdots & h_{1 N} d_{N} \\h_{21} d_{1} & h_{22} d_{2} & \cdots & h_{2 N} d_{N} \\\vdots & \vdots & \ddots & \vdots \\h_{k 1} d_{1} & h_{k 2} d_{2} & \cdots & h_{k N} d_{N}\end{array}\right)\left(\begin{array}{cccc}h_{1} & h_{2} & \cdots & h_{k}\end{array}\right) \\&=\left(\begin{array}{ccccc}\sum\limits _{i=1}^{N} h_{1 i}^{2} d_{i} & \sum\limits_{i=1}^{N} h_{1 i} h_{2 i} d_{i} & \cdots & \sum\limits_{i=1}^{N} h_{1 i} h_{k i j} d_{i} \\\sum\limits_{i=1}^{N} h_{2 i} h_{1 i} d_{i} & \sum\limits_{i=1}^{N} h_{2 i}^{2} d_{i} & \cdots & \sum\limits_{i=1}^{N} h_{2 i} h_{k i} d_{i} \\\vdots & \vdots&\ddots & \vdots \\\sum\limits_{i=1}^{N} h_{k i} h_{1 i} d_{i} & \sum\limits_{i=1}^{N}h_{k i} h_{2 i} d_{i} & \cdots & \sum\limits_{i=1}^{N} h_{k i}^{2} d_{i}\end{array}\right)\end{array}$
對於對角線元素
$ h_i^TDh_i = \sum\limits_{j=1}^{n}h_{ij}^2d_j =\frac{1}{vol(A_i)}\sum\limits_{j \in A_i}d_j= \frac{1}{vol(A_i)}vol(A_i) =1$
由於 $h_{m i}$ 和 $h_{n i}$ 不可能同時非零（$v_i \in V_m \quad and \quad v_i \in V_n$ ），因此對於非對角線元素有:
$\sum \limits _{i=1}^{N} h_{m i} h_{n i} d_{i}=\sum \limits _{i=1}^{N} 0 \cdot d_{i}=0$
此時我們的優化目標最終為：
$\begin {array} {c} \underset{H}{arg\;min} \; tr(H^TLH) \\s.t.\;H^TDH=I \end {array} $
此時我們的 $H$ 中的指示向量 $h$ 並不是標準正交基，所以在 RatioCut裡面的降維思想不能直接用。怎麼辦呢？其實只需要將指示向量矩陣 $H$ 做一個小小的轉化即可。
我們令$H = D^{-1/2}F$, 則：
$H^TLH = F^TD^{-1/2}LD^{-1/2}F$
$H^TDH=F^TF = I$
此時優化目標為:
$\begin {array}{c}\underset{F}{\operatorname{arg\;\;min}}\;\; \operatorname{tr}\left(F^{T} D^{-1 / 2} L D^{-1 / 2} F\right) \\\text { s.t. } F^{T} F=I\end {array}$
可以發現這個式子和 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$ 做了一次標準化，即
$\left(D^{-1 / 2} L D^{-1 / 2}\right)_{i j}=\frac{L_{i j}}{\sqrt{d_{i} * d_{j}}} $
5 譜聚類演算法流程
5.1 Ncut 譜聚類演算法流程
下面 Ncut 譜聚類演算法流程。
輸入：樣本集D=$(x_1,x_2,...,x_n)$，相似矩陣的生成方式, 降維後的維度$k_1$, 聚類方法，聚類後的維度$k_2$
輸出： 簇劃分$C(c_1,c_2,...c_{k_2})$.
- 根據輸入的相似矩陣的生成方式構建樣本的相似矩陣 $S$ ；
- 根據相似矩陣 $S$ 構建鄰接矩陣 $W$，構建度矩陣 $D$；
- 計算出拉普拉斯矩陣 $L$；
- 構建標準化後的拉普拉斯矩陣$D^{-1/2}LD^{-1/2}$；
- 計算$D^{-1/2}LD^{-1/2}$最小的 $k_1$ 個特徵值所各自對應的特徵向量 $f$；
- 將各自對應的特徵向量$f$組成的矩陣按行標準化，最終組成 $n \times k_1$ 維的特徵矩陣 $F$；
- 對F中的每一行作為一個$k_1$維的樣本，共 $n$ 個樣本，用輸入的聚類方法進行聚類，聚類維數為$k_2$ ；
- 得到簇劃分$C(c_1,c_2,...c_{k_2})$.
5.2 譜聚類演算法的優缺點
譜聚類演算法的主要優點有：
- 譜聚類只需要資料之間的相似度矩陣，因此對於處理稀疏資料的聚類很有效。這點傳統聚類演算法比如K-Means很難做到
- 由於使用了降維，因此在處理高維資料聚類時的複雜度比傳統聚類演算法好。
譜聚類演算法的主要缺點有：
- 如果最終聚類的維度非常高，則由於降維的幅度不夠，譜聚類的執行速度和最後的聚類效果均不好。
- 聚類效果依賴於相似矩陣，不同的相似矩陣得到的最終聚類效果可能很不同。