原文:Global spectral clustering in dynamic networks
作者:Fuchen Liua(a), David Choib(b), Lu Xie(a), and Kathryn Roedera(a,c)
a:Department of Statistics and Data Science, Carnegie Mellon University, Pittsburgh, PA 15213; b:Heinz College, Carnegie Mellon University, Pittsburgh, PA15213; c: Computational Biology Department, Carnegie Mellon University, Pittsburgh, PA 15213
1. 问题
这篇文章提出了一个新的基于谱聚类的动态社区检测方法。
网络或图可以用表示显示复杂系统内的连接。网络中的顶点往往会表现出聚类(簇)的形式,同一个簇的顶点之间有较多的边相连接,而不同簇之间的顶点连接较为稀疏。这样的集群或社区可能源自网络的不同组成部分,例如基因调控细胞过程。目前社区检测比较成熟的统计(机器学习)理论主要集中在静态网络上(一段时间或一个发展时期的单一快照)。而在现实中,网络通常是动态的,并且对它们的进化进行可视化和建模是非常有意义的(我感觉比研究静态的网络更有意义,因为动态的网络才是现实的真实体现,而现实中少有静止不变的网络)。并且它的应用比比皆是,例如Twitter中的社交网络、物理学中的动态扩散网络和大脑发展的基因共表达网络。在所有这些领域中,为说明网络节点之间的关系的结构以及它们如何随着时间的推移而变化,社区检测是至关重要的。尽管静态网络中的统计推断方法已经很好的建立起来了,但是相对来说研究者还没有充分理解如何在动态网络中结合信息。最近的工作试图将社区检测扩展到动态网络、中心性和动态数据的聚类。
当网络结构具有不确定性时,社区检测是具有挑战性的。动态网络呈现额外的挑战,但也增加了跨时间段的信息。作者提出了一个全局社区检测方法,即持久社区的特征向量平滑方法(PISCES),将信息通过一系列网络,纵向结合,以加强对每个时期的网络结构推断。我们的方法从进化谱聚类和程度校正方法推导出来,并且为调节参数提供了数据驱动的解决方案。在模拟中,我们发现PISCES的性能优于设计用于低信噪比的同类方法。最近从恒河猴脑获得的基因表达数据提供了精细划分脑的样本。广泛的时间跨度包括出生前和出生后的时期。我们感兴趣的是基因群落是如何在空间上发展的。然而,一旦数据被分为均匀的空间和时间周期,样本量非常小,使得推断颇具挑战性。PISCES应用于从近孕到成年的恒河猴大脑皮层的内侧前额叶,显示出持续、融合和发散的、致密群落说明如何动态社区检测可以产生有趣的见解,如脑发育过程。
2. 解决方法
我们的方法,持久社区的特征向量平滑方法(PISCES),实现度校正谱聚类,用平滑项来促进跨时间段的相似性,迭代直到达到收敛为止。具体来说,这个全局谱聚类方法将当前网络与前面和未来结果的主特征向量组合,该组合形成为可以优化的问题。在适度平滑的情况下在全局范围内解决社区的数量已知的问题。我们发现,重要的是选择合适的平滑和模型阶数,
以及平衡规则与“让数据说话”。我们使用数据驱动的方法来实现这一点。
具体的解决方法如下。
...
静态网络的谱方法:
谱聚类是一类较流行的在静态网络中寻找社区的方法,在文献(19—22)中已经讨论了 它的许多变化。原型法由参考文献6给出。给定对称n×n邻接矩阵A和固定社区数K, 该方法计算度归一化邻接矩阵(或“拉普拉斯化”邻接矩阵L):
...
然后,该方法返回由对L的K个绝对值最大的特征根对应的特征向量进行k-均值聚类得 到的聚类结果。选个聚类个数K 的方法见文献23-26.
...
动态网络的特征向量平滑方法:
...
思想:
动态网络中,当r(变化率)不是很大的时候,网络前后的结构肯定是有很强的联系的。 如果某两个结点在一段时间内都有联系,那么把它们聚成一类是合理的。通过结合当前 网络极其前后网络的结构,能更好地利用我们所拥有的信息,得到更有功效的统计推断。
3. 结果
作者进行了模拟和对真实数据进行了分析,并于其他方法进行了比较。在模拟中,能明显看出PISCES要优于其他方法。而在真实数据(恒河猴大脑皮层的内侧前额叶的基因表达样本)中,作者得到非常有说服力的结果。
以下是论文中部分结果的图片。
4. 算法优缺点
优点:提出了一个新的对动态网络进行社区检测的方法。算法简单,容易实现。时间复杂度不高,有理论保证可以收敛到最优解。
缺点:谱聚类要求每个节点的度(或者权重和)都大于0 ,但是在文章中这种用邻接矩阵表示结点之间的联系时,有时候会出现结点的度为0的情况。在一般的谱聚类中,我们(我感觉)可以把这种度为0的结点直接删去(或直接划为一类)。但是在文章遇到的这种情况中,结点个数n和社团个数K是不变的,所以不能直接删去。但是作者并没有说遇到这种情况应该如何处理(而且这种情况很容易遇到,设想某个结点很有可能在某个时期和所有社团都没有联系,但在之后的时期和某个社团产生联系)。此外,真正的动态网络肯定是n和K都会发生变化的,作者限定n和K不变,适用范围太窄,没有很大的实用价值。但是他这种处理方式可以给我们提供一种研究思路。此外算法中还有一个表示参数a要选取,如何选取论文正文没有给出明确解答。
5. 复现
对每一个条件,均做一百次模拟(与论文一样),计算平均的ARI。
可以看出,在大部分条件下(尤其是easy和medium以及r较小时,以及时间段的中间部分)PISCES的结果明显好于静态谱聚类的结果。
1.n = 100, K = 2, pout = 0.1
r = 0, pin = (0.2,0.2) T=1,2,3,4,5 | |||||
Static | 0.4033438 | 0.4059342 | 0.3958336 | 0.4400773 | 0.4123917 |
PISCES | 0.3947865 | 0.4368355 | 0.4298128 | 0.4604113 | 0.3964476 |
r = 0, pin = (0.2,0.25) | |||||
Static | 0.5861688 | 0.5978856 | 0.6078584 | 0.6145001 | 0.5781638 |
PISCES | 0.6516057 | 0.6791103 | 0.6879314 | 0.6932404 | 0.6479088 |
r = 0, pin = (0.25,0.25) | |||||
Static | 0.7682986 | 0.7725866 | 0.7822470 | 0.7725070 | 0.8002579 |
PISCES | 0.8419036 | 0.8657115 | 0.8903100 | 0.8760034 | 0.8596214 |
r = 0.2, pin = (0.2,0.2) T=1,2,3,4,5 | |||||
Static | 0.4055519 | 0.4322142 | 0.4071909 | 0.4237366 | 0.4078626 |
PISCES | 0.4130687 | 0.4375518 | 0.3858658 | 0.4246821 | 0.3976516 |
r = 0.2, pin = (0.2,0.25) | |||||
Static | 0.6142186 | 0.6237895 | 0.6227808 | 0.6418403 | 0.6508497 |
PISCES | 0.6839912 | 0.6938021 | 0.6863698 | 0.7133095 | 0.6880474 |
r = 0.2, pin = (0.25,0.25) | |||||
Static | 0.7471508 | 0.7673495 | 0.7453637 | 0.7617394 | 0.7630894 |
PISCES | 0.8191363 | 0.8425363 | 0.8626521 | 0.8670415 | 0.8292460 |
r = 0.5, pin = (0.2,0.2) T=1,2,3,4,5 | |||||
Static | 0.4107960 | 0.4007530 | 0.4211592 | 0.3862258 | 0.4473993 |
PISCES | 0.3771513 | 0.3973456 | 0.4205767 | 0.3946420 | 0.4292825 |
r = 0.5, pin = (0.2,0.25) | |||||
Static | 0.6190488 | 0.6476665 | 0.6408698 | 0.5925187 | 0.6382355 |
PISCES | 0.6532756 | 0.6888740 | 0.6971124 | 0.6637767 | 0.6425344 |
r = 0.5, pin = (0.25,0.25) | |||||
Static | 0.7439139 | 0.7777161 | 0.7774293 | 0.7743039 | 0.7625964 |
PISCES | 0.8121530 | 0.8282229 | 0.8329703 | 0.8389092 | 0.8302620 |
6. 总结
模拟的时候遇到的问题:由于是随机生成的数据,运气不好的时候会遇到有些节点的度为0的情况,这样矩阵D就不可逆,于是没法使用谱聚类。鉴于在使用谱聚类的时时候,应该把度为0的结点给删掉,所以遇到这种问题,我的解决方法就是加上一个很小的数,使得矩阵可逆。
由于时间有限,我没有像论文那样,同时比较了PisCES和三个benchmark聚类算法的效果,只比较了PisCES和静态(普通)谱聚类算法的效果。不知道是由于代码的差异,我并没有达到论文中的效果,可是也可以得出PisCES的效果比静态谱聚类效果要好的结论。
这是我第一次实现论文中的算法,花了接近六个小时,比看论文的时间还要长。“纸上得来终觉浅,觉知此事要躬行。”看论文的时候,作者的思路、想法都很好理解,证明的部分比较难,但是跟着推导一遍也能懂。而让我觉得收获最大的,还是动手去实现算法。即使是看起来比较简单的算法,实现起来也会遇到各种问题,而这些问题是之前看算法的时候完全没有考虑到的,这加深了我对算法的认识,让我意识算法的一些优缺点以及算法原来没有考虑过的问题。
7. 代码
PisCESdata <-function(Time=5, K=2, n=100, pin=c(0.2,0.25), pout=0.1, r=0){
A <- array(0, c(n,n,Time))
B <- array(0, c(K,K,Time))
piT <- array(0, c(n,Time))
fT <- array(0, c(n,Time))
ZT <- array(0, c(n,Time+1))
ZT[,1]<-sample(1:K, n, replace = TRUE, prob=rep(1/K, K))
for(tt in 1:Time){
B[1:K,1:K,tt] <- rep(pout,K*K)
diag(B[ , ,tt]) <- runif(K,pin[1],pin[2])
piT[,tt] <- sample(1:n,n,prob=rep(1/n,n))
fT[,tt] <- 1/2 + piT[,tt]/n
zt = sample(1:K, n, replace = TRUE, prob=rep(1/K, K))
it = runif(n)
ZT[,tt+1] <-(zt - ZT[,tt])*I(it<r) + ZT[,tt]
S <- (fT[,tt]%*%t(fT[,tt]))*B[ZT[,tt+1],ZT[,tt+1],tt]
for(i in 1:(n-1)){
A[i,i,tt] <-1
for(j in (i+1):n){
A[i,j,tt] <- rbinom(1,1,S[i,j])
A[j,i,tt] <- A[i,j,tt]
}
}
}
return(list(A,ZT[,2:(Time+1)]))
}
PisCES <- function(A, Time, n, K, a, h){
U <- array(0, c(n,n,Time))
V <- array(0, c(n,K,Time))
for(tt in 1:Time){
d <- apply(A[,,tt],1,sum)
if(prod(d)==0){d=d+0.001}
D <- diag(d,n)
L <- solve(sqrt(D))%*%A[,,tt]%*%solve(sqrt(D))
eigenL <- eigen(L)
c <- order(abs(eigenL$values))[(n-K+1):n]
V[,,tt] <- eigenL$vectors[,c]
U[,,tt] <- V[,,tt]%*%t(V[,,tt])
}
U2 <- U
for(k in 1:h){
U1 <- U2
eigenU <- eigen(U[,,1] +a*U1[,,2])
c <- order(abs(eigenU$values))[(n-K+1):n]
v <- eigenU$vectors[,c]
U2[,,1] <- v%*%t(v)
for(tt in 2:(Time-1)){
eigenU <- eigen(a*U1[,,tt-1] + U[,,tt] +a*U1[,,tt+1])
c <- order(abs(eigenU$values))[(n-K+1):n]
v <- eigenU$vectors[,c]
U2[,,tt] <- v%*%t(v)
}
eigenU <- eigen(a*U1[,,Time-1] + U[,,Time])
c <- order(abs(eigenU$values))[(n-K+1):n]
v <- eigenU$vectors[,c]
U2[,,Time] <- v%*%t(v)
}
return(list(U1,U2))
}
Compare <- function(N, Time, n, K, pin=c(0.25,0.25), r=0.5){
ARI1 <- array(0, c(Time,N))
ARI2 <- array(0, c(Time,N))
for(i in 1:N){
PisCESdata1 <- PisCESdata(pin=pin,r=r)
A <- PisCESdata1[[1]]
ZT <- PisCESdata1[[2]]
PisCES1 <- PisCES(A, 5, 100, 2, 0.1, 100)
for(j in 1:Time){
model1 <- specc(A[,,j], 2)
SL1 <- as.factor(model1)
pam <- table(ZT[,j],SL1)
ARI1[j,i] <- randIndex(pam)
eigenU <- eigen(PisCES1[[2]][,,j])
c <- order(abs(eigenU$values))[(n-1):n]
v <- eigenU$vectors[,c]
model2 <- kmeans(v, 2)
SL2 <- as.factor(model2$cluster)
pam <- table(ZT[,j],SL2)
ARI2[j,i] <- randIndex(pam)
}
}
result1 = apply(ARI1, 1, mean)
result2 = apply(ARI2, 1, mean)
print(result1)
print(result2)
return(list(result1,result2))
}
N =100
Time =5
#1. n = 100
n=100
#1.hard
K = 2
pin = c(0.2, 0.2)
r = 0
result1 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.2, 0.25)
r = 0
result2 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.25, 0.25)
r = 0
result3 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.2, 0.2)
r = 0.2
result4 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.2, 0.25)
r = 0.2
result5 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.25, 0.25)
r = 0.2
result6 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.2, 0.2)
r = 0.5
result7 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.2, 0.25)
r = 0.5
result8 = Compare(N, Time, n, K, pin, r)
K = 2
pin = c(0.25, 0.25)
r = 0.5
result9 = Compare(N, Time, n, K, pin, r)
#2. n = 500
K = 10
pin = c(0.2, 0.2)
r = 0
result10 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.2, 0.25)
r = 0
result11 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.25, 0.25)
r = 0
result12 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.2, 0.2)
r = 0.2
result13 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.2, 0.25)
r = 0.2
result14 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.25, 0.25)
r = 0.2
result15 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.2, 0.2)
r = 0.5
result16= Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.2, 0.25)
r = 0.5
result17 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.25, 0.25)
r = 0.5
result18 = Compare(N, Time, n, K, pin, r)
#3. n = 500
pout = 0.3
K = 10
pin = c(0.55, 0.55)
r = 0
result19 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.6, 0.6)
r = 0
result20 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.4, 0.8)
r = 0
result21 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.55, 0.55)
r = 0.2
result22 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.6, 0.6)
r = 0.2
result23 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.4, 0.8)
r = 0.2
result24 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.55, 0.55)
r = 0.5
result25 = Compare(N, Time, n, K, pin, r, pout)
K = 10
pin = c(0.6, 0.6)
r = 0.5
result26 = Compare(N, Time, n, K, pin, r)
K = 10
pin = c(0.4, 0.8)
r = 0.5
result27 = Compare(N,