Louvain社区发现算法出自2008年的论文《Fast unfolding of communities in large networks》,其名字是根据作者所在的城市来命名的。它基于模块度优化来实现社区划分。
准备知识
模块度(modularity)是用来衡量社区内部的链接密度相比社区之间的链接密度的介于-1和1的分数。对于带权网络,其定义如下:
Q = 1 2 m ∑ i , j [ A i j − k i k j 2 m ] δ ( c i , c j ) ( 1 ) Q = \frac{1}{2m} \sum_{i,j} \left[A_{ij} - \frac{k_i k_j}{2m} \right] \delta(c_i, c_j) \qquad (1) Q=2m1i,j∑[Aij−2mkikj]δ(ci,cj)(1)
上式中的 A i j A_{ij} Aij是节点i和节点j之间的权重, k i = ∑ j A i j k_i=\sum_j A_{ij} ki=∑jAij是与节点i相连的边的权重之和, c i c_i ci是节点i所属的社区,如果u=v,则 δ ( u , v ) \delta(u,v) δ(u,v)为1否则为0。 m = 1 2 ∑ i j A i j m = \frac{1}{2} \sum_{ij} A_{ij} m=21∑ijAij(图中所有边的权重之和,因为一条边的权重被计算了两次,所以要除以2;当权重为1时,m即图中边的条数)。
根据上述定义,社区c的模块度可被计算为:
Q c = 1 2 m ∑ i ∑ j A i j 1 { c i = c j = c } − ( ∑ i k i 2 m 1 { c i = c } ) = ∑ i n 2 m − ( ∑ t o t 2 m ) 2 ( 2 ) \begin{aligned} Q_c &= \frac{1}{2m} \sum_i \sum_j A_{ij} \mathbb{1}\{c_i = c_j = c \} - \left( \sum_i \frac{k_i}{2m} \mathbb{1} \{c_i =c \} \right) \\ &= \frac{\sum_{in}}{2m} - \left(\frac{\sum_{tot}}{2m} \right)^2 \qquad (2) \end{aligned} Qc=2m1i∑j∑Aij1{ci=cj=c}−(i∑2mki1{ci=c})=2m∑in−(2m∑tot)2(2)
上式中 ∑ i n \sum_{in} ∑in是社区c内部的边的权重之和(每条边会被计算两次), ∑ t o t \sum_{tot} ∑tot是所有与c内节点相连的边的权重之和(包括与其他社区相连的边)。网络的模块度也记为:
Q = ∑ c Q c ( 3 ) Q = \sum_c Q_c \qquad (3) Q=c∑Qc(3)
Louvain算法详情
Louvain算法分为两阶段,第一阶段是模块度优化,第二阶段是社区合并。不断重复这个阶段直到网络的模块度不再发生变化或者达到最大迭代次数。
第一阶段:
在初始化时,将每个节点都当做一个独立的社区,接着对每一个节点i,计算将i从原来的社区移到其邻居节点j对应的社区后的模块度增益,将i的社区变成增益最大的邻居对应的社区,注意这里的增益只考虑正增益,如果最大增益为负则i的社区保持不变。这个过程对所有节点不断地重复直到无法因为改变某个节点的社区得到正模块度增益。
Lovain算法比较高效是因为将孤立节点i加入到社区C的模块度增益计算可按下式进行(注:下式中的 k i , i n k_{i,in} ki,in 网上资料有两个版本,有的前面有系数2,有的没有,arxiv上和networkx都是使用下面式子的版本,如果从前面定义的社区c的模块度来看,应该是下面的式子,式(4)的第一部分是加入节点i后的社群C的模块度,第二部分的前两项是社区c的模块度,第二项是节点i的模块度),
Δ Q = [ ∑ i n + k i , i n 2 m − ( ∑ t o t + k i 2 m ) 2 ] − [ ∑ i n 2 m − ( ∑ t o t 2 m ) 2 − ( k i 2 m ) 2 ] ( 4 ) = k i , i n 2 m − ∑ t o t ⋅ k i 2 m 2 ( 5 ) \begin{aligned} \Delta Q &= \left[ \frac{\sum_{in} + k_{i, in}}{2m} - \left( \frac{\sum_{tot} + k_i}{2m} \right)^2 \right] - \left[ \frac{\sum_{in} }{2m} - \left(\frac{\sum_{tot} }{2m}\right)^2 - \left(\frac{k_i }{2m}\right)^2 \right] \qquad (4)\\ &= \frac{k_{i,in}}{2m} - \frac{\sum_{tot} \cdot k_i}{2m^2} \qquad (5) \end{aligned} ΔQ=[2m∑in+ki,in−(2m∑tot+ki)2]−[2m∑in−(2m∑tot)2−(2mki)2](4)=2mki,in−2m2∑tot⋅ki(5)
∑ i n \sum_{in} ∑in是社区C内部的边的权重之和, ∑ t o t \sum_{tot} ∑tot是所有与C内节点相连的边的权重之和(包括与其他社区相连的边), k i k_i ki是与节点i相连的边的权重之和, k i , i n k_{i, in} ki,in是节点i与社区C内的节点之间的边的权重之和,m是图中所有边的权重之和。
每个节点的最终社区取决于节点遍历顺序,作者说实验结果显示遍历顺序对最终模块度没有显著影响,但是会影响计算时间。因此实际实现时都会对节点进行随机排序后再进行遍历。
第二阶段:
第二阶段将前一阶段得到每个社区合并成一个新的节点构成新的图。新节点之间的边的权重为对应的两个社区里节点的边的权重之和,而原来社区内部节点之间的边当做新节点自环(self-loop)其权重为原来边的权重之和。
Louvain算法的优点:
- 时间复杂度低,计算快适合大规模的网络。
- 社区划分结果稳定,可用模块度来评估社区划分好坏。
- 算法流程使其划分结果有层次化,如果将每一次迭代的社区划分结果保存下来,可使用其中间结果。
- 因为其层次化特性可以规避分辨率问题,因为在初始化时每个节点都是一个独立的社区。
Louvain算法的缺点:
- 计算过程存储在内存中,如果图很大,算法对内存要求高。
- 不容易进行分布式实现(面有针对分布式的改进Louvain算法)
Louvain算法实现
Networkx中实现了Louvain算法,使用很简单:
import networkx as nx
G = nx.karate_club_graph()
nx.community.louvain_communities(G, seed=123)
csdn blog实现的代码理解算法更容易
import collections
import random
import time
import networkx as nx
import matplotlib.pyplot as pltdef 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 = 1.0 # 数据集有权重的话则读取数据集中的权重G[v_i][v_j] = wG[v_j][v_i] = wreturn G# 节点类 存储社区与节点编号信息
class Vertex:def __init__(self, vid, cid, nodes, k_in=0):# 节点编号self._vid = vid# 社区编号self._cid = cidself._nodes = nodesself._kin = k_in # 结点内部的边的权重class Louvain:def __init__(self, G, seed=123):self._G = Gself.seed = seedself._m = 0 # 边数量 图会凝聚动态变化self._cid_vertices = {} # 需维护的关于社区的信息(社区编号,其中包含的结点编号的集合)self._vid_vertex = {} # 需维护的关于结点的信息(结点编号,相应的Vertex实例)for vid in self._G.keys():# 刚开始每个点作为一个社区self._cid_vertices[vid] = {vid}# 刚开始社区编号就是节点编号self._vid_vertex[vid] = Vertex(vid, vid, {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.Random(self.seed).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节点的权重(度数) 内部与外部边权重之和k_v = sum(self._G[v_vid].values()) + \self._vid_vertex[v_vid]._kin# 存储模块度增益大于0的社区编号cid_Q = {}# 遍历节点的邻居for w_vid in self._G[v_vid].keys():# 获得该邻居的社区编号w_cid = self._vid_vertex[w_vid]._cidif w_cid in cid_Q:continueelse:# tot是关联到社区C中的节点的链路上的权重的总和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是从节点i连接到C中的节点的链路的总和k_v_in = sum([v for k, v in self._G[v_vid].items() if k in self._cid_vertices[w_cid]])# 由于只需要知道delta_Q的正负,所以少乘了1/(2*self._m)delta_Q = k_v_in - k_v * tot / self._mcid_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 = Falsemod_inc = Trueif can_stop:breakreturn mod_inc# 网络凝聚阶段def second_stage(self):cid_vertices = {}vid_vertex = {}# 遍历社区和社区内的节点for cid, vertices in self._cid_vertices.items():if len(vertices) == 0:continuenew_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# k,v为邻居和它们之间边的权重 计算kin社区内部总权重 这里遍历vid的每一个在社区内的邻居 因为边被两点共享后面还会计算 所以权重/2for k, v in self._G[vid].items():if k in vertices:new_vertex._kin += v / 2.0# 新的社区与节点编号cid_vertices[cid] = {cid}vid_vertex[cid] = new_vertexG = collections.defaultdict(dict)# 遍历现在不为空的社区编号 求社区之间边的权重for cid1, vertices1 in self._cid_vertices.items():if len(vertices1) == 0:continuefor cid2, vertices2 in self._cid_vertices.items():# 找到cid后另一个不为空的社区if cid2 <= cid1 or len(vertices2) == 0:continueedge_weight = 0.0# 遍历 cid1社区中的点for vid in vertices1:# 遍历该点在社区2的邻居已经之间边的权重(即两个社区之间边的总权重 将多条边看做一条边)for k, v in self._G[vid].items():if k in vertices2:edge_weight += vif edge_weight != 0:G[cid1][cid2] = edge_weightG[cid2][cid1] = edge_weight# 更新社区和点 每个社区看做一个点self._cid_vertices = cid_verticesself._vid_vertex = vid_vertexself._G = Gdef 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(list(c))return communitiesdef execute(self):iter_time = 1while True:iter_time += 1# 反复迭代,直到网络中任何节点的移动都不能再改善总的 modularity 值为止mod_inc = self.first_stage()if mod_inc:self.second_stage()else:breakreturn self.get_communities()# 可视化划分结果
def showCommunity(G, partition, pos):# 划分在同一个社区的用一个符号表示,不同社区之间的边用黑色粗体cluster = {}labels = {}for index, item in enumerate(partition):for nodeID in item:labels[nodeID] = r'$' + str(nodeID) + '$' # 设置可视化labelcluster[nodeID] = index # 节点分区号# 可视化节点colors = ['r', 'g', 'b', 'y', 'm']shapes = ['v', 'D', 'o', '^', '<']for index, item in enumerate(partition):nx.draw_networkx_nodes(G, pos, nodelist=item,node_color=colors[index],node_shape=shapes[index],node_size=350,alpha=1)# 可视化边edges = {len(partition): []}for link in G.edges():# cluster间的linkif cluster[link[0]] != cluster[link[1]]:edges[len(partition)].append(link)else:# cluster内的linkif cluster[link[0]] not in edges:edges[cluster[link[0]]] = [link]else:edges[cluster[link[0]]].append(link)for index, edgelist in enumerate(edges.values()):# cluster内if index < len(partition):nx.draw_networkx_edges(G, pos,edgelist=edgelist,width=1, alpha=0.8, edge_color=colors[index])else:# cluster间nx.draw_networkx_edges(G, pos,edgelist=edgelist,width=3, alpha=0.8, edge_color=colors[index])# 可视化labelnx.draw_networkx_labels(G, pos, labels, font_size=12)plt.axis('off')plt.show()def cal_Q(partition, G): # 计算Q# 如果为真,则返回3元组(u、v、ddict)中的边缘属性dict。如果为false,则返回2元组(u,v)m = len(G.edges(None, False))# print(G.edges(None,False))# print("=======6666666")a = []e = []for community in partition: # 把每一个联通子图拿出来t = 0.0for node in community: # 找出联通子图的每一个顶点# G.neighbors(node)找node节点的邻接节点t += len([x for x in G.neighbors(node)])a.append(t / (2 * m))# self.zidian[t/(2*m)]=communityfor community in partition:t = 0.0for i in range(len(community)):for j in range(len(community)):if (G.has_edge(community[i], community[j])):t += 1.0e.append(t / (2 * m))q = 0.0for ei, ai in zip(e, a):q += (ei - ai ** 2)return qclass Graph:graph = nx.DiGraph()def __init__(self):self.graph = nx.DiGraph()def createGraph(self, filename):file = open(filename, 'r')for line in file.readlines():nodes = line.split()edge = (int(nodes[0]), int(nodes[1]))self.graph.add_edge(*edge)return self.graphif __name__ == '__main__':G = load_graph('data/club.txt')G1 = nx.karate_club_graph()pos = nx.spring_layout(G1)start_time = time.time()algorithm = Louvain(G)communities = algorithm.execute()end_time = time.time()# 按照社区大小从大到小排序输出communities = sorted(communities, key=lambda b: -len(b)) # 按社区大小排序count = 0for communitie in communities:count += 1print("社区", count, " ", communitie)print(cal_Q(communities, G1))print(f'算法执行时间{end_time - start_time}')# 可视化结果showCommunity(G1, communities, pos)
其他
另外,Louvain作者在2023年11月又写了一篇论文《Fast unfolding of communities in large networks: 15 years later》,介绍Louvain算法的实现与相关改进方法。
参考资料
-
Fast unfolding of communities in large networks
-
csdn blog, blog对应的python实现
-
networkx里实现的Louvain