Spectral Clustering 从入门到入门

0x00 Introduction

写这篇文章的由头是因为自己前一阵读了Thomas N. Kipf, Max Welling大佬的Semi-Supervised Classification with Graph Convolutional Networks 这篇文章。 文章核心在于怎么用数学的方式减小用在Graph数据上神经网络的运算量。因为之前所接触的神经网络限于Andrew NG的两门课，故而没接触过处理Graph的情况。这里的Graph指点集，边集构成的图。这篇文章并非paper的讲解。只是记录一下我关于Graph Clustering这个领域的入门心得。

Figure 1

0x01 Graph Clustering

直觉上: Graph Clustering是把一个大的图分割成多个子图。如图Figure 1 所示，我们希望cluster之间的相似度尽可能小，希望同一个cluster里点的相似度尽可能大。
数学上: 我们可以写出这样的函数:

$Cut(A_1,A_2,...,A_k) :=  \frac{1}{2}  \sum_{i=1}^{k}W(A_i,\bar A_i)$

这里的$\bar A_i$表示集合A的补集，  $W(A,B)$ 代表对所有属于A集合的点i 到所有属于B集合的点j 相似度求和。

$W(A,B) := \sum_{i\in A,j\in B}w_{ij}$

于是目标就变成了最小化Cut函数。

但是可以想象如果把每个节点都设为一个cluster的话Cut函数一定会取得最小值，但这并不是我们想要的。
改进：
直觉上：我们希望每个cluster尽可能大。
数学上：我们可以加上一个约束项，对于'大'的不同理解我们有两种改进方式。

$RatioCut(A_1,A_2,...,A_k) := \frac{1}{2}\sum_{i=1}^{k}\frac{W(A_i,\bar A_i)}{\left | A_i \right |} = \sum_{i=1}^{k}\frac{cut(A_i,\bar A_i)}{\left | A_i \right |}$

$Ncut(A_1,A_2,...,A_k) := \frac{1}{2}\sum_{i=1}^{k}\frac{W(A_i,\bar A_i)}{vol(A_i)} = \sum_{i=1}^{k}\frac{cut(A_i,\bar A_i)}{vol(A_i)}$

RatioCut 认为集合里面包含的点的数量多叫做'大'；

Ncut认为集合里面边权重和越大，那么集合越'大'；

0x02 Laplacian Matrix

需要先把刚才的内容放在一边，生硬的先讲拉普拉斯矩阵。根据定义 $L= D - A$ ，L就是拉普拉斯矩阵。

Figure 2

定义里面D 代表一个图的Degree Matrix，A代表图的Adajacency Matrix.

Degree Matrix, D:

$\begin{bmatrix}  2 &0  & 0 & 0 & 0 &0 \\  0& 3 & 0 & 0 & 0 & 0\\  0&  0& 2 & 0 & 0 &0 \\  0& 0 & 0 &3  & 0 &0 \\  0&  0&  0& 0 &3  &0 \\  0&  0& 0 &0  & 0 &1 \end{bmatrix}$

Adajacency Matrix, A:

$\begin{bmatrix}  0 & 1 & 0 & 0 & 1 &0 \\  1& 0 & 1 & 0 & 1 & 0\\  0&  1& 0 & 1 & 0 &0 \\  0& 0 & 1 &0  & 1 &1 \\  1&  1&  0& 1 &0  &0 \\  0&  0& 0 & 1 & 0 &0 \end{bmatrix}$

Then D - A = L , L:
$\begin{bmatrix}  2 & -1 & 0 & 0 & -1 &0 \\  -1& 3 & -1 & 0 & -1 & 0\\  0&  -1& 2 & -1 & 0 &0 \\  0& 0 & -1 &3  & -1 &-1 \\  -1&  -1&  0& -1 &3  &0 \\  0&  0& 0 & -1 & 0 &1 \end{bmatrix}$

拉普拉斯矩阵有如下性质：

• L是对称半正定矩阵 (symmetric and positive semi-definite matrix)，对于每个 $f \in \mathbb{R}^n$ ，下面式子成立:

$f'L f=\frac{1}{2}\sum_{i,j=1}^{N}w_{ij}(f_i-f_j)^2$

$Proof$ :

             $f'Lf$

        $= f'Df - f'Wf$

        $= \sum_{i=1}^{n}d_if_i^2 - \sum_{i,j=1}^{n}f_if_jw_{ij}$

        $=\frac{1}{2}(\sum_{i=1}^{n}d_if_i^2 - 2\sum_{i,j=1}^{n}f_if_jw_{ij} + \sum_{j=1}^{n}d_jf_j^2 )$

        $=\frac{1}{2}\sum_{i,j=1}^{n}w_{ij}(f_i-f_j)^2$

• L有n个非负，实数特征值，最小的为0, 最小的特征值$0$所对应的特征向量是$\vec{1}$

$0=\lambda_1 \leq \lambda_2 ... \leq \lambda_n$

0x03 Eigenvalue & Eigenvectors 特征值 & 特征向量

这里A是一个矩阵，$\vec{v}$是一个向量，$\lambda$是一个实数

$A\vec{v}  = \lambda\vec{v}$

满足上式的$\vec{v}$是特征向量，$\lambda$为其对应的特征值

0x04 Connect RatioCut to Laplacian Matrix

目标是 min $RatioCut(A_1,A_2...A_j...A_k)$

为此我们引入{$A_1,A_2...A_j...A_k$}的指示向量 $h_j = \{h_1,h_2...h_i...h_n\}$ ,  j = 1...k, i = 1... N; j表示子集下标, i 表示样本下标, 具体的值为:

$h_{j,i}=\begin{cases} \frac{1}{\sqrt{|A_j|}} & \text{ if } v_i \in A_j\\ 0 & \text{ if } v_i  \notin A_j \end{cases}$

每一个子集$A_j$对应一个指示向量$h_j$，每个$h_j$里面有N个元素分别代表N个样本点的指示结果。如果在原始数据中第i个样本被分割到了子集$A_j$中，那么$h_j$的第i个元素值为$\frac{1}{\sqrt{|A_j|}}$,否则为0.
这里$h_j$是一个长度为N的向量

$\begin{align*} h_j^TLh_j &=h_j^T(D-W)h_j \\  &=h_j^TDh_j - h_j^TWh_j \\  &= \sum_{m=1}\sum_{n=1}h_{j,m}h_{j,n}D_{m,n} -\sum_{m=1}\sum_{n=1}h_{j,m}h_{j,n}w_{m,n} \\  &=\sum_{m=1}h_{j,m}^2D_{m,m} - \sum_{m=1}\sum_{n=1}h_{j,m}h_{j,n}w_{m,n} \\  &= \frac{1}{2}(\sum_{m=1}h_{j,m}^2D_{m,m} - 2 \sum_{m=1}\sum_{n=1}h_{j,m}h_{j,n}w_{m,n} + \sum_{n=1}h_{j,n}^2D_{n,n}) \end{align*}$

步骤解释：
第一第二步：拉普拉斯变换
第三步：把矩阵的乘法变成element-wise乘法
第四步：由于D是对角阵，非主对角线元素为0。
第五步：代数变换

然后：

$\begin{align*} h_j^TLh_j &= \frac{1}{2}(\sum_{m=1}h_{j,m}^2D_{m,m} - 2 \sum_{m=1}\sum_{n=1}h_{j,m}h_{j,n}w_{m,n} + \sum_{n=1}h_{j,n}^2D_{n,n}) \\ &=\frac{1}{2}(\sum_{m=1}\sum_{n=1}h_{j,m}^2w_{m,n} - 2\sum_{m=1}\sum_{n=1}h_{j,m}h_{j,n}w_{m,n} + \sum_{n=1}\sum_{m=1}h_{j,n}^2w_{n,m})\\ &=\frac{1}{2}\sum_{m=1}\sum_{n=1}(h_{j,m}-h_{j,n})^2 \end{align*}$

步骤解释：
第一步到第二步：$D_{m,m} = \sum_{n=1}w_{m,n}$，即把Degree 矩阵还原成邻接矩阵来表达
第二步到第三步：对于$w$矩阵而言，里面的元素非0即1；二次展开
进一步, 把$h_{j,i}$ 代入上式:
$\begin{align*} h_j^T L h_j &=\frac{1}{2}\sum_{m=1}\sum_{n=1}(h_{j,m}-h_{j,n})^2 \\  &=\frac{1}{2}(\sum_{m \in A_j,n \in \bar A_j}w_{m,n}(\frac{1}{\sqrt{|A_j|}} - 0)^2 +  \sum_{m \in \bar A_j,n \in  A_j}w_{m,n}( 0 - \frac{1}{\sqrt{|A_j|}})^2)  \\  &=\frac{1}{2}(\sum_{m \in A_j,n \in \bar A_j}w_{m,n} \frac{1}{|A_j|} + \sum_{m \in \bar A_j,n \in  A_j}w_{m,n} \frac{1}{|A_j|}) \\  &= \frac{1}{2} (cut(A_j,\bar A_j) \frac{1}{|A_j|}+ cut(\bar A_j,A_j)\frac{1}{|A_j|})\\  &= \frac{cut(A_j,\bar A_j)}{|A_j|} \\  &= RatioCut(A_j,\bar A_j) \end{align*}$
推到这里，很显然$h_{j,i}$的奇怪的值是为了得到这个结果而凑出来的。(摊手)
进一步，因为在引入指示变量的时候，每一个元素都只能放在一个类里面，所以可以得到这个性质：$h_i \ast h_j = 0, i \neq j$ , 且 $h_i \ast h_i = 1$，即指示向量之间正交，于是：
$$\begin{align*} RatioCut(A_1,A_2...A_k) &= \sum_{i=1}^kh_i^T L h_i\\  &=\sum_{i=1}^k(H^TLH)_{ii} \\  &= Tr(H^TLH) \end{align*}$$
这里随后一步解释一下，Tr ()表示求矩阵的Trace的和，$H^TLH$ 是一个k*k 的矩阵。一个矩阵的Trace等于这个矩阵特征值的总和。[5]
另外一个需要解释的点在于：
$Lh_j = \lambda h_j$ 这里$h_j, \lambda$可以理解为L的一组特征值特征向量，在等式两个都左乘一个$h_j^{'}$ 即 $h_j$的转置：
$h_j^{'} L h_j = \lambda h_j^{'} h_j \\     (Here : h_j^{'} h_j = 1 ) \\     h_j^{'} L h_j = \lambda$
所以我们可以用特征值求和，也就是Trace，来代表$RatioCut(A_1,A_2...A_k)$

所以现在的问题从 $arg\min_{H}Ratiocut(A_1,A_2...A_k)$ 变为 $arg\min_{H}Tr(H^TLH) , s.t. H^TH = I$
在之前的推导过程中我们人为设定了一组特征向量, 即指示向量$h_j$，其中的元素要么是0，要么是$\frac{1}{\sqrt{|A_j|}}$，如果我们用同样的方式去求H，那为了得到最优解，需要遍历$(C_k^1)^N$种情况。

于是我们想不通过特征向量，直接求特征值。

既然我们要分成k类，那我们就对L求前k个特征值，这里说的前后指的是把特征值按大小排序，取前k小的。然后求和以后就是最有的RatioCut解，然后我们根据这个特征值再去得到特征向量。根据特征向量来决定聚类。
但是这时候对应的特征向量, 比如$h_j$里面的元素就不是离散的值了，而是连续的值。

解决这个问题的一个方式是，我们得到了$h_1,h_2...h_k$ 之后用对于$h_j$ 里面的元素，看它距离0或者$\frac{1}{\sqrt{|A_j|}}$哪个近，就改成哪个值。在变换完之后，我们相当于得到了k个中心，然后把N个点和k个中心 进行一次K-means。

类似的对于N-Cut的推导在[3]里面可以看到。

ref:

[1] https://blog.csdn.net/v_JULY_v/article/details/40738211
[2] https://arxiv.org/pdf/0711.0189.pdf
[3] https://blog.csdn.net/yc_1993/article/details/52997074
[4] https://www.cnblogs.com/xingshansi/p/6702174.html
[5] https://zh.wikipedia.org/wiki/%E8%B7%A1