一、任务简介
通过编程实现两个社区发现算法，并利用实例进行验证，最后对算法的复杂度进行分析。
本次大作业选取的算法为GN算法和Louvain算法
二、概念简介
复杂网络
在网络理论的研究中，复杂网络是由数量巨大的节点和节点之间错综复杂的关系共同构成的网络结构。用数学的语言来说，就是一个有着足够复杂的拓扑结构特征的图。复杂网络具有简单网络，如规则网络、随机图等结构所不具备的特性，而这些特性往往出现在真实世界的网络结构中。
在自然界和人类社会中，复杂网络无处不在，如自然界中的生态网络、新陈代谢网络、蛋白质网络等；人类社会中的社交网络、交通网络、电子邮件网络、科学家合作网、电子电路网络等，这些复杂系统都可以用复杂网络的相关特性进行描述和分析，网络中的结点表示系统个体，边表示个体之间的关系。如今，复杂网络的研究是现今科学研究中的一个热点，它已渗透到各个领域，并成为这些领域的重要研究方向
社区发现
随着复杂网络理论的发展，继1998年Watts和Strogtz提出了WS小世界模型(Small-world network)，1999年Barabasi和Albert提出无标度网络模型(Scale-free network)之后，人们发现复杂网络具有一定的社区结构。相同类型节点之间连接较多，构成一个一个的小社区，不同类型节点之间连接较少，但成为沟通不同社区的重要通道，这种连接的不均匀性表明，网络内部存在一定的自然分化。
对于社区（community），目前没有明确的定义，常见的是Newman和Gievan提出的：社区是图的一个子图，相比于图的其他部分，其中包含着密集的节点，或者等价地说，如果图一个子图内部的链接数量高于这些子图之间的链接数，我们就可以说这个图含有社区。如下图所示，图中网络包含三个社区结构，分别对应图中的三个虚线圈内的节点，同一社区内的节点和节点之间连接很紧密，而社区与社区之间的连接比较稀疏。换而言之，复杂网络的社区结构是复杂网络中社区之间的连接相对稀疏，而社区内部的连接相对稠密的特性。
介数和边介数
2004年，Newman和Girvan的论文中提出了利用介数中心性（betweenness centrality）的方法，称之为Girvan-Newman算法，这篇论文也是开创了社区发现这一领域的文章。本文主要针对无权和有权网络，对此类算法进行一定的研究分析。
在图论中，介数中心性（betweenness centrality）是基于最短路径的图的中心的度量方法。对于连通图中的每对顶点，在顶点之间至少存在一条最短路径，而一个节点的介数一般是指通过该节点的最短路径数。
节点介数中心性的定义如下：
$$
g(v)=\sum_{s \neq v \neq t \in V} \frac{\sigma_{s t(v)}}{\sigma_{s t}}
$$
其中$\sigma_{s t(v)}$表示通过节点$v$的节点$s\Rightarrow t$最短路径的数量，$s、t$表示图中的节点，$\sigma_{st}$表示节点$s\Rightarrow t$的最短路径的数量。直观上来说，节点的介数中心性反映了节点$v$作为”桥梁”的重要程度。
通过节点介数我们可以引出边介数的概念，边介数即为经过该边的最短路径的数量，因此可以引出边介数中心性的概念，它可以用来衡量图中一条边作为”桥梁”的重要程度。
模块度
模块度，是指复杂网络中连接社区内部节点的边所占的比例，与另一个随机网络中连接社区内部节点的边所占的比例的期望值的差值。其中，这个随机网络的构造方法为：保持每个节点的社区属性不变，根据节点度对节点间的边进行随机连接。如果社区结构比较明显，则社区内部连接的边的比例应改高于随机网络的期望值。
模块度定义如下：
$$
Q=\frac{1}{2m}\sum_{i,j}[A_{ij}-\frac{k_ik_j}{2m}]\delta(c_i,c_j)
$$
其中，m表示网络中边的数量，A为邻接矩阵，$A_{ij}$就是i和j的之间的连边权重，$k_{i}$表示所有指向i节点的连接边权重和，$\delta (c_i,c_j)$满足：
社区发现效果的评价指标——NMI指数
NMI是社区发现(community detection)在有标准ground-truth的情况下的重要衡量指标，基本可以比较客观地评价出一个社区划分与标准划分之间相比的准确度。NMI的值域是0到1，越高代表划分得越准。NMI指数的计算公式如下：
$$
NMI(X;Y)=2\frac{I(X;Y)}{H(X)+H(Y)}
$$
其中，$X,Y$表示实验社区发现的结果以及实际情况，$H(X)$表示$X$的信息熵，计算公式如下：
$$
H(X)=-\sum_{i}p(x_i)logp(x_i)
$$
$I(X;Y)$表示$X,Y$的互信息，它可以看成是一个随机变量中包含的关于另一个随机变量的信息量，或者说是一个随机变量由于已知另一个随机变量而减少的不肯定性。计算公式如下：
$$
I(X;Y)=\sum_{x}\sum_{y}p(x,y)log\frac{p(x,y)}{p(x)p(y)}
$$
$p(x),p(x,y)$分别表示随机变量$X,Y$的边缘分布和联合分布。
三、算法原理介绍
基于边介数的社区发现算法——Girvan-Newman算法
直观来看，在社区内部节点之间相互连接的边密度较大，因此，通过边介数（连接社区的边的边介数大）来识别社区是一种较为直观的社区发现算法。Girvan-Newman算法即是一种基于介数的社区发现算法，其基本思想是根据边介数中心性（edge betweenness）从大到小的顺序不断的将边从网络中移除直到整个网络分解为各个社区。因此，Girvan-Newman算法实际上是一种分裂方法。
该算法的基本流程如下：
基于模块度的社区发现算法——Louvain算法
Louvain算法是一种基于多层次优化Modularit（模块度）的算法，它的优点是快速、准确，被认为是性能最好的社区发现算法之一。Modularity函数最初被用于衡量社区发现算法结果的质量，它能够刻画发现的社区的紧密程度。那么既然能刻画社区的紧密程度，也就能够被用来当作一个优化函数，如果能够提升当前社区结构的modularity，即将结点加入它的某个邻居所在的社区中。
如果当前结点所在的社区只有它自己，那么在计算将它加入到其它社区时的modularity的变化有个技巧来加速计算，Louvain的高效性也在一定程度上受益于此，它具体如下:
$$
\begin{align*}
\Delta Q &=\left [ \frac{\sum_{in}^{}{+k_{i,in}} }{2m} -(\frac{\sum_{tot}^{}{+k_{i}} }{2m})^2 \right ]-\left [ \frac{\sum_{in}^{}}{2m} -(\frac{\sum_{tot}^{}}{2m})^2-(\frac{k_i}{2m})^2 \right ]\
&=\frac{1}{2m}(k_{i,in}-\frac{\sum_{tot}^{}k_i }{m} )
\end{align*}
$$
Louvain算法主要包括两个阶段：
遍历网络中的所有结点，尝试将单个结点加入能使modulatity提升最大的社区中，直到所有结点不再变化
将上一步的结果生成的一个个小社区归并为一个超结点来重新构造网络，这时边的权重为两个结点内所有原始结点的边权重之和。
算法基本流程如下：
四、算法实现
GN算法
将用于测试的网络数据存放于network.dat文件中，格式如下：
比如前两行就表示0，1这两个节点相连，利用networkx库进行无向图的绘制，输出如下：
经过GN算法，删除掉“桥梁”边之后，得到的图像如下：
可见总的节点被分为14个社团，任选其中一个社团放大，如下：
可以看出社团内部之间的连接很紧密，具体的社区发现结果我们保存在predict.txt文件中，该文件第一列表示节点，第二列表示社团编号，在这里截取一部分展示。
Louvain算法
同样是读取network.dat文件，然后利用Louvain算法进行社区发现，结果图如下：
放大右下角，如下：
得到的结果存放于Louvain_predict.txt文件中，截取一部分结果进行展示
五、效果评价部分
对于网络network.dat，它的实际社区分布保存在community.dat中，我们的社区发现算法结果都保存在txt文件中，利用上述文件计算两种算法的NMI指数：
GN算法
运行结果如下：
Louvain算法
运行结果如下：
简要分析
可以看出：两种算法的NMI指数都很高，说明两种社区发现算法的准确率都很高，与实际情况较为相符。其中GN算法的NMI指数略高于Louvain算法的NMI指数，可能的原因如下：
Louvain算法是基于模块度的社区发现算法，模块度也可以作为社区发现的评价指标，不过模块度主要是用于不知道真实社区分布情况下的评价指标（类似于无监督学习评价指标），而NMI指数主要用于已知真实情况下的评价。Louvain可能是过于侧重模块度的优化，注重的方面是社区更集中，而不是跟真实情况更相近，所以可能导致算法的准确率未达到90%，而GN算法考虑网络全局，准确率就较高。
六、算法复杂度分析
GN算法
GN算法计算边界数的时间复杂度为$\mathcal{O}(m\times n)$，其中m为边数，n为节点数。
总时间复杂度在m条边和n个节点的网络下为 $\mathcal{O}(m^2n)$。
Louvain算法
在Louvain算法的第一步中，扫描所有节点并将其加入能带来最高模块度收益的邻居节点所在的社区，每一次迭代的时间复杂度为$\mathcal{O}(n)$，其中的n为输入数据中的边的数量；
在第二步中完美将第一步形成的社区进行折叠，把每个社区看成一个超结点进行下一步的边权重计算，将结果用于下一次迭代的步骤一，每一次迭代中，步骤二的复杂度为$\mathcal{O}(m+n)$，m为本轮迭代中超结点的个数。
七、算法评价
GN算法
优点
- GN算法每次迭代删去最大边介数所在的边时，都在考虑网络全局，所以划分社区的准确率很高
缺点
- 因为每次迭代都要考虑整个网络，所以时间复杂度较高
- 对网络社区结构缺少量的定义，得事先知道社区的个数
Louvain算法
优点
- 简单、易于理解
- 非监督
- 时间复杂度低，运算速度很快，是最快的模块化算法之一
缺点
- 它会将较小的社区合并为较大的社区
- 当存在具有类似模块性的多个候选分区时，可能会出现平台，形成局部最大值并阻止进展
八、代码展示
GN算法
import networkx as nx
import matplotlib.pyplot as plt
def load_graph(path):
g = nx.Graph()
with open(path, 'r', encoding='utf8') as f:
lines = f.readlines()
for line in lines:
line = line.replace('\n', '')
nodex, nodey = line.split('\t')
# 因为我们是无权图，所以这里将所有的边权重设为1.0
g.add_edge(nodex, nodey, weight=1.0)
return g
# 移除边界数最大的边
def removing_based_on_betweeness(g):
init_ncomp = nx.number_connected_components(g) # 连通图的数量
curr_ncomp = init_ncomp
while curr_ncomp <= init_ncomp:
# 因为我们删除边界数大的边就是为了不连通
# 所以只要边界数大的边删除后，连通数量大于原来的数量就停止
bws = nx.edge_betweenness_centrality(g, weight='weight')
max_bw = max(bws.values())
for nodes, bw in bws.items():
if bw == max_bw:
g.remove_edge(*nodes) # 相当于删除当前边界数大的边
curr_ncomp = nx.number_connected_components(g) # 更新连通数量
def get_deg(g):
# 统计每个队和别人打的次数
adj = nx.adj_matrix(g)
nodes = g.nodes()
t = adj.sum(axis=1)
return {node: t[i, 0] for i, node in enumerate(nodes)}
def get_modularity(g, init_deg, m):
deg = get_deg(g)
modularity = 0
for comp in nx.connected_components(g):
for node in comp:
modularity += (deg[node] - init_deg[node] \
** 2 / (2 * m))
return modularity / (2 * m)
def gn(g):
init_n = g.number_of_edges()
print("总共的边数(相当于两两比赛的总次数):", init_n)
if init_n == 0:
return None
# 简单的一个预处理
m = nx.adj_matrix(g) # 按第一列从小到大 整理我们的数据
m = m.sum() / 2 # 总的比赛数
init_deg = get_deg(g)
i = 1
while g.number_of_edges() > 420:
removing_based_on_betweeness(g)
# 迭代5次计算一下模块度
if i % 5 == 0:
print("迭代次数:{} 模块度的值:{} 边数:{}".format\
(i, get_modularity(g, init_deg, m), \
g.number_of_edges()))
i += 1
print("模块度最大: {}".format(get_modularity(g, \
init_deg, m)))
# 社团中止时的图
nx.draw_networkx(g)
plt.savefig("gn_2.jpg",dpi = 600)
plt.show()
return nx.connected_components(g)
if __name__ == "__main__":
# 1. 加载数据
graph = load_graph("data/network.dat")
# 2. 画出我们的网络结构图
nx.draw_networkx(graph)
plt.savefig("gn_1.jpg",dpi = 600)
plt.show()
# 3. 进行社团检测
set_list = gn(graph)
total = []
for temp in set_list:
total.append(list(temp))
# 4.写入文件 这种写入文件的数据与标签文件的样子
# 相同 但是我们检测预测的准确率
result = []
for temp in total:
for i in temp:
result.append([i, total.index(temp)])
f = open('./data/predict.txt', 'w')
for r in result:
f.write(str(r[0]) + '\t')
f.write(str(r[1]) + '\n')
f.close()Louvain算法
import collections
import random
import networkx as nx
import matplotlib.pyplot as plt
def load_graph(path):
G = collections.defaultdict(dict)
with open(path) as text:
for line in text:
vertices = line.strip().split()
v_i = int(vertices[0])
v_j = int(vertices[1])
# w = float(vertices[2])
G[v_i][v_j] = 1
return G
class Vertex:
def __init__(self, vid, cid, nodes, k_in=0):
self._vid = vid
self._cid = cid
self._nodes = nodes
self._kin = k_in # 结点内部的边的权重
class Louvain:
def __init__(self, G):
self._G = G
self._m = 0 # 边数量
# 需维护的关于社区的信息(社区编号,其中包含的结点编号的集合)
self._cid_vertices = {}
# 需维护的关于结点的信息(结点编号，相应的Vertex实例)
self._vid_vertex = {}
for vid in self._G.keys():
self._cid_vertices[vid] = set([vid])
self._vid_vertex[vid] = Vertex(vid, vid, set([vid]))
self._m += sum([1 for neighbor in self._G[vid].keys() \
if neighbor > vid])
def first_stage(self):
mod_inc = False # 用于判断算法是否可终止
visit_sequence = self._G.keys()
random.shuffle(list(visit_sequence))
while True:
can_stop = True # 第一阶段是否可终止
for v_vid in visit_sequence:
v_cid = self._vid_vertex[v_vid]._cid
k_v = sum(self._G[v_vid].values()) + \
self._vid_vertex[v_vid]._kin
cid_Q = {}
for w_vid in self._G[v_vid].keys():
w_cid = self._vid_vertex[w_vid]._cid
if w_cid in cid_Q:
continue
else:
tot = sum(
[sum(self._G[k].values()) + \
self._vid_vertex[k]._kin for \
k in self._cid_vertices[w_cid]])
if w_cid == v_cid:
tot -= k_v
k_v_in = sum([v for k, v in self._G[v_vid]\
.items() if k in self._cid_vertices[w_cid]])
delta_Q = k_v_in - k_v * tot / self._m # 由于只需要知道delta_Q的正负，所以少乘了1/(2*self._m)
cid_Q[w_cid] = delta_Q
cid, max_delta_Q = sorted(cid_Q.items(), key=\
lambda item: item[1], reverse=True)[0]
if max_delta_Q > 0.0 and cid != v_cid:
self._vid_vertex[v_vid]._cid = cid
self._cid_vertices[cid].add(v_vid)
self._cid_vertices[v_cid].remove(v_vid)
can_stop = False
mod_inc = True
if can_stop:
break
return mod_inc
def second_stage(self):
cid_vertices = {}
vid_vertex = {}
for cid, vertices in self._cid_vertices.items():
if len(vertices) == 0:
continue
new_vertex = Vertex(cid, cid, set())
for vid in vertices:
new_vertex._nodes.update(self._vid_vertex[vid]._nodes)
new_vertex._kin += self._vid_vertex[vid]._kin
for k, v in self._G[vid].items():
if k in vertices:
new_vertex._kin += v / 2.0
cid_vertices[cid] = set([cid])
vid_vertex[cid] = new_vertex
G = collections.defaultdict(dict)
for cid1, vertices1 in self._cid_vertices.items():
if len(vertices1) == 0:
continue
for cid2, vertices2 in self._cid_vertices.items():
if cid2 <= cid1 or len(vertices2) == 0:
continue
edge_weight = 0.0
for vid in vertices1:
for k, v in self._G[vid].items():
if k in vertices2:
edge_weight += v
if edge_weight != 0:
G[cid1][cid2] = edge_weight
G[cid2][cid1] = edge_weight
self._cid_vertices = cid_vertices
self._vid_vertex = vid_vertex
self._G = G
def get_communities(self):
communities = []
for vertices in self._cid_vertices.values():
if len(vertices) != 0:
c = set()
for vid in vertices:
c.update(self._vid_vertex[vid]._nodes)
communities.append(c)
return communities
def execute(self):
iter_time = 1
while True:
iter_time += 1
mod_inc = self.first_stage()
if mod_inc:
self.second_stage()
else:
break
return self.get_communities()
def load(path):
g = nx.Graph()
with open(path, 'r', encoding='utf8') as f:
lines = f.readlines()
for line in lines:
line = line.replace('\n', '')
nodex, nodey = line.split('\t')
# 因为我们是无权图，所以这里将所有的边权重设为1.0
g.add_edge(nodex, nodey, weight=1.0)
return g
if __name__ == '__main__':
G = load_graph('./data/network.dat')
algorithm = Louvain(G)
communities = algorithm.execute()
# 按照社区大小从大到小排序输出
communities = sorted(communities, key=lambda b: -len(b)) # 按社区大小排序
print(G)
total = []
for temp in communities:
total.append(list(temp))
# 4.写入文件 这种写入文件的数据与标签文件的样子相同 但是我们检测预测的准确率
result = []
for temp in total:
for i in temp:
result.append([i, total.index(temp)])
f = open('./data/Louvain_predict.txt', 'w')
for r in result:
f.write(str(r[0]) + '\t')
f.write(str(r[1]) + '\n')
f.close()
graph = load("./data/Louvain_predict.txt")
# 5. 画出我们的网络结构图
nx.draw_networkx(graph)
plt.show()计算NMI指数
import numpy as np
import math
def NMI(A, B):
# 社团检测评估算法 NMI 基于归一化互信息向量的匹配方法
# 1. len(A) 应该等于 len(B)
total = len(A)
A_ids = set(A)
B_ids = set(B)
# 2. 互信息
MI = 0
eps = 1.4e-45
for idA in A_ids:
for idB in B_ids:
idAOccur = np.where(A == idA) # 找出A 在A_ids中的编号
idBOccur = np.where(B == idB) # 找出B 在B_ids中的编号
idABOccur = np.intersect1d(idAOccur, idBOccur) # 返回两个数组中相同值的个数
px = 1.0 * len(idAOccur[0]) / total # 每种类别所占比例
py = 1.0 * len(idBOccur[0]) / total # 每种类别所占比例
pxy = 1.0 * len(idABOccur) / total # 两种类别直接的联合概率
MI = MI + pxy * math.log(pxy / (px*py) + eps, 2)
# 3.归一化互信息
Hx = 0
for idA in A_ids:
idAOccurCount = 1.0 * len(np.where(A == idA)[0])
Hx = Hx - (idAOccurCount / total) * math.log(idAOccurCount / total + eps, 2)
Hy = 0
for idB in B_ids:
idBOccurCount = 1.0 * len(np.where(B == idB)[0])
Hy = Hy - (idBOccurCount / total) * math.log(idBOccurCount / total + eps, 2)
MIhat = 2.0 * MI / (Hx + Hy)
return MIhat
def load_data(path):
data = {}
with open(path, 'r', encoding='utf8') as f:
lines = f.readlines()
for line in lines:
g, t = line.replace('\n', '').split('\t')
data[g] = t
return data
if __name__ == '__main__':
path_real = './data/community.dat'
path_predict = './data/Louvain_predict.txt'
data_real = load_data(path_real)
data_predict = load_data(path_predict)
# print("真实数据:", data_real)
# print("预测数据:", data_predict)
# print(len(data_real))
# print(len(data_predict))
real, predict = [], []
for i in range(0, 115):
i = str(i)
real.append(data_real.get(i))
predict.append(data_predict.get(i))
real = np.array(real)
predict = np.array(predict)
print("Louvain算法准确率:\n", NMI(real, predict))
