从谱聚类说起

谱聚类(spectral clustering)是一种针对图结构的聚类方法,它跟其他聚类算法的区别在于,他将每个点都看作是一个图结构上的点,所以,判断两个点是否属于同一类的依据就是,两个点在图结构上是否有边相连,可以是直接相连也可以是间接相连。举个例子,一个紧凑的子图(如完全图)一定比一个松散的子图更容易聚成一类。
图卷积神经网络(Graph Convolutional Network, GCN)
那谱聚类为什么叫谱而不是图聚类呢?这个spectral是什么东西?我们知道一个图是可以用一个邻接矩阵A来表示的。而矩阵的谱(spectral)就是指矩阵的特征值,那么这个特征值跟图的矩阵到底有什么深刻的联系呢?
那么首先,图的聚类是什么?我们可以将聚类问题简化为一个分割问题,如果图的结点被分割成A,B这两个集合,那么我们自然是希望在集合中的结点的相互连接更加紧密比如团,而使得子图之间更加尽可能松散。
为了建立这个联系,我们构造一个laplace matrix:

L=DA L=D-A

D是一个对角矩阵,每个对角元素Dii\displaystyle D_{ii}表示第i个结点的度。A则是这个图邻接矩阵。为什么要这样去构造一个矩阵呢?因为研究图的一些性质的时候,我们常常用到一个类似于下式的目标函数:

xTMx={u,v}E(xuxv)2 \mathbf{x^{T}} M\mathbf{x} =\sum _{\{u,v\} \in E}( x_{u} -x_{v})^{2}

这个目标函数可以定义图上的很多问题,比如最小图分割问题,就是要找到一个方法将图分成两块的使得切割的边最少(如果边有权重那就是切割的权重最小)。如下图,你不能找到一个比切两条边更少的分割方法了。
图卷积神经网络(Graph Convolutional Network, GCN)

而这个优化问题其实等价于当x{0,1}V\displaystyle x\in \{0,1\}^{V}的时候:

min{u,v}E(xuxv)2=uA,v̸A(10)2+uA,v̸A(01)2=2cut(A,A) \min\sum _{\{u,v\} \in E}( x_{u} -x_{v})^{2} =\sum _{u\in A,v\not{\in }\overline{A}}( 1-0)^{2} +\sum _{u\in A,v\not{\in }\overline{A}}( 0-1)^{2} =2cut\left( A,\overline{A}\right)

而这个方程不正是一个二次型吗。为了让二次型得到这个结果。我们发现,当M=dIA=DA\displaystyle M=dI-A=D-A的时候就可以了。验证一下:

xT(dIA)x=dxTxxTAx=vdxv22{u,v}Exuxv={u,v}E(xuxv)2 \mathbf{x}^{T}( dI-A)\mathbf{x} =d\mathbf{x}^{T}\mathbf{x} -\mathbf{x}^{T} A\mathbf{x} =\sum _{v} dx^{2}_{v} -2\sum _{\{u,v\} \in E} x_{u} x_{v} =\sum _{\{u,v\} \in E}( x_{u} -x_{v})^{2}

此外,当A不是邻接矩阵而是权重矩阵W的时候,于是d就从度推广到权重的求和,那么这个公式还可以推广为:

xTLx=xT(dIW)x=i,jωij(xuxv)2 x^{T} Lx\mathbf{=x}^{T}( dI-W)\mathbf{x} =\sum _{i,j} \omega _{ij}( x_{u} -x_{v})^{2}

RatioCut 切图聚类

现在,我们可以尝试将这个目标函数与Ratio切图聚类的目标函数建立起联系,建立联系有什么好处呢?好处就是如果我们发现切图的目标函数是这个二次型,那么我们只要优化这个二次型,不就可以用连续的方法来解决一个离散的问题吗?
RatioCut考虑最小化cut(A1,A2,...,Ak)\displaystyle cut( A_{1} ,A_{2} ,...,A_{k}),同时最大化每个子图的个数即:

RatioCut(A1,A2,...Ak)=12i=1kcut(Ai,Ai)Ai RatioCut(A_{1} ,A_{2} ,...A_{k} )=\frac{1}{2}\sum\limits ^{k}_{i=1}\frac{cut(A_{i} ,\overline{A}_{i} )}{|A_{i} |}

其中cut(Ai,Ai)\displaystyle cut(A_{i} ,\overline{A}_{i} )表示两个子图之间的距离(两个子图结点之间距离的求和):

cut(Ai,Ai)=iA,jAiwij cut(A_{i} ,\overline{A}_{i} )=\sum\limits _{i\in A,j\in \overline{A}_{i}} w_{ij}

这里公式里的是A与A的补集的切图权重(切的边权重的求和),也就是说我们希望子图A与其余的图分离的代价最小,比如我只要切掉一条微不足道的边就能将两个复杂的图(比如两个团)分离开,那么就可以认为这是一个好的切割。

现在我们仿照上面的x,将其推广到多个簇,于是我们用一个指示函数(one-hot)来表达每个结点属于哪个子图,这样就将切割问题跟二次型建立起了联系

我们引入指示向量hj{h1,h2,..hk} j=1,2,...kh_{j} \in \{h_{1} ,h_{2} ,..h_{k} \}\ j=1,2,...k,表示有k个子图,对于任意一个向量hj\displaystyle h_{j}, 它是一个|V|-维向量(|V|为结点数,用来标记哪个结点属于哪个子图,类似于one-hot),我们定义hij\displaystyle h_{ij}为:
hij={0vjAi1AivjAi h_{ij} =\begin{cases} 0 & v_{j} \notin A_{i}\\ \frac{1}{\sqrt{|A_{i} |}} & v_{j} \in A_{i} \end{cases}

那么,对于每一个子图都有:
hiTLhi=12m=1Vn=1Vwmn(himhin)2=12(mAi,nAiwmn(1Ai0)2+mAi,nAiwmn(01Ai)2=12(mAi,nAiwmn1Ai+mAi,nAiwmn1Ai=12(cut(Ai,Ai)1Ai+cut(Ai,Ai)1Ai)=cut(Ai,Ai)Ai \begin{aligned} h^{T}_{i} Lh_{i} & =\frac{1}{2}\sum\limits ^{|V|}_{m=1}\sum\limits ^{|V|}_{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{aligned}

其原理在于,因为当mAi,nAi\displaystyle m\in A_{i} ,n\notin A_{i}时,因为结点vm\displaystyle v_{m}属于子图i,结点vn\displaystyle v_{n}不属于子图i,于是himhin=1Ai0\displaystyle h_{im} -h_{in} =\frac{1}{\sqrt{|A_{i} |}} -0,同理,当vm\displaystyle v_{m},vn\displaystyle v_{n}都属于子图的时候himhin=1Ai1Ai=0\displaystyle h_{im} -h_{in} =\frac{1}{\sqrt{|A_{i} |}} -\frac{1}{\sqrt{|A_{i} |}} =0

上述是第i个子图的式子,我们将k个子图的h合并成一个H,于是式子变成:

RatioCut(A1,A2,...Ak)=i=1khiTLhi=i=1k(HTLH)ii=tr(HTLH)s.t. hiThi=1, i=1,2,...,k RatioCut(A_{1} ,A_{2} ,...A_{k} )=\sum\limits ^{k}_{i=1} h^{T}_{i} Lh_{i} =\sum\limits ^{k}_{i=1} (H^{T} LH)_{ii} =tr(H^{T} LH)\\ s.t.\ h^{T}_{i} h_{i} =1,\ i=1,2,...,k

也就是说Ratiocut本质上就是在最小化tr(HTLH)\displaystyle tr(H^{T} LH)这个东西。那么怎么优化呢?注意到每个hi\displaystyle h_{i}都是相互正交的,因为一个结点不能同时属于多个类别,因此H是一个正交矩阵,又因为L是一个对称矩阵,那么可以证明,H是L的特征向量的时候,恰好是这个优化问题的解,我们需要要找到那么特征值比较小的特征向量,就可以找到一种代价最小的切割方法。我们可以来证明一下,特征向量恰好是他的极值:

h(hTLhλ(1hTh))=htr(hTLhλ(1hTh))=htr(hTLh)λhtr(hhT)=utr(uuTL)λhtr(hEhTE)=utr(uEuTL)λhtr(hEhTE)=Lu+LTuλ(h+h)=2Lu2λh=0Lu=λh \begin{aligned} \nabla _{h}\left( h^{T} Lh-\lambda \left( 1-h^{T} h\right)\right) & =\nabla _{h} tr\left( h^{T} Lh-\lambda \left( 1-h^{T} h\right)\right)\\ & =\nabla _{h} tr\left( h^{T} Lh\right) -\lambda \nabla _{h} tr\left( hh^{T}\right)\\ & =\nabla _{u} tr(uu^{T} L)-\lambda \nabla _{h} tr(hEh^{T} E)\\ & =\nabla _{u} tr(uEu^{T} L)-\lambda \nabla _{h} tr(hEh^{T} E)\\ & =Lu+L^{T} u-\lambda ( h+h)\\ & =2Lu-2\lambda h\\ & =0\\ & \Longrightarrow Lu=\lambda h \end{aligned}

这里用到了一些最优化求导常用公式技巧,其实这个推导跟PCA是一样的,只不过PCA找的是最大特征值(PCA中L是协方差矩阵,目标是找到一个向量最大化方差),这里是找最小特征值,我们目标是找到一个向量最小化这个二次型矩阵。
最后,通过找到L的最小的k个特征值,可以得到对应的k个特征向量,这k个特征向量组成一个nxk维度的矩阵,即为我们的H。一般需要对H矩阵按行做标准化,即

hij=hij(t=1khit2)1/2 h^{*}_{ij} =\frac{h_{ij}}{(\sum\limits ^{k}_{t=1} h^{2}_{it} )^{1/2}}

由于我们在使用维度规约的时候损失了少量信息,导致得到的优化后的指示向量h对应的H不能完全指示各样本的归属(因为是连续的优化,不可能恰到得到一个one-hot向量),因此一般在得到nxk维度的矩阵H后还需要对每一行进行一次传统的聚类,比如使用K-Means聚类,从而得到一个真正的one-hot指示向量。

所以谱聚类的流程可以总结如下:

  1. 计算标准化后的lapace矩阵
  2. 求解标准化lapace矩阵的特征值与特征向量
  3. 取最小的k1个特征向量
  4. 对这k1个特征向量聚类,聚类数为k2
  5. 得到k2个簇,就是对应k2个划分。

GCN

图卷积神经网络,顾名思义就是在图上使用卷积运算,然而图上的卷积运算是什么东西?为了解决这个问题题,我们可以利用图上的傅里叶变换,再使用卷积定理,这样就可以通过两个傅里叶变换的乘积来表示这个卷积的操作。那么为了介绍图上的傅里叶变换,我接来下从最原始的傅里叶级数开始讲起。

从傅里叶级数到傅里叶变换

此部分主要参考了马同学的两篇文章:

  1. 从傅立叶级数到傅立叶变换
  2. 如何理解傅立叶级数公式?

傅里叶级数的直观意义

如下图,傅里叶级数其实就是用一组sin,cos的函数来逼近一个周期函数,那么每个sin,cos函数就是一组基,这组基上的系数就是频域,你会发现随着频域越来越多(基越来越多),函数的拟合就越准确。
图卷积神经网络(Graph Convolutional Network, GCN)

傅里叶变换推导

要讲傅里叶变换的推导,我们要先从傅里叶级数讲起,考虑一周期等于T,现定义于区间[-T/2,T/2]的周期函数f(x),傅里叶级数近似的表达式如下:

f(x)=C+n=1(ancos(2πnTx)+bnsin(2πnTx)),CR {\displaystyle f(x)=C+\sum ^{\infty }_{n=1}\left( a_{n} cos(\frac{2\pi n}{T} x)+b_{n} sin(\frac{2\pi n}{T} x)\right) ,C\in \mathbb{R}}

利用偶函数*奇函数=奇函数的性质可以计算出ak\displaystyle a_{k}bk\displaystyle b_{k}

an=T/2T/2f(x)cos(2πnTx)dxT/2T/2cos2(2πnTx)dx=2TT/2T/2f(x)cos(2πnTx)dxbn=T/2T/2f(x)sin(2πnTx)dxT/2T/2sin2(2πnTx)dx=2TT/2T/2f(x)sin(2πnTx)dx a_{n} =\frac{\int ^{T/2}_{-T/2} f(x)cos(\frac{2\pi n}{T} x)dx}{\int ^{T/2}_{-T/2} cos^{2} (\frac{2\pi n}{T} x)dx} =\frac{2}{T}\int ^{T/2}_{-T/2} f(x)cos(\frac{2\pi n}{T} x)dx\\ b_{n} =\frac{\int ^{T/2}_{-T/2} f(x)sin(\frac{2\pi n}{T} x)dx}{\int ^{T/2}_{-T/2} sin^{2} (\frac{2\pi n}{T} x)dx} =\frac{2}{T}\int ^{T/2}_{-T/2} f(x)sin(\frac{2\pi n}{T} x)dx

利用欧拉公式eix=cosx+isine^{ix} =\cos x+i\sinx,我们发现cosx,sinx\displaystyle \cos x,\sin x 可表示成

cosx=eix+eix2,sinx=eixeix2i {\displaystyle \cos x=\frac{e^{ix} +e^{-ix}}{2} ,\sin x=\frac{e^{ix} -e^{-ix}}{2i} ,}

再将傅立叶级数f(x)中cos(2πnTx)\cos (\frac{2\pi n}{T} x)sin(2πnTx)\sin (\frac{2\pi n}{T} x)的线性组合式改写如下:

ancos(2πnTx)+bnsin(2πnTx)=an(ei2πnTx+ei2πnTx2)+bk(ei2πnTxei2πnTx2i)=(anibn2)ei2πnTx+(an+ibn2)ei2πnTx=cnei2πnTx+cnei2πnTx \begin{aligned} a_{n}\cos (\frac{2\pi n}{T} x)+b_{n}\sin (\frac{2\pi n}{T} x) & =a_{n}\left(\frac{e^{i\frac{2\pi n}{T} x} +e^{-i\frac{2\pi n}{T} x}}{2}\right)+b_{k}\left(\frac{e^{i\frac{2\pi n}{T} x} -e^{-i\frac{2\pi n}{T} x}}{2i}\right)\\ & =\left(\frac{a_{n} -ib_{n}}{2}\right) e^{i\frac{2\pi n}{T} x} +\left(\frac{a_{n} +ib_{n}}{2}\right) e^{-i\frac{2\pi n}{T} x}\\ & =c_{n} e^{i\frac{2\pi n}{T} x} +c_{-n} e^{-i\frac{2\pi n}{T} x} \end{aligned}

可以验证cn=anibn2=an+ibn2\displaystyle c_{-n} =\frac{a_{-n} -ib_{-n}}{2} =\frac{a_{n} +ib_{n}}{2},这是因为an是一个偶函数,bn是一个奇函数。此外,若n=0,就有c0=a0/2c_{0} =a_{0} /2。将以上结果代回f(x)的傅立叶级数即得指数傅立叶级数:

f(x)=n=cnei2πnxT {\displaystyle f(x)=\sum ^{\infty }_{n=-\infty }\underbrace{c_{n}}_{基的坐标} \cdot \underbrace{e^{i\tfrac{2\pi nx}{T}}}_{正交基}}

现在我们知道cn=anibn2\displaystyle c_{n} =\frac{a_{n} -ib_{n}}{2},将an,bn\displaystyle a_{n} ,b_{n}的结果代进去可以得到:
cn=1TT/2T/2f(x)(cos(2πnTx)isin(2πnTx))dx=1TT/2T/2f(x)ei2πnTxdx c_{n} =\frac{1}{T}\int ^{T/2}_{-T/2} f(x)(\cos (\frac{2\pi n}{T} x)-i\sin (\frac{2\pi n}{T} x))dx=\frac{1}{T}\int ^{T/2}_{-T/2} f(x)e^{-i \frac{2\pi n}{T}} x dx

公式用频率替换:Δω=2πT\displaystyle \Delta \omega =\frac{2\pi }{T},再令ωn=ωn\displaystyle \omega _{n} =\omega n现在我们可以写出全新的傅里叶级数:

f(x)=n=Δω2πT/2T/2f(x)eiωnxdxeiωnx {\displaystyle f(x)=\sum ^{\infty }_{n=-\infty }\frac{\Delta \omega }{2\pi }\int ^{T/2}_{-T/2} f(x)e^{-i\omega _{n} x} dx\cdot } e^{i\omega _{n} x}

现在令TΔω0\displaystyle T\rightarrow \infty ,\Delta \omega \rightarrow 0,并设F(ω)=limTT/2T/2f(x)eiωxdx\displaystyle F{\displaystyle ( \omega ) =\lim _{T\rightarrow \infty }\int ^{T/2}_{-T/2} f(x)e^{-i\omega x} dx}
f(x)=n=Δω2πF(ωn)eiωnx=12πn=F(ωn)eiωnxΔω=12π+F(ω)eiωxdω \begin{aligned} {\displaystyle f(x)} & ={\displaystyle \sum ^{\infty }_{n=-\infty }\frac{\Delta \omega }{2\pi } F( \omega _{n}) \cdot } e^{i\omega _{n} x}\\ & ={\displaystyle \frac{1}{2\pi }\sum ^{\infty }_{n=-\infty } F( \omega _{n}) \cdot } e^{i\omega _{n} x} \Delta \omega \\ & ={\displaystyle \frac{1}{2\pi }\int ^{+\infty }_{-\infty } F( \omega ) \cdot } e^{i\omega x} d\omega \end{aligned}

于是得到了傅里叶变换就是

F(ω)=+f(x)eiωxdx {\displaystyle F( \omega ) =\int ^{+\infty }_{-\infty } f(x)e^{-i\omega x} dx}

Signal Processing on Graph

在将图的傅里叶变换之前,我们先介绍一下图信号是什么。我们在传统概率图中,考虑每个图上的结点都是一个feature,对应数据的每一列,但是图信号不一样,这里每个结点不是随机变量,相反它是一个object。也就是说,他描绘概率图下每个样本之间的图联系,可以理解为刻画了不满足iid假设的一般情形。
图卷积神经网络(Graph Convolutional Network, GCN)

图上的傅里叶变换

那么我们要怎么将传统的傅里叶变换推广到图结构中去?回忆一下,传统对f作傅里叶变换的方法:

f^(ξ):=<f,e2πiξt>=Rf(t)e2πiξtdt \hat{f} (\xi ):=\left< f,e^{2\pi i\xi t}\right> =\int _{\mathbb{R}} f(t)e^{-2\pi i\xi t} dt

我们换了种写法,其实我们发现这个傅里叶变换本质上是一个内积。这个e2πiξt\displaystyle e^{-2\pi i\xi t}其实是lapace算子的一个特征函数,可以理解为一种特殊形式的特征向量:

Δ(e2πiξt)=2t2e2πiξt=(2πξ)2e2πiξt -\Delta \left( e^{2\pi i\xi t}\right) =-\frac{\partial ^{2}}{\partial t^{2}} e^{2\pi i\xi t} =(2\pi \xi )^{2} e^{2\pi i\xi t}

注意,这里导数本质上是一个线性变换,因为它满足线性算子的两个性质,T(x+y)=T(x)+T(y), cT(x)=T(cx)。可以看到e2πiξt\displaystyle e^{2\pi i\xi t}是laplace算子的特征向量,而(2πξ)2\displaystyle (2\pi \xi )^{2}则是lapace算子的特征值。那么在图上我们的laplace矩阵就是离散化的lapace算子,而这个算子在图上的基显然就是特征向量了!

因此,只要意识到传统的傅里叶变换本质上求的是与正交基的内积(比如基e2πiξt\displaystyle e^{2\pi i\xi t})上的系数,而推广到图上的正交基很显然就是laplace矩阵的特征值,于是对于laplace矩阵的傅里叶变换就可以表达为:

f^(λl):=<f,ul>=i=1Nf(i)ul(i) \hat{f}( \lambda _{l}) :=< \mathbf{f} ,\mathbf{u}_{l}> =\sum ^{N}_{i=1} f(i)u^{*}_{l} (i)

显然这个变换就是在求解特征向量的系数,也就是特征值,因此,可以理解为图上的经过傅里叶变换后的函数f^\displaystyle \hat{f}就是一个计算特征值的函数。

更一般的,图上的傅里叶变换可以写成以下内积的形式,其中U是laplace矩阵的特征向量矩阵:
傅里叶变换:

x^=UTx \hat{x} =U^{T} x

傅里叶逆变换:

x=Ux^ x=U\hat{x}

因此,我们就可以定义图上的卷积,因为它就是简单的两个变换的乘积而已:
比如,x,y的卷积,就是他们傅里叶变换相乘

yx=UTyUx y\star x=U^{T} yU^{\top } x

如果我们将y参数化,设gθ=diag(θ)\displaystyle g_{\theta } =\operatorname{diag} (\theta ),就可以训练一个卷积核:

gθx=UgθUx g_{\theta } \star x=Ug_{\theta } U^{\top } x

然而计算U的代价太高了,因此要想办法去近似它,有人提出,

gθ(Λ)k=0KθkTk(Λ~) g_{\theta ^{\prime }} (\Lambda )\approx \sum ^{K}_{k=0} \theta ^{\prime }_{k} T_{k} (\tilde{\Lambda } )
其中Λ~=2λmaxΛIN\displaystyle \tilde{\Lambda } =\frac{2}{\lambda _{\max}} \Lambda -I_{N},现在假设λmax2\displaystyle \lambda _{\max} \approx 2,并且k=1,于是

gθxθ0x+θ1(LIN)x=θ0xθ1D12AD12x g_{\theta ^{\prime }} \star x\approx \theta ^{\prime }_{0} x+\theta ^{\prime }_{1}( L-I_{N}) x=\theta ^{\prime }_{0} x-\theta ^{\prime }_{1} D^{-\frac{1}{2}} AD^{-\frac{1}{2}} x

最后再假设这两个参数是共享的,可以得到:

gθxθ(IN+D12AD12)x g_{\theta } \star x\approx \theta \left( I_{N} +D^{-\frac{1}{2}} AD^{-\frac{1}{2}}\right) x

最后,再将中间的项用一个trick变成:IN+D12AD12D~12A~D~12\displaystyle I_{N} +D^{-\frac{1}{2}} AD^{-\frac{1}{2}}\rightarrow \tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}},其中A~=A+IN\displaystyle \tilde{A} =A+I_{N},D~ii=jA~ij\displaystyle \tilde{D}_{ii} =\sum _{j}\tilde{A}_{ij},最后的最后,终于得到了这样的近似卷积公式:

Z=D~12A~D~12XΘ Z=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}} X\Theta

这样我们就可以直接用神经网络训练了:

H(l+1)=σ(D~12A~D~12H(l)W(l)) H^{(l+1)} =\sigma \left(\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}} H^{(l)} W^{(l)}\right)

参考资料

https://towardsdatascience.com/spectral-clustering-82d3cff3d3b7
https://www.cnblogs.com/pinard/p/6221564.html
https://tkipf.github.io/graph-convolutional-networks/
https://ccjou.wordpress.com/2012/04/03/傅立葉級數-下/
https://www.matongxue.com/madocs/712/
https://www.youtube.com/watch?v=Q99ZPGnUBAQ

相关文章: