一、前言
本文将介绍一种较为突出的聚类算法——谱聚类。这是一种基于图模型的算法。有过图论或者数据结构基础的人估计会相对来说觉得容易。能力有限,推导并不严谨,请见谅。代码谱聚类代码实现
二、相关知识
2.1、无向图
无向图,由两个要素构成——节点和边。以上面的无向图为例
节点,我们定义为 V = { A , B , C , D , E } V=\{A,B,C,D,E\} V={A,B,C,D,E}。边定义为 E = { w 1 , w 2 , w 3 , w 4 , w 5 } E=\{w_1,w_2,w_3,w_4,w_5\} E={w1,w2,w3,w4,w5}对应图中的w1…。
所谓边,就是两个点之间的关联性,一般情况下,关联性越高,所连接的边 w w w就越高。
在谱聚类中,我们的一个样本点就是一个节点,那么如果我们有n个样本,那就有n个节点。样本之间的相关性就是用 w w w来衡量。既然是聚类,那自然就是相关性越高,则越能聚成一类。
2.2、邻接矩阵
前面我们讲到,某个节点与某个节点之间连接的值用
w
w
w来表示。那有没有一种东西,能够充分表达所有点之间的关系?有,就是邻接矩阵。比如图中就可以表示为(无连接我们暂时用0表示)
W
=
[
A
B
C
D
E
A
0
w
1
w
2
0
0
B
w
1
0
0
w
3
0
C
w
2
0
0
0
w
4
D
0
w
3
0
0
w
5
E
0
0
0
w
4
w
5
]
5
×
5
W=\left[\begin{array}{c|ccccc} &A &B &C &D &E\\\hline A & 0 & w_1 & w_2 & 0 & 0\\ B & w_1 & 0 & 0 & w_3 & 0\\ C & w_2 & 0 & 0 & 0 & w_4\\ D & 0 & w_3 & 0 & 0 & w_5\\ E & 0 & 0 & 0 & w_4 & w_5 \end{array}\right]_{5\times5}
W=
ABCDEA0w1w200Bw100w30Cw20000D0w300w4E00w4w5w5
5×5
比如节点A,它与B,C有连接,所以在B,C列有值,其余都为0。以此类推。
如果有n个样本,那么它就是一个
n
×
n
n\times n
n×n的方阵。所以我们可以表示成这样
W
=
[
w
11
w
12
⋯
w
1
n
w
21
w
22
⋯
w
2
n
⋮
⋮
⋱
⋮
w
n
1
w
n
2
⋯
w
n
n
]
n
×
n
W=\begin{bmatrix} w_{11} & w_{12} & \cdots & w_{1n}\\ w_{21} & w_{22} & \cdots & w_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ w_{n1} & w_{n2} & \cdots & w_{nn} \end{bmatrix}_{n\times n}
W=
w11w21⋮wn1w12w22⋮wn2⋯⋯⋱⋯w1nw2n⋮wnn
n×n
2.3、度
度,意为某个节点所有相连的边的总和。比如图中的节点A,它的度
d
=
w
1
+
w
2
d=w_1+w_2
d=w1+w2。如果有n个样本,那么第
i
i
i个样本的度可以表示为
d
i
=
∑
j
=
1
n
w
i
j
d_i=\sum\limits_{j=1}^nw_{ij}
di=j=1∑nwij
其实就等于邻接矩阵的第
i
i
i行求和
那我们又该如何表示所有的节点的度呢?我们用
D
D
D表示
D
=
[
d
1
0
0
⋯
0
d
2
0
⋯
⋮
⋮
⋱
⋮
0
0
⋯
d
n
]
n
×
n
D=\begin{bmatrix} d_1 & 0 & 0 & \cdots \\ 0 & d_2 & 0 & \cdots \\ \vdots & \vdots &\ddots & \vdots\\ 0 & 0 & \cdots & d_n \end{bmatrix}_{n\times n}
D=
d10⋮00d2⋮000⋱⋯⋯⋯⋮dn
n×n
其实就是一个对角矩阵(主对角线有值,其余都为0)。所以每个样本的都对应矩阵
D
D
D中的一个度。
2.4、核函数
我们如何去构造邻接矩阵 W W W,构造点与点之间的联系?
我们仅介绍一种较为普遍的方法
\boxed{我们仅介绍一种较为普遍的方法}
我们仅介绍一种较为普遍的方法——全连接法,即两个点中间的关系表示为
w
i
j
=
exp
{
−
∣
∣
x
i
−
x
j
∣
∣
2
2
2
σ
2
}
w_{ij}=\exp\left\{-\frac{||x_i-x_j||_2^2}{2\sigma^2}\right\}
wij=exp{−2σ2∣∣xi−xj∣∣22}
分子处为L2范数,其实这就是高斯核函数的定义。同样的,你也可以选择其他核函数。
通过此法,我们便可以构造出邻接矩阵 W W W,由 W W W可得出矩阵 D D D。
2.5、拉普拉斯矩阵
定义 L = D − W L=D-W L=D−W。由这个算出来的矩阵 L L L被称为拉普拉斯矩阵。
拉普拉斯矩阵的性质:
①拉普拉斯矩阵是半正定矩阵。
②特征值中0出现的次数就是图连通区域的个数。
③最小特征值是0,因为拉普拉斯矩阵每一行的和均为0。
并且,拉普拉斯矩阵还有一个性质—— 对任意一个 n 维 y 向量,有 y T L y = 1 2 ∑ i = 1 n ∑ j = 1 n w i j ( y i − y j ) 2 \boxed{对任意一个n维y向量,有y^TLy=\frac{1}{2}\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{ij}(y_i-y_j)^2} 对任意一个n维y向量,有yTLy=21i=1∑nj=1∑nwij(yi−yj)2
证明:
y
T
L
y
=
y
T
D
y
−
y
T
W
y
=
∑
i
=
1
n
d
i
y
i
2
−
∑
i
=
1
n
∑
j
=
1
n
w
i
j
y
i
y
j
=
1
2
(
∑
i
=
1
n
d
i
y
i
2
−
2
∑
i
=
1
n
∑
j
=
1
n
w
i
j
y
i
y
j
+
∑
j
=
1
n
d
j
y
j
2
)
=
1
2
(
∑
i
=
1
n
∑
j
=
1
n
w
i
j
y
i
2
−
2
∑
i
=
1
n
∑
j
=
1
n
w
i
j
y
i
y
j
+
∑
j
=
1
n
∑
i
=
1
n
w
j
i
y
j
2
)
=
1
2
∑
i
=
1
n
∑
j
=
1
n
w
i
j
(
y
i
−
y
j
)
2
\begin{align} y^TLy=&y^TDy-y^TWy\tag{1} \\=&\sum\limits_{i=1}^nd_iy_i^2-\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{ij}y_iy_j\tag{2} \\=&\frac{1}{2}\left( \sum\limits_{i=1}^nd_iy_i^2-2\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{ij}y_iy_j+\sum\limits_{j=1}^nd_jy_j^2 \right)\tag{3} \\=&\frac{1}{2}\left(\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{_{ij}}y_i^2-2\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{ij}y_iy_j+\sum\limits_{j=1}^n\sum\limits_{i=1}^nw_{_{ji}}y_j^2\tag{4} \right) \\=&\frac{1}{2}\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{ij}(y_i-y_j)^2\tag{5} \end{align}
yTLy=====yTDy−yTWyi=1∑ndiyi2−i=1∑nj=1∑nwijyiyj21(i=1∑ndiyi2−2i=1∑nj=1∑nwijyiyj+j=1∑ndjyj2)21(i=1∑nj=1∑nwijyi2−2i=1∑nj=1∑nwijyiyj+j=1∑ni=1∑nwjiyj2)21i=1∑nj=1∑nwij(yi−yj)2(1)(2)(3)(4)(5)
式(3)到式(4)用到
d
i
=
∑
j
=
1
n
w
i
j
d_i=\sum\limits_{j=1}^nw_{ij}
di=j=1∑nwij。式(4)到式(5)用到平方差公式。
三、原理推导
3.1、目标
一个节点就是一个样本。我们要在这张无向图中切一刀,然后得到两部分。比如上图中,用虚线切开,得到两部分。我们就认为上半部分为一类。下半部分为一类。
显然,如何砍这一刀是至关重要的。我们定义砍一刀的代价
目标是砍一刀的代价最小
\boxed{目标是砍一刀的代价最小}
目标是砍一刀的代价最小
W
(
A
,
B
)
=
∑
i
∈
A
∑
j
∉
B
w
i
j
W(A,B)=\sum\limits_{i\in A}\sum\limits_{j\notin B}w_{ij}
W(A,B)=i∈A∑j∈/B∑wij
比如按照上面的无向图,上半部分记为A集合,下半部分记为B集合。那么代价就是a节点与B集合所有点的w(仅有w2,其余为0)加上b节点与B集合所有点的w(仅有w3,其余为0)。所以代价就是
w
2
+
w
3
w_2+w_3
w2+w3。
以此类推,如果要聚成k类,定义代价函数(
1
2
\frac{1}{2}
21是因为不同类别之间会重复计算代价,故如此)
C
u
t
(
A
1
,
A
2
,
⋯
,
A
k
)
=
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
Cut(A_1,A_2,\cdots,A_k)=\frac{1}{2}\sum\limits_{i=1}^kW(A_i,\bar{A_i})
Cut(A1,A2,⋯,Ak)=21i=1∑kW(Ai,Aiˉ)
A
i
ˉ
\bar{A_i}
Aiˉ表示不属于第i类的节点集合。
然而,这种形式的代价有问题,求解的时候算法会趋向于将权重w小的单个节点作为一类,这是我们很不想要的。为此我们必须做出整改。
有两种方案: ① R a t i o C u t \boxed{①RatioCut} ①RatioCut、 N C u t \boxed{NCut} NCut。
3.2、RatioCut
定义 ∣ A i ∣ |A_i| ∣Ai∣表示 i i i集合的节点个数。比如上图中的上半部分集合A,只有节点 a , b a,b a,b,我们就得到 ∣ A ∣ = 2 |A|=2 ∣A∣=2。
我们将代价函数变成这样
R
a
t
i
o
C
u
t
(
A
1
,
A
2
,
⋯
,
A
k
)
=
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
∣
A
i
∣
RatioCut(A_1,A_2,\cdots,A_k)=\frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar{A_i})}{|A_i|}
RatioCut(A1,A2,⋯,Ak)=21i=1∑k∣Ai∣W(Ai,Aiˉ)
也就是要最小化这个代价函数
arg
min
A
i
R
a
t
i
o
C
u
t
(
A
1
,
A
2
,
⋯
,
A
k
)
=
arg
min
A
i
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
∣
A
i
∣
\arg\min\limits_{A_i} RatioCut(A_1,A_2,\cdots,A_k)=\arg\min\limits_{A_i} \frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar{A_i})}{|A_i|}
argAiminRatioCut(A1,A2,⋯,Ak)=argAimin21i=1∑k∣Ai∣W(Ai,Aiˉ)
这种用集合来表达的形式我们是没办法去求解的。因此,我们换一种表达方式。对于每一个样本,我们都定义一个k维的指示向量。
y
=
{
1
∣
A
i
∣
j
∈
A
i
0
j
∉
A
i
j
∈
{
1
,
2
,
⋯
,
k
}
y=\left\{\begin{matrix} \frac{1}{\sqrt|A_i|} & j\in A_i \\ 0 & j\notin A_i &j\in\{1,2,\cdots,k\} \end{matrix}\right.
y={∣Ai∣10j∈Aij∈/Aij∈{1,2,⋯,k}
意思是如果这个样本属于第
i
i
i类,那么对应位置就有值,否则为0。
比如对于第1个样本,假设它属于第二类:第二个样本属于第三类。则有
y
1
=
(
0
1
∣
A
2
∣
0
⋮
0
)
k
×
1
;
y
2
=
(
0
0
1
∣
A
3
∣
⋮
0
)
k
×
1
;
y_1=\begin{pmatrix} 0 \\ \frac{1}{\sqrt|A_2|} \\ 0 \\ \vdots \\ 0 \end{pmatrix}_{k\times1}; y_2=\begin{pmatrix} 0 \\ 0 \\ \frac{1}{\sqrt|A_3|} \\ \vdots \\ 0 \end{pmatrix}_{k\times1};
y1=
0∣A2∣10⋮0
k×1;y2=
00∣A3∣1⋮0
k×1;
以此类推,所有的样本,用
Y
Y
Y表示
Y
=
(
y
1
T
y
2
T
⋮
y
n
T
)
=
[
y
11
y
12
⋯
y
1
k
y
21
y
22
⋯
y
2
k
⋮
⋮
⋱
⋮
y
n
1
y
n
2
⋯
y
n
k
]
n
×
k
Y=\begin{pmatrix} y_1^T \\ y_2^T \\ \vdots \\y_n^T \end{pmatrix}=\begin{bmatrix} y_{11} & y_{12} & \cdots & y_{1k} \\ y_{21} & y_{22} & \cdots & y_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ y_{n1} & y_{n2} & \cdots & y_{nk} \end{bmatrix}_{n\times k}
Y=
y1Ty2T⋮ynT
=
y11y21⋮yn1y12y22⋮yn2⋯⋯⋱⋯y1ky2k⋮ynk
n×k
所以我们的目标就是
arg
min
Y
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
∣
A
i
∣
\arg\min\limits_{Y} \frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar{A_i})}{|A_i|}
argYmin21i=1∑k∣Ai∣W(Ai,Aiˉ)
同样的,里面的目标函数我们也要换成Y的表达形式
对于矩阵
Y
Y
Y,第
m
m
m列中不为0的部分,就是属于第
m
m
m类的样本。我们用
y
m
y_{m}
ym表示
Y
Y
Y的第
m
m
m列
y
m
T
L
y
m
=
1
2
∑
i
=
1
n
∑
j
=
1
n
w
i
j
(
y
m
i
−
y
m
j
)
2
=
1
2
(
∑
i
∈
A
m
∑
j
∉
A
m
w
i
j
(
1
∣
A
m
∣
−
0
)
2
+
∑
i
∉
A
m
∑
j
∈
A
m
w
i
j
(
0
−
1
∣
A
m
∣
)
2
+
∑
i
∈
A
m
∑
j
∈
A
m
w
i
j
(
0
−
0
)
2
)
=
1
2
(
∑
i
∈
A
m
∑
j
∉
A
m
w
i
j
1
∣
A
m
∣
+
∑
i
∉
A
m
∑
j
∈
A
m
w
i
j
1
∣
A
m
∣
)
=
1
2
(
C
u
t
(
A
m
,
A
ˉ
m
)
∣
A
m
∣
+
C
u
t
(
A
m
,
A
ˉ
m
)
∣
A
m
∣
)
=
C
u
t
(
A
m
,
A
ˉ
m
)
∣
A
m
∣
\begin{align} y_{m}^TLy_{m}=&\frac{1}{2}\sum\limits_{i=1}^n\sum\limits_{j=1}^nw_{ij}(y_{mi}-y_{mj})^2\tag{1} \\=&\frac{1}{2}\left( \sum\limits_{i \in A_{m}}\sum\limits_{j \notin A_{m}}w_{ij}(\frac{1}{\sqrt{|A_{m}|}}-0)^2+\sum\limits_{i \notin A_{m}}\sum\limits_{j \in A_{m}}w_{ij}(0-\frac{1}{\sqrt{|A_{m}|}})^2+\sum\limits_{i \in A_{m}}\sum\limits_{j \in A_{m}}w_{ij}(0-0)^2 \right)\tag{2} \\=&\frac{1}{2}\left( \sum\limits_{i \in A_{m}}\sum\limits_{j \notin A_{m}}w_{ij}\frac{1}{|A_{m}|}+\sum\limits_{i \notin A_{m}}\sum\limits_{j \in A_{m}}w_{ij}\frac{1}{|A_{m}|} \right)\tag{3} \\=&\frac{1}{2}\left( \frac{Cut(A_m,\bar{A}_m)}{|A_m|}+\frac{Cut(A_m,\bar{A}_m)}{|A_m|} \right)\tag{4} \\=&\frac{Cut(A_m,\bar{A}_m)}{|A_m|}\nonumber \end{align}
ymTLym=====21i=1∑nj=1∑nwij(ymi−ymj)221
i∈Am∑j∈/Am∑wij(∣Am∣1−0)2+i∈/Am∑j∈Am∑wij(0−∣Am∣1)2+i∈Am∑j∈Am∑wij(0−0)2
21
i∈Am∑j∈/Am∑wij∣Am∣1+i∈/Am∑j∈Am∑wij∣Am∣1
21(∣Am∣Cut(Am,Aˉm)+∣Am∣Cut(Am,Aˉm))∣Am∣Cut(Am,Aˉm)(1)(2)(3)(4)
式(1)用到了上面拉普拉斯矩阵的性质。式(2)用到了指示向量。式(3)到式(4)用到了代价函数
发现了吗,它刚好等于我们第 m m m类的代价。
我们再将原目标函数转化成矩阵的形式
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
∣
A
i
∣
=
∑
i
=
1
k
C
u
t
(
A
i
,
A
ˉ
i
)
∣
A
i
∣
=
T
r
[
C
u
t
(
A
1
,
A
1
ˉ
)
∣
A
1
∣
0
⋯
0
0
C
u
t
(
A
2
,
A
2
ˉ
)
∣
A
2
∣
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
C
u
t
(
A
k
,
A
k
ˉ
)
∣
A
k
∣
]
k
×
k
=
T
r
[
y
1
T
L
y
1
0
⋯
0
0
y
2
T
L
y
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
y
k
T
L
y
k
]
k
×
k
(5)
\begin{align} \frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar{A_i})}{|A_i|} =&\sum\limits_{i=1}^k\frac{Cut(A_i,\bar A_i)}{|A_i|} \\=&\mathbb{Tr}\begin{bmatrix} \frac{Cut(A_1,\bar{A_1})}{|A_1|} & 0 & \cdots & 0\\ 0 & \frac{Cut(A_2,\bar{A_2})}{|A_2|} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \frac{Cut(A_k,\bar{A_k})}{|A_k|} \end{bmatrix}_{k\times k} \\=&\mathbb{Tr}\begin{bmatrix} y_{1}^TLy_{1} & 0 & \cdots & 0\\ 0 & y_{2}^TLy_{2} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & y_{k}^TLy_{k} \end{bmatrix}_{k\times k} \end{align}\tag{5}
21i=1∑k∣Ai∣W(Ai,Aiˉ)===i=1∑k∣Ai∣Cut(Ai,Aˉi)Tr
∣A1∣Cut(A1,A1ˉ)0⋮00∣A2∣Cut(A2,A2ˉ)⋮0⋯⋯⋱⋯00⋮∣Ak∣Cut(Ak,Akˉ)
k×kTr
y1TLy10⋮00y2TLy2⋮0⋯⋯⋱⋯00⋮ykTLyk
k×k(5)
T
r
\mathbb{Tr}
Tr表示求矩阵的迹(也就是对角线求和)。里面是一个对角矩阵。
又因为(对这里不熟悉的,读者可自行尝试,在这里就不展开了)
Y
T
L
Y
=
[
y
1
T
L
y
1
y
1
T
L
y
2
⋯
y
1
T
L
y
k
y
2
T
L
y
1
y
2
T
L
y
2
⋯
y
2
T
L
y
k
⋮
⋮
⋱
⋮
y
k
T
L
y
1
y
k
T
L
y
2
⋯
y
k
T
L
y
k
]
k
×
k
(6)
\begin{align} Y^TLY=&\begin{bmatrix} y_1^TLy_1 & y_1^TLy_2 & \cdots & y_1^TLy_k\\ y_2^TLy_1 & y_2^TLy_2 & \cdots & y_2^TLy_k \\ \vdots & \vdots & \ddots & \vdots \\ y_k^TLy_1 & y_k^TLy_2 & \cdots &y_k^TLy_k \end{bmatrix}_{k\times k} \end{align}\tag{6}
YTLY=
y1TLy1y2TLy1⋮ykTLy1y1TLy2y2TLy2⋮ykTLy2⋯⋯⋱⋯y1TLyky2TLyk⋮ykTLyk
k×k(6)
所以,很容易看到,式(5)和式(6)的对角线是一致的,
那么我们就得到目标函数
:
\boxed{那么我们就得到目标函数:}
那么我们就得到目标函数:
arg
min
Y
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
∣
A
i
∣
=
arg
min
Y
T
r
(
Y
T
L
Y
)
\begin{align} \arg\min\limits_{Y}\frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar{A_i})}{|A_i|}=\arg\min\limits_{Y}\mathbb{Tr}(Y^TLY)\nonumber \end{align}
argYmin21i=1∑k∣Ai∣W(Ai,Aiˉ)=argYminTr(YTLY)
3.3、求解
有了目标表达式,接下来就是求解,但是还有个问题,就是矩阵
Y
Y
Y的每一行是一个指示向量,由这个可得
Y
T
Y
=
I
Y^TY=\mathbb{I}
YTY=I
I
\mathbb{I}
I是单位矩阵,比如我以开头的图构造矩阵
Y
Y
Y,我们切的那刀分成了两份,可得
Y
T
Y
=
[
1
2
1
2
0
0
0
0
0
1
3
1
3
1
3
]
[
1
2
0
1
2
0
0
1
3
0
1
3
0
1
3
]
=
[
1
0
0
1
]
Y^TY=\begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} \end{bmatrix}\begin{bmatrix} \frac{1}{\sqrt{2}} & 0 \\ \frac{1}{\sqrt{2}} & 0 \\ 0 & \frac{1}{\sqrt{3}} \\ 0 & \frac{1}{\sqrt{3}} \\0 & \frac{1}{\sqrt{3}} \end{bmatrix}=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
YTY=[210210031031031]
212100000313131
=[1001]
并且,因为指示向量是离散的。我们是没办法去求导,一般情况下,我们会忽略掉这个指示向量这个约束,而保留
Y
T
Y
=
I
Y^TY=\mathbb{I}
YTY=I
所以我们的问题就变成了
\boxed{所以我们的问题就变成了}
所以我们的问题就变成了
arg
min
Y
T
r
(
Y
T
L
Y
)
s
.
t
.
Y
T
Y
=
I
\begin{aligned} \arg\min\limits_{Y}\mathbb{Tr}&(Y^TLY) \\s.t. \hspace{1cm}& Y^TY=\mathbb{I} \end{aligned}
argYminTrs.t.(YTLY)YTY=I
3.3.1、拉格朗日乘数法
由于 Y Y Y难以求解,所以我们对 Y Y Y的每一列进行求导求解,又因为 T r ( Y T L Y ) = ∑ i k y i T L y i \mathbb{Tr}(Y^TLY)=\sum\limits_{i}^ky_i^TLy_i Tr(YTLY)=i∑kyiTLyi,所以只要每一个最小,那么最后求和也就是最小。但是这种情况其实仅是一个近似解罢了,实际上并不一样,只是由于 Y Y Y太难求了,所以退而求其次。
所以最终的目标函数为
arg
min
y
i
(
y
i
T
L
y
i
)
s
.
t
.
y
i
T
y
i
=
1
\begin{aligned} \arg\min\limits_{y_i}&(y_i^TLy_i) \\s.t. \hspace{1cm}& y_i^Ty_i=1 \end{aligned}
argyimins.t.(yiTLyi)yiTyi=1
约束条件由
Y
T
Y
=
I
→
y
i
T
y
i
=
1
Y^TY=\mathbb{I} \rightarrow y_i^Ty_i=1
YTY=I→yiTyi=1得到
构造拉格朗日函数
L
(
y
i
,
λ
)
=
y
i
T
L
y
i
+
λ
(
1
−
y
i
T
y
i
)
L(y_i,\lambda)=y_i^TLy_i+\lambda(1-y_i^Ty_i)
L(yi,λ)=yiTLyi+λ(1−yiTyi)
在这里给出两个求导公式,不懂的可直接背,对其原理感兴趣自行百度查阅
d
x
T
A
x
x
=
2
A
x
,其中
A
为对称阵,
x
为列向量
d
x
T
x
x
=
2
x
,
x
为列向量
\boxed{\frac{dx^TAx}{x}=2Ax,其中A为对称阵,x为列向量}\\ \boxed{\frac{dx^Tx}{x}=2x,x为列向量}
xdxTAx=2Ax,其中A为对称阵,x为列向量xdxTx=2x,x为列向量
对
y
i
y_i
yi求导
∂
L
(
y
i
,
λ
)
∂
y
i
=
2
L
y
i
−
2
λ
y
i
=
L
y
i
−
λ
y
i
\begin{aligned} \frac{\partial L(y_i,\lambda)}{\partial y_i}=&2Ly_i-2\lambda y_i \\=&Ly_i-\lambda y_i \end{aligned}
∂yi∂L(yi,λ)==2Lyi−2λyiLyi−λyi
令其为0得
L
y
i
=
λ
y
i
Ly_i=\lambda y_i
Lyi=λyi
如果你对特征值和特征向量敏感,便一眼可以看出来,此时的 λ \lambda λ就是特征值,而 y i y_i yi就是特征向量。
式子左右左乘
y
i
T
y_i^T
yiT得
y
i
T
L
y
i
=
λ
y
i
T
y
i
=
λ
y_i^TLy_i=\lambda y_i^Ty_i=\lambda
yiTLyi=λyiTyi=λ
左边不就是我们的目标函数吗?要最小化它,不就是最小化
λ
\lambda
λ?而
λ
\lambda
λ又是拉普拉斯矩阵的特征值。
所以,我们是否可以这么想,对于 Y Y Y的每一列,我们都用拉普拉斯矩阵的特征向量组成,而特征向量的组成便是由前k个特征值(升序)所对应的特征向量。
3.3.2、重离散化
前面我们提到,我们是忽略了指示向量那个离散约束,使得y的取值为全体实数。所以最后我们得出来的 Y Y Y值它也是在全体实数当中的,难以判定它属于哪一类。
所以,我们可以直接在矩阵 Y Y Y当中按行作为样本,使用传统的聚类算法(比如k-mean),实现最终的聚类。
3.4、NCut
弄懂了RatioCut,NCut其实也懂了一半了。
同样的,我们定义
v
o
l
(
A
i
)
=
∑
i
∈
A
i
d
i
vol(A_i)=\sum\limits_{i \in A_i}d_i
vol(Ai)=i∈Ai∑di
也就是对于类别为
A
i
A_i
Ai的集合,集合内所有的节点的度的总和。
而我们的代价函数则定义为
N
C
u
t
(
A
1
,
A
2
,
⋯
,
A
k
)
=
1
2
∑
i
=
1
k
W
(
A
i
,
A
i
ˉ
)
v
o
l
(
A
i
)
=
∑
i
=
1
k
C
u
t
(
A
i
,
A
i
ˉ
)
v
o
l
(
A
i
)
NCut(A_1,A_2,\cdots,A_k)=\frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar{A_i})}{vol(A_i)}=\sum\limits_{i=1}^k\frac{Cut(A_i,\bar{A_i})}{vol(A_i)}
NCut(A1,A2,⋯,Ak)=21i=1∑kvol(Ai)W(Ai,Aiˉ)=i=1∑kvol(Ai)Cut(Ai,Aiˉ)
对于我们每一个样本的指示向量,定义为
y
=
{
1
v
o
l
(
A
i
)
j
∈
A
i
0
j
∉
A
i
j
∈
{
1
,
2
,
⋯
,
k
}
y=\left\{\begin{matrix} \frac{1}{\sqrt{vol(A_i)}} & j\in A_i \\ 0 & j\notin A_i &j\in\{1,2,\cdots,k\} \end{matrix}\right.
y={vol(Ai)10j∈Aij∈/Aij∈{1,2,⋯,k}
同样有(推导过程与RatioCut同,不再赘述)
y
i
T
L
y
i
=
C
u
t
(
A
i
,
A
i
ˉ
)
v
o
l
(
A
i
)
y_i^TLy_i=\frac{Cut(A_i,\bar{A_i})}{vol(A_i)}
yiTLyi=vol(Ai)Cut(Ai,Aiˉ)
所以,同样的都有
N
C
u
t
(
A
1
,
A
2
,
⋯
,
A
k
)
=
T
r
(
Y
T
L
Y
)
NCut(A_1,A_2,\cdots,A_k)=\mathbb{Tr}(Y^TLY)
NCut(A1,A2,⋯,Ak)=Tr(YTLY)
但是,此时
Y
T
Y
≠
I
Y^TY\neq\mathbb{I}
YTY=I,而是
Y
T
D
Y
=
I
Y^TDY=\mathbb{I}
YTDY=I
证明
y
i
仍然代表
Y
的第
i
列
\boxed{y_i仍然代表Y的第i列}
yi仍然代表Y的第i列:
y
i
T
D
y
i
=
∑
j
=
1
n
y
i
j
2
d
j
=
∑
j
∈
A
i
1
v
o
l
(
A
i
)
d
j
=
1
v
o
l
(
A
i
)
∑
j
∈
A
i
d
i
=
1
v
o
l
(
A
i
)
v
o
l
(
A
i
)
=
1
\begin{align} y_i^TDy_i=&\sum\limits_{j=1}^ny_{ij}^2d_j\tag{1} \\=&\sum\limits_{j\in A_i}\frac{1}{vol(A_i)}d_j\tag{2} \\=&\frac{1}{vol(A_i)}\sum\limits_{j\in A_i}d_i\tag{3} \\=&\frac{1}{vol(A_i)}vol(A_i)\tag{4} \\=&1\tag{5} \end{align}
yiTDyi=====j=1∑nyij2djj∈Ai∑vol(Ai)1djvol(Ai)1j∈Ai∑divol(Ai)1vol(Ai)1(1)(2)(3)(4)(5)
式子(1)将矩阵展开可得,式子(2)利用指示向量,式子(4)利用了前面定义的
v
o
l
(
A
i
)
vol(A_i)
vol(Ai)。
因为 y i T D y i y_i^TDy_i yiTDyi和指示向量的定义,可得 Y T D Y = I Y^TDY=\mathbb{I} YTDY=I。
所以,目标函数就变成了
arg
min
Y
T
r
(
Y
T
L
Y
)
s
.
t
.
Y
T
D
Y
=
I
\begin{aligned} &\arg\min\limits_{Y}\mathbb{Tr}(Y^TLY) \\&s.t.\hspace{1cm}Y^TDY=\mathbb{I} \end{aligned}
argYminTr(YTLY)s.t.YTDY=I
由于
y
i
T
y
i
≠
1
y_i^Ty_i\neq1
yiTyi=1,所以我们没办法像RatioCut那样去解决。所以我们不妨对这个式子进行一下转化。
令
Y
=
D
−
1
2
F
Y=D^{-\frac{1}{2}}F
Y=D−21F,所以
arg
min
F
(
Y
T
L
Y
)
=
arg
min
F
T
r
(
F
T
D
−
1
2
L
D
−
1
2
F
)
s
.
t
.
F
T
F
=
I
\begin{aligned} \arg\min_{F} (Y^TLY)=&\arg\min_{F}\mathbb{Tr}(F^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}F) \\&s.t. \hspace{1cm} F^TF=\mathbb{I} \end{aligned}
argFmin(YTLY)=argFminTr(FTD−21LD−21F)s.t.FTF=I
所以,我们只需要将
D
−
1
2
L
D
−
1
2
D^{-\frac{1}{2}}LD^{-\frac{1}{2}}
D−21LD−21当作原来的
L
L
L即可,那么这样一切便都是照常即可。
同样我们一样一列一列地看(假设某一列为f)
arg
min
f
(
f
T
D
−
1
2
L
D
−
1
2
f
)
s
.
t
.
f
T
f
=
1
\begin{aligned} &\arg\min_{f}(f^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}f) \\&s.t. \hspace{1cm} f^Tf=1 \end{aligned}
argfmin(fTD−21LD−21f)s.t.fTf=1
构造拉格朗日函数
L
(
f
,
λ
)
=
f
T
D
−
1
2
L
D
−
1
2
f
+
λ
(
1
−
f
T
f
)
L(f,\lambda)=f^TD^{-\frac{1}{2}}LD^{-\frac{1}{2}}f+\lambda(1-f^Tf)
L(f,λ)=fTD−21LD−21f+λ(1−fTf)
对f求导
∂
L
(
f
,
λ
)
∂
f
=
2
(
D
−
1
2
L
D
−
1
2
)
f
−
2
λ
f
\begin{aligned} \frac{\partial L(f,\lambda)}{\partial f}=2(D^{-\frac{1}{2}}LD^{-\frac{1}{2}})f-2\lambda f \end{aligned}
∂f∂L(f,λ)=2(D−21LD−21)f−2λf
令其等于0,得
(
D
−
1
2
L
D
−
1
2
)
f
=
λ
f
(D^{-\frac{1}{2}}LD^{-\frac{1}{2}})f=\lambda f
(D−21LD−21)f=λf
所以,一样的,我们只需要求出
D
−
1
2
L
D
−
1
2
D^{-\frac{1}{2}}LD^{-\frac{1}{2}}
D−21LD−21的前k个最小值。然后组成矩阵F。不同的是,
N
C
u
t
需要将得到的
F
按行按行归一化
\boxed{NCut需要将得到的F按行按行归一化}
NCut需要将得到的F按行按行归一化,再用传统聚类方法聚类即可。
3.4.1、标准化(一点小补充)
其实, D − 1 2 L D − 1 2 D^{-\frac{1}{2}}LD^{-\frac{1}{2}} D−21LD−21也被叫做对拉普拉斯矩阵L的归一化。同样的还有其他的归一化方法,感兴趣可自行查阅。
4、流程
以RatioCut为例:假设保留前 k 1 k_1 k1个特征,聚成 k 2 k_2 k2类
①求出邻接矩阵w,度矩阵D
②计算拉普拉斯矩阵L
③求出L的特征值,并取出前 k 1 k_1 k1个特征值所对应的特征向量组成Y。
④以Y作为样本,使用传统的聚类方法进行聚类,比如K_mean聚类,聚成 k 2 k_2 k2类
5、结束
以上,便是谱聚类的原理的简单推导了,由于笔者能力有限,还有很多的细节都并未深入解读。如有问题,还望指出,阿里嘎多。