前言
社区切分或者说社区监测，是一类图聚类算法，它主要作用是将图数据划分为不同的社区，社区内的节点都是连接紧密或者相似的，而社区与社区之间的节点连接则是稀疏的或者是社区与社区之间的节点并不相似。社区切分算法目前已经广泛运用在各个领域，例如在金融反欺诈的场景中，社区切分算法可以根据资金交易图谱划分出可疑的欺诈团伙或者洗钱团伙。再就是在营销场景中，社区切分算法可以对不同的人群进行自然分类，帮助业务人员发现不同特性的人群，然后进行精准营销活动的投放。还有在智能运维场景中去发现高风险网络架构等等。社区切分算法是很重要的图谱数据分析算法，是所有进行图数据挖掘领域从业人员必须要掌握的技术。本篇文章将重点介绍一种基础的社区切分算法Louvain，它其实也是一种半监督的图聚类算法，通过模块度增益来作为约束条件衡量划分效果，而且也有着不错的计算性能，是目前广泛运用的图聚类算法之一。
参考资料
[1] M. E. J. Newman (2003). “Mixing patterns in networks”. Physical Review E. 67 (2): 026126. https://arxiv.org/abs/cond-mat/0209450
[2]M. E. J. Newman (2002). “Assortative mixing in networks”. Physical Review Letters. 89 (20): 208701.https://arxiv.org/abs/cond-mat/0205405
[3] Newman, M. E. J. (2006). “Modularity and community structure in networks”. Proceedings of the National Academy of Sciences of the United States of America. 103 (23): 8577–8696. https://arxiv.org/abs/physics/0602124
[4]M.E.J, Newman; M.Girvan (2004). “Finding and evaluating community structure in networks”. Phys. Rev. E. 69 (2): 026113.https://arxiv.org/abs/cond-mat/0308217
[5] Blondel, Vincent D; Guillaume, Jean-Loup; Lambiotte, Renaud; Lefebvre, Etienne (9 October 2008). “Fast unfolding of communities in large networks”. Journal of Statistical Mechanics: Theory and Experiment. 2008 (10): P10008. https://arxiv.org/abs/0803.0476
[6] Modularity (networks). https://en.wikipedia.org/wiki/Modularity_(networks)
[7] 大龙1993. “网络算法系列之社区发现（二）：模块度”. https://juejin.cn/post/6844904120202035213
[8] 薄德瑜、王啸. “异配图神经网络小结：突破同配性，拥抱异配性”. 北邮GAMMA Lab. https://www.sohu.com/a/500262605_121119001
同配系数（assortative coefficient）
在讲解模块度定义之前，会有两个补充知识点，因为模块度的定义也是基于这两个知识点发展来的，其中一个就是同配系数。我们先来聊一聊什么是图的同配混合性质（assortative mixing），它指的是经过研究发现，现实生活中很多的关系其实会呈现一种同配性质（assortative），也就是具有相同特征的节点之间会倾向于相互连接。那比如说出国的留学生们会发现，除了平时上课和做作业时，相同国籍且有相同兴趣爱好的人其实会更喜欢形成小圈子，所以会经常看见韩国学生周末聚在一起喝酒，中国学生周末聚在一起打麻将打游戏，印度学生周末聚在一起打板球…………
那么，相对的另外一些现实世界上的网络则是呈现非同配性质（disassortative），例如相亲网站中的用户多数带有比较强的目的性，配对成功的情侣首先性别肯定不一样，其次就是学历和收入也会出现较大差异，所以我们可以发现在这类社交网络中，节点与节点之间的特征会出现较大不同，这些异质的节点则倾向于相互连接。
为了去衡量网络中同配性质的程度，Newman在《Mixing patterns in networks》论文中提出了针对节点中离散特征的同配系数公式，用来衡量网络的同配性质的程度：
r = ∑ i e i i − ∑ i a i b i 1 − ∑ i a i b i = T r ( e ) − ∥ e 2 ∥ 1 − ∥ e 2 ∥ （ 1 ） r= \frac{\sum_{i}e_{ii}-\sum_{i}a_ib_i}{1-\sum_{i}a_ib_i} =\frac{Tr(e) -\begin{Vmatrix}e^2\end{Vmatrix}}{1-\begin{Vmatrix}e^2\end{Vmatrix}} \qquad （1） r=1−∑iaibi∑ieii−∑iaibi=1− e2 Tr(e)− e2 （1）
其中， e i j e_{ij} eij表示连接类型为 i i i的节点和类型为 j j j的节点之间的边数量与总的边数量的比值， ∑ i j e i j = 1 \sum_{ij}e_{ij}=1 ∑ijeij=1。
a i = ∑ j e i j a_{i}=\sum_{j}e_{ij} ai=∑jeij 指的是连接到类型为 i i i的节点的所有边数量与总的边数量的比值。这里类型为 i i i的节点是Source节点。
b j = ∑ i e i j b_{j}=\sum_{i}e_{ij} bj=∑ieij 指的是连接到类型为 j j j的节点的所有边数量与总的边数量的比值。这里类型为 j j j的节点是Target节点。
e e e 指的是一个 i × j i \times j i×j的矩阵，矩阵中的元素是 e i j e_{ij} eij。
T r ( e ) Tr(e) Tr(e)表示的是矩阵的迹，也就是主对角线上元素的相加。 ∥ e 2 ∥ \begin{Vmatrix}e^2\end{Vmatrix} e2 指的是矩阵中的所有元素先平方然后再相加。
当r = 0时， 表示整个网络都没有同配混合性，此时 e i j = a i b j e_{ij} =a_ib_j eij=aibj，或者说 ∑ i e i i = ∑ i a i b i \sum_{i}e_{ii} =\sum_{i}a_ib_i ∑ieii=∑iaibi。
当r = 1 时，表示此时整个网络是具有完全的同配混合性， ∑ i e i i = 1 \sum_{i}e_{ii}=1 ∑ieii=1，表示所有的点只和同类型的点相连。
当 r m i n = − ∑ i a i b i 1 − ∑ i a i b i r_{min}= \frac{-\sum_{i}a_ib_i}{1-\sum_{i}a_ib_i} rmin=1−∑iaibi−∑iaibi时，也就是 − 1 ≤ r < 0 -1 \leq r<0 −1≤r<0, 则认为整个网络是具有完全的不同配混合性，表示所有的点只和非同类型的点相连，或者说是每一条边连接的两个点都是类型不同的。
对于带有连续型特征变量的网络，它的同配系数公式则可以写成如下形式：
r = ∑ x y x y ( e x y − a x b y ) σ a σ b （ 2 ） r= \frac{\sum_{xy}xy(e_{xy}-a_xb_y)}{\sigma_a\sigma_b} \qquad （2） r=σaσb∑xyxy(exy−axby)（2）
其中， e x y e_{xy} exy表示的是拥有连续型变量 x x x的所有节点和拥有连续型变量 y y y的所有节点之间相连的边数量占总的边数量的比值。
类比离散型变量网络的同配系数，我们同样有以下规则， ∑ x y e x y = 1 \sum_{xy}e_{xy}=1 ∑xyexy=1， ∑ y e x y = a x \sum_{y}e_{xy}=a_x ∑yexy=ax， ∑ x e x y = b y \sum_{x}e_{xy}=b_y ∑xexy=by。 a x a_x ax就代表着拥有连续型变量值 x x x的Source节点连接其它变量值的节点的边总数与总的边数量的比值。同理， b y b_y by。
σ a σ b \sigma_a\sigma_b σaσb分别代表 a x a_x ax和 b y b_y by的标准差。所以我们可以看出用来衡量拥有连续型变量网络的同配系数的公式本质上也是节点之间度数的皮尔逊相关系数( Pearson correlation coefficient），这里就不做过多展开了，感兴趣的小伙伴可以再去读一下原论文。
另外，关于网络的同配性质需要注意的是，正如《异配图神经网络小结：突破同配性，拥抱异配性》这篇文章提到的，GNN模型在具有同配性的网络中效果较好，但是在具有异配性质的网络中，GNN模型的效果还不如利用图结构信息的多层感知机 。所以，我们在使用GNN模型时，最好也要分析网络的同配性系数。
节点之间边数量的数学期望值
我们上面已经讲解完同配系数了，另外一个需要补充的知识点就是节点之间边连接数量的数学期望值。假设有两个节点，一个是 v v v,另一个是 w w w, 他们的节点度数分别是 k v k_v kv和 k w k_w kw。另外，我们还需要再引入一个概念，叫节点的存根（stub），它来自《Random Graphs and Complex Networks》这篇论文中提出的configuration model，其实就是将网络中所有的边都切成两个半边，每个半边就是一个stub。然后每个stub都会和另外一个stub随机连接，以此来生成随机图。后续我们所提到的随机图都是基于configuration model来讨论。
所以，节点 v v v拥有 k v k_v kv个stub。每个stub我们都用一个 I i ( v , w ) I_i^{(v,w)} Ii(v,w)来表示当前 v v v节点的第 i i i个stub是否有和 w w w节点的某个stub相连，如果相连，就是 I i ( v , w ) = 1 I_i^{(v,w)}=1 Ii(v,w)=1，否则 I i ( v , w ) = 0 I_i^{(v,w)}=0 Ii(v,w)=0。其中， i = 1 , 2 , . . . . . . , k v i=1,2,......,k_v i=1,2,......,kv。因为节点 v v v某个stub与其余所有stub连接的概率都为 1 2 m − 1 \frac{1}{2m-1} 2m−11（2m表示整个网络中stub的总数）, 又因为节点 v v v与节点 w w w可以相连的stub数量为 k w k_w kw，所以我们可以得出下面的两个公式：
p ( I i ( v , w ) = 1 ) = E [ I i ( v , w ) ] = k w 2 m − 1 （ 3 ） p(I_i^{(v,w)}=1)=E[I_i^{(v,w)}]=\frac{k_w}{2m-1}\qquad （3） p(Ii(v,w)=1)=E[Ii(v,w)]=2m−1kw（3）
J v w = ∑ i = 1 k v I i ( v , w ) （ 4 ） J_{vw}=\sum_{i=1}^{k_v}I_i^{(v,w)}\qquad\qquad\qquad\qquad （4） Jvw=∑i=1kvIi(v,w)（4）
J v w J_{vw} Jvw代表着节点 v v v和节点 w w w之间所有边的总数量，通过公式(3)和公式(4)我们可以得出 J v w J_{vw} Jvw的数学期望为：
E [ J v w ] = E [ ∑ i = 1 k v I i ( v , w ) ] = ∑ i = 1 k v E [ I i ( v , w ) ] = ∑ i = 1 k v k w 2 m − 1 = k v k w 2 m − 1 ≈ k v k w 2 m （ 5 ） E[J_{vw}]=E[\sum_{i=1}^{k_v}I_i^{(v,w)}] = \sum_{i=1}^{k_v}E[I_i^{(v,w)}]=\sum_{i=1}^{k_v}\frac{k_w}{2m-1}=\frac{k_vk_w}{2m-1}\approx\frac{k_vk_w}{2m}\qquad （5） E[Jvw]=E[i=1∑kvIi(v,w)]=i=1∑kvE[Ii(v,w)]=i=1∑kv2m−1kw=2m−1kvkw≈2mkvkw（5）
在大型的随机网络中，可以假设任意两个节点之间最多有一条边， k v k w 2 m \frac{k_vk_w}{2m} 2mkvkw也可以近似为两个节点 v , w v,w v,w之间存在边的概率。
Newman的第一版模块度定义
正如之前所说的，Louvain是一种半监督的聚类算法，它通过模块度的增益来作为社区划分的约束条件。模块度的定义最早是由Newman和Girvan在《Finding and evaluating community structure in networks》这篇论文中提出来的，它可以用来评估社区切分质量的”好坏“。模块度的定义是从我们之前说的同配系数的基础上发展来的，它的公式定义如下：
Q = ∑ i ( e i i − a i 2 ) = T r ( e ) − ∥ e 2 ∥ （ 6 ） Q= \sum_{i}(e_{ii}-a_i^2)=Tr(e)-\begin{Vmatrix}e^2\end{Vmatrix} \qquad（6） Q=i∑(eii−ai2)=Tr(e)− e2 （6）
对于划分为k个社区的网络， e e e代表的是一个 k × k k \times k k×k的矩阵，里面的元素 e i j e_{ij} eij代表所有连接社区 i i i节点到社区 j j j节点的边的数量与总的边数量的比值。
矩阵 e e e的迹 T r ( e ) = ∑ i e i i Tr(e) = \sum_{i}e_{ii} Tr(e)=∑ieii, 且 e i i e_{ii} eii代表着两头都连接属于同社区 i i i内的点的边数量与总边数的比值。所以，很容易联想到一个好的社区划分，矩阵 e e e的迹 T r ( e ) Tr(e) Tr(e)就越高。但是当 T r ( e ) = 1 Tr(e)=1 Tr(e)=1时，也有可能是将所有节点都分到了同一个社区里，这也是不合理的。所以上面的公式中又引入了 a i = ∑ j e i j a_i = \sum_{j}e_{ij} ai=∑jeij , 代表所有连接到社区 i i i中的节点的边数量与总边数的比值。
原文中还有一句定义：“In a network in which edges fall between vertices without regard for the communities they belong to, we would have e i j = a i a j e_{ij} = a_ia_j eij=aiaj”，这句话我目前不太理解，翻阅了很多资料也没有看见进一步的解释，有些晦涩难懂，这里还想请教一下看懂的读者，如果您理解其中的含义欢迎在评论区积极留言一起讨论。根据类比之前Newman对同配系数的定义，我的理解是当节点不属于任何社区时， e i j = a i a j e_{ij} = a_ia_j eij=aiaj也就是代表着 e i j 、 a i 、 a j e_{ij}、a_i、a_j eij、ai、aj都为1，因为也相当于所有节点都在一个社区里。
模块度公式（6）其实也代表着所有连接相同社区节点的边在全图中的比例，减去一个拥有相同社区划分、相同节点度数的随机网络中社区之间边连接数量的数学期望值。当相同社区内的边数量与全图边数量的比值不如随机网络生成的边数量的数学期望，那么 Q = 0 Q=0 Q=0。如果 Q = 1 Q=1 Q=1，则代表着这个网络有着很好的社区划分结构。一般从实践来看，很多网络的模块度处于 0.3 < Q < 0.7 0.3 < Q < 0.7 0.3<Q<0.7区间。
另外，需要特别强调的是，当一个网络被划分为多个社区时，连接不同社区的边其实会被计数两次。所以， e i i = 2 × 社区 i 内部边的总数量 2 × 边的总数量 m e_{ii}=\frac{2\times 社区i内部边的总数量}{2 \times 边的总数量m} eii=2×边的总数量m2×社区i内部边的总数量， e i j = 连接社区 i 到社区 j 的边的总数量 2 × 边的总数量 m e_{ij}=\frac{连接社区i到社区j的边的总数量}{2 \times 边的总数量m} eij=2×边的总数量m连接社区i到社区j的边的总数量。至于社区 i i i到社区 j j j的边的总数量为什么不乘以2，显而易见是因为 e j i e_{ji} eji还会再计数一次。这也是模块度被算错的主要原因，大家要特别注意。
下面我们通过下面一个10个节点和12条边的网络来演示如何计算模块度。
上面的图显示出已经被两条黄色的线划分成了三个社区，绿色的节点属于社区A，红色的节点属于社区B，蓝色的节点属于社区C。那么三个社区的矩阵 e e e如下所示：
那么 Q = 2 ∗ 3 2 ∗ 12 − ( 2 ∗ 3 + 1 2 ∗ 12 ) 2 + 2 ∗ 4 2 ∗ 12 − ( 1 + 2 ∗ 4 + 1 2 ∗ 12 ) 2 + 2 ∗ 3 2 ∗ 12 − ( 1 + 2 ∗ 3 2 ∗ 12 ) 2 = 0.4896 Q=\frac{2*3}{2*12}-(\frac{2*3+1}{2*12})^2+\frac{2*4}{2*12}-(\frac{1+2*4+1}{2*12})^2+\frac{2*3}{2*12}-(\frac{1+2*3}{2*12})^2=0.4896 Q=2∗122∗3−(2∗122∗3+1)2+2∗122∗4−(2∗121+2∗4+1)2+2∗122∗3−(2∗121+2∗3)2=0.4896
Newman的第二版模块度定义
为了使得社区检测更加高效，Newman老爷子在《Modularity and community structure in networks》这篇论文中提出了第二版模块度的定义，引入了模块度矩阵的概念，也是为了使用谱图算法去实现社区检测。
不过，我们在这里不对上述论文中的社区检测算法展开详述，感兴趣的朋友可以自己去读一下论文。我们只对这个新的模块度定义进行讲解，因为后面要讲解的Louvain算法主要也是基于这一版的模块度定义去实现的。
那么，我们先来看公式,下面是当节点划分为两个社区时的模块度定义：
Q = 1 2 m ∑ v w [ A v w − k v k w 2 m ] s v s w + 1 2 ( 7 ) Q=\frac{1}{2m}\sum_{vw}[A_{vw}-\frac{k_vk_w}{2m} ]\frac{s_vs_w+1}{2} \qquad(7) Q=2m1vw∑[Avw−2mkvkw]2svsw+1(7)
其中， A v w A_{vw} Avw表示的是节点 v v v和节点 w w w之间实际的边数量；根据我们上面公式（5）的定义，两节点之间存在的边数量的数学期望值可以表示为 k v k w 2 m \frac{k_vk_w}{2m} 2mkvkw；其中 k v k_v kv和 k w k_w kw分别是节点 v v v和 w w w的度数。 s v s_v sv有-1和1两种取值，表示节点 v v v属于社区1时， s v = 1 s_v=1 sv=1，节点属于社区2时， s v = − 1 s_v=-1 sv=−1。 2 m 2m 2m表示的是全图所有stub的数量，也就是 2 m = ∑ v k v 2m=\sum_{v}k_v 2m=∑vkv。
当我们将社区划分为2个以上时，假设有 i i i个社区，那么公式（7）可以写成：
Q = 1 2 m ∑ v w [ A v w − k v k w 2 m ] δ ( c v , c w ) ( 8 ) Q=\frac{1}{2m}\sum_{vw}[A_{vw}-\frac{k_vk_w}{2m} ]\delta(c_v,c_w) \qquad(8) Q=2m1vw∑[Avw−2mkvkw]δ(cv,cw)(8)
δ ( c v , c w ) \delta(c_v,c_w) δ(cv,cw)是一个函数，当节点 v v v和节点 w w w在同一社区时，即 c v = c w c_v=c_w cv=cw，那么 δ ( c v , c w ) = 1 \delta(c_v,c_w)=1 δ(cv,cw)=1，否则属于不同社区，那么就为0。
我们再将公式（8）做一下变换，将有如下形式：
Q = 1 2 m [ ∑ v w A v w − ∑ v k v ∑ w k w 2 m ] δ ( c v , c w ) Q=\frac{1}{2m}[\sum_{vw}A_{vw}-\frac{\sum_{v}k_v\sum_{w}k_w}{2m} ]\delta(c_v,c_w) Q=2m1[∑vwAvw−2m∑vkv∑wkw]δ(cv,cw)
= ∑ v w A v w δ ( c v , c w ) 2 m − ∑ v k v ∑ w k w δ ( c v , c w ) ( 2 m ) 2 \quad=\frac{\sum_{vw}A_{vw}\delta(c_v,c_w)}{2m}-\frac{\sum_{v}k_v\sum_{w}k_w\delta(c_v,c_w)}{(2m)^2} =2m∑vwAvwδ(cv,cw)−(2m)2∑vkv∑wkwδ(cv,cw)
= ∑ i [ ∑ i n 2 m − ( ∑ t o t 2 m ) 2 ] \quad=\sum_{i}[\frac{\sum_{in}}{2m}-(\frac{\sum_{tot}}{2m})^2] =∑i[2m∑in−(2m∑tot)2]
= ∑ i ( e i i − a i 2 ) ( 9 ) \quad=\sum_{i}(e_{ii}-a_i^2) \qquad(9) =∑i(eii−ai2)(9)
上述公式中， ∑ i n \sum_{in} ∑in表示的是当前社区 i i i中内部边的总数量 × 2 \times2 ×2 （社区内部的边指的是边的两头连接的节点都属于同一个社区。 × 2 \times2 ×2 是因为每条边连接着两个stub，或者说一条边会被计数两次。）， ∑ t o t \sum_{tot} ∑tot 则表示的是当前社区 i i i中所有节点的总度数（或者说是社区 i i i中内部边,以及连接其它社区的外部边的总数量 × 2 \times2 ×2）
经过公式（9）的推导，我们发现第二版本的模块度定义公式在多社区划分的场景下和第一版本的模块度定义公式达到了一致。
Louvain社区切分算法
大家学习了同配系数、节点之间边数量的数学期望值、模块度等知识点后，接下来学习Louvain算法原理将会非常轻松，也将了解到Louvain如何通过模块度增益来作为约束条件对图进行层次性的社区划分。
Louvain算法来自《Fast unfolding of communities in large networks》这篇论文，作者并没有给这个算法命名，后来大家就以第一作者所在的Louvain大学来命名该社区切分算法。Louvain算法的优势在于可以用于海量数据的大图上，且计算的效率很高，1.18亿节点的图在单台主机上用了152分钟（其实挺慢的，所以我后面讲解的源码是基于Spark的版本）。Louvain算法主要由两个步骤组成。
Step1: 首先，对于一个有N个节点的图，初始化时每个节点都属于一个不包含其它节点的独立社区，也就是有N个社区。接下来，每个节点 i i i都会去考虑是否加入邻接节点 j j j的社区，如果加入某个邻接节点 j j j的社区获得最大的模块度增益，且是正数，那么节点 i i i就会加入节点 j j j的社区；如果无论加入哪一个邻接节点 j j j的社区都是负的模块度增益，则节点 i i i这一轮就不加入任何社区。
计算一个节点 i i i加入社区 C C C的模块度增益可以由以下公式方便地计算出来：
△ Q \triangle Q △Q =节点 i i i加入社区C后的模块度 - （节点 i i i加入前的社区C的模块度+节点 i i i离开原本社区后的模块度变化值）
节点i加入社区C后的模块度= [ ∑ i n + k i , i n + ∑ i n k i 2 m − ( ∑ t o t + k i 2 m ) 2 ] [\frac{\sum_{in}+k_{i,in}+\sum_{in}^{k_i}}{2m}-(\frac{\sum_{tot}+k_i}{2m})^2] [2m∑in+ki,in+∑inki−(2m∑tot+ki)2]。其中， ∑ i n \sum_{in} ∑in表示的是社区C的所有内部边的权重和； k i , i n k_{i,in} ki,in表示的是节点 i i i连接到社区C内节点的边的权重和； ∑ t o t \sum_{tot} ∑tot指的是连接到社区C中所有节点的边的权重和,这里既包括了社区C的内部边，也包括了社区C连接到其它社区的外部边； k i k_i ki指的是连接到节点 i i i的所有边的权重和； ∑ i n k i 2 m \frac{\sum_{in}^{k_i}}{2m} 2m∑inki表示的是节点 i i i所在社区内部边的权重和； m m m代表着全图边的权重和。
大家应该注意到这里有两个问题，一是不再是统计边的总数量，而是统计边的权重和；二是初始化时节点 i i i所在的社区只有 i i i节点，为什么会有内部边？这其实都和Louvain算法的第二个步骤相关，后续会介绍到。
节点 i i i加入前的社区C的模块度= [ ∑ i n 2 m − ( ∑ t o t 2 m ) 2 ] [\frac{\sum_{in}}{2m}-(\frac{\sum_{tot}}{2m})^2] [2m∑in−(2m∑tot)2]
节点 i i i离开原本社区后的模块度变化值= [ ∑ i n k i 2 m − ( k i 2 m ) 2 ] [\frac{\sum_{in}^{k_i}}{2m}-(\frac{k_i}{2m})^2] [2m∑inki−(2mki)2]
根据上面的定义，计算一个节点 i i i加入社区 C C C的模块度增益就可以写成：
△ 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 ] ( 10 ) \triangle Q =[\frac{\sum_{in}+k_{i,in}}{2m}-(\frac{\sum_{tot}+k_i}{2m})^2]-[\frac{\sum_{in}}{2m}-(\frac{\sum_{tot}}{2m})^2-(\frac{k_i}{2m})^2]\qquad(10) △Q=[2m∑in+ki,in−(2m∑tot+ki)2]−[2m∑in−(2m∑tot)2−(2mki)2](10)
公式（10）中对 ∑ i n k i 2 m \frac{\sum_{in}^{k_i}}{2m} 2m∑inki进行了消解，同时公式（10）也是论文中对模块度增益计算公式的定义。我们继续对公式（10）进行展开可以得到以下化简公式：
△ Q = 1 2 m ( k i , i n − ∑ t o t k i m ) ( 11 ) \triangle Q =\frac{1}{2m}(k_{i,in}-\frac{\sum_{tot}k_i}{m})\qquad(11) △Q=2m1(ki,in−m∑totki)(11)
整个过程会一直循环直到所有节点所属社区不再变化，也就是模块度不再有进一步的提升，那么当前step1阶段就结束了。
Step2: 根据第一阶段为每个节点分配的社区生成一个全新的网络，网络中的节点是上一个阶段所生成的社区，节点之间的边权重是对应两个社区之间所连接的边的总数量。可以简单地理解为将step1阶段划分的社区内的所有节点，在step2都合并为一个新的节点，新节点之间的连接其实也是上一个阶段社区之间的连接。生成新的网络后再继续重复Step1的步骤，一直到节点所在的社区不再变化，也代表着模块度达到最大值，这时候就停止迭代。
至此，社区切分的算法原理已经全部讲解完了，下面我们将进入到源码解读阶段。
Louvain算法源码解读
接下来，我们将从一个开源框架 distributed-graph-analytics去学习如何用代码实现Louvain社区切分。因为我们平时在处理工业界的图数据时，都是上亿的节点和边，所以这也是为什么去选择一个分布式的图计算框架来讲解。
因为该框架是基于Spark Graphx去实现的图计算方法，所以Louvain算法的开发语言是Scala。这里建议大家具备一些Spark和Scala的基础知识，不然下面的源码解读比较难以理解。
另外，下面的讲解中，提到的代码行数对应的是源码文件中实际的代码行数。
Louvain算法Step1源码解读
Louvain算法的第一步就是先初始化，每一个节点都分配一个只包含自己的社区，也就是N个节点的图初始化后有N个社区。我们先来看核心类 LouvainCore.scala，里面就实现了初始化的步骤。
代码块1:
/**
* Generates a new graph of type Graph[VertexState,Long] based on an input graph of type.
* Graph[VD,Long]. The resulting graph can be used for louvain computation.
*
*/
def createLouvainGraph[VD: ClassTag](graph: Graph[VD,Long]) : Graph[VertexState,Long]= {
// Create the initial Louvain graph.
val nodeWeightMapFunc = (e:EdgeTriplet[VD,Long]) => Iterator((e.srcId,e.attr), (e.dstId,e.attr))
val nodeWeightReduceFunc = (e1:Long,e2:Long) => e1+e2
val nodeWeights = graph.mapReduceTriplets(nodeWeightMapFunc,nodeWeightReduceFunc)
val louvainGraph = graph.outerJoinVertices(nodeWeights)((vid,data,weightOption)=> {
val weight = weightOption.getOrElse(0L)
val state = new VertexState()
state.community = vid
state.changed = false
state.communitySigmaTot = weight
state.internalWeight = 0L
state.nodeWeight = weight
state
}).partitionBy(PartitionStrategy.EdgePartition2D).groupEdges(_+_)
return louvainGraph
}
上面的createLouvainGraph方法就是初始化图谱并用于后续的Louvain计算。
第30-32行，通过graphx的mapreduce方法计算图上所有节点连接的边的权重和，也就是计算 k i k_i ki。
第34-42行则是通过outerJoinVertices方法将节点的连接边权重和结果合并回节点自身，并创建了一个VertexState对象来保存节点的信息。这里面的community属性记录的是节点所属的社区id，初始化阶段社区id就是节点本身的id，且只包含节点自己。communitySigmaTot属性记录了节点所在社区的所有内部边和外部边的权重和，也就是记录 ∑ t o t \sum_{tot} ∑tot，这里初始化时用的是前面计算出来的节点的weight。internalWeight记录的是当前节点所属社区的内部边的权重和，也就是记录 ∑ i n \sum_{in} ∑in或者 ∑ i n k i \sum_{in}^{k_i} ∑inki的值。nodeWeight就是记录之前计算的节点的连接边权重和，也是节点自己本身的权重。
第43行则是合并节点之间的多条边，并将边上的权重进行相加，赋予合并后的边。这里可以减少需要计算的图的大小。
接下来，我们看到 LouvainCore.scala中的louvain方法。
代码块2:
val graphWeight = louvainGraph.vertices.values.map(vdata=> vdata.internalWeight+vdata.nodeWeight).reduce(_+_)
var totalGraphWeight = sc.broadcast(graphWeight)
上面的代码块来自louvain方法中的第67行-68行，其作用是计算全图的模块度，并作为Spark的广播变量。
代码块3：
var msgRDD = louvainGraph.mapReduceTriplets(sendMsg,mergeMsg)
代码块4:
/**
* Creates the messages passed between each vertex to convey neighborhood community data.
*/
private def sendMsg(et:EdgeTriplet[VertexState,Long]) = {
val m1 = (et.dstId,Map((et.srcAttr.community,et.srcAttr.communitySigmaTot)->et.attr))
val m2 = (et.srcId,Map((et.dstAttr.community,et.dstAttr.communitySigmaTot)->et.attr))
Iterator(m1, m2)
}
代码块5:
/**
* Merge neighborhood community data into a single message for each vertex
*/
private def mergeMsg(m1:Map[(Long,Long),Long],m2:Map[(Long,Long),Long]) ={
val newMap = scala.collection.mutable.HashMap[(Long,Long),Long]()
m1.foreach({case (k,v)=>
if (newMap.contains(k)) newMap(k) = newMap(k) + v
else newMap(k) = v
})
m2.foreach({case (k,v)=>
if (newMap.contains(k)) newMap(k) = newMap(k) + v
else newMap(k) = v
})
newMap.toMap
}
代码块3来自louvain方法中的第72行，利用Graphx的消息传播机制，并通过sendMsg和mergeMsg两个自定义方法将节点周围的邻居节点信息传递过来，这里包括周围节点所属的社区ID、节点所属社区的所有内外部边的权重和（ ∑ t o t \sum_{tot} ∑tot）、以及周围节点连接的边的权重和。"et.attr"是当前边的权重，mergeMsg将所有相同社区的节点的边权重进行聚合相加，并返回一个map对象保存每个社区的边权重相加后的结果。
代码块6:
// label each vertex with its best community based on neighboring community information
val labeledVerts = louvainVertJoin(louvainGraph,msgRDD,totalGraphWeight,even).cache()
代码块7:
/**
* Join vertices with community data form their neighborhood and select the best community for each vertex to maximize change in modularity.
* Returns a new set of vertices with the updated vertex state.
*/
private def louvainVertJoin(louvainGraph:Graph[VertexState,Long], msgRDD:VertexRDD[Map[(Long,Long),Long]], totalEdgeWeight:Broadcast[Long], even:Boolean) = {
louvainGraph.vertices.innerJoin(msgRDD)( (vid, vdata, msgs)=> {
var bestCommunity = vdata.community
var startingCommunityId = bestCommunity
var maxDeltaQ = BigDecimal(0.0);
var bestSigmaTot = 0L
msgs.foreach({ case( (communityId,sigmaTotal),communityEdgeWeight ) =>
val deltaQ = q(startingCommunityId, communityId, sigmaTotal, communityEdgeWeight, vdata.nodeWeight, vdata.internalWeight,totalEdgeWeight.value)
//println(" communtiy: "+communityId+" sigma:"+sigmaTotal+" edgeweight:"+communityEdgeWeight+" q:"+deltaQ)
if (deltaQ > maxDeltaQ || (deltaQ > 0 && (deltaQ == maxDeltaQ && communityId > bestCommunity))){
maxDeltaQ = deltaQ
bestCommunity = communityId
bestSigmaTot = sigmaTotal
}
})
// only allow changes from low to high communties on even cyces and high to low on odd cycles
if ( vdata.community != bestCommunity && ( (even && vdata.community > bestCommunity) || (!even && vdata.community < bestCommunity) ) ){
//println(" "+vid+" SWITCHED from "+vdata.community+" to "+bestCommunity)
vdata.community = bestCommunity
vdata.communitySigmaTot = bestSigmaTot
vdata.changed = true
}
else{
vdata.changed = false
}
vdata
})
}
当完成邻居节点的信息聚合后，就会进入louvain方法的循环体中，执行模块度计算，并进行必要的节点合并。上面的代码块6来自 LouvainCore.scala的第86行，主要是根据周围节点的信息计算模块度增益，然后更新节点所属的社区信息，这里通过调用louvainVertJoin这个方法实现。
代码块7是louvainVertJoin的方法实现，主要是通过.innerJoin(msgRDD)将当前节点的周围节点社区情况关联起来，然后调用q这个方法进行模块度的增益计算。 msgs.foreach这个算子会遍历计算当前节点加入到周围某一个节点社区中时的模块度增益，如果增益度是正数，且和加入其它节点相比，它的增益为最大，那么就加入该社区，同时将当前节点的社区信息则更新为所加入的社区的信息。
代码块8:
/**
* Returns the change in modularity that would result from a vertex moving to a specified community.
*/
private def q(currCommunityId:Long, testCommunityId:Long, testSigmaTot:Long, edgeWeightInCommunity:Long, nodeWeight:Long, internalWeight:Long, totalEdgeWeight:Long) : BigDecimal = {
val isCurrentCommunity = (currCommunityId.equals(testCommunityId));
val M = BigDecimal(totalEdgeWeight);
val k_i_in_L = if (isCurrentCommunity) edgeWeightInCommunity + internalWeight else edgeWeightInCommunity;
val k_i_in = BigDecimal(k_i_in_L);
val k_i = BigDecimal(nodeWeight + internalWeight);
val sigma_tot = if (isCurrentCommunity) BigDecimal(testSigmaTot) - k_i else BigDecimal(testSigmaTot);
var deltaQ = BigDecimal(0.0);
if (!(isCurrentCommunity && sigma_tot.equals(0.0))) {
deltaQ = k_i_in - ( k_i * sigma_tot / M)
//println(s" $deltaQ = $k_i_in - ( $k_i * $sigma_tot / $M")
}
return deltaQ;
}
代码块8就是计算节点加入某一个具体的社区时模块度的增益。q方法的输入参数分别为：当前节点的社区id、要加入的社区id、要加入社区的内外边的权重和（ ∑ t o t \sum_{tot} ∑tot）、节点连接到要加入社区的边的权重和（ k i , i n k_{i,in} ki,in）、当前节点外部边的权重和、当前节点内部边的权重和（ ∑ i n k i \sum_{in}^{k_i} ∑inki）。我们可以看到 k i k_i ki是由节点 i i i的外部边权重和以及内部边的权重和相加得出的。“deltaQ = k_i_in - ( k_i * sigma_tot / M)”这行代码其实就是对应着我们上面提到的模块度增益计算公式（11）。
代码块9:
do {
count += 1
even = ! even
// label each vertex with its best community based on neighboring community information
val labeledVerts = louvainVertJoin(louvainGraph,msgRDD,totalGraphWeight,even).cache()
// calculate new sigma total value for each community (total weight of each community)
val communtiyUpdate = labeledVerts
.map( {case (vid,vdata) => (vdata.community,vdata.nodeWeight+vdata.internalWeight)})
.reduceByKey(_+_).cache()
// map each vertex ID to its updated community information
val communityMapping = labeledVerts
.map( {case (vid,vdata) => (vdata.community,vid)})
.join(communtiyUpdate)
.map({case (community,(vid,sigmaTot)) => (vid,(community,sigmaTot)) })
.cache()
// join the community labeled vertices with the updated community info
val updatedVerts = labeledVerts.join(communityMapping).map({ case (vid,(vdata,communityTuple) ) =>
vdata.community = communityTuple._1
vdata.communitySigmaTot = communityTuple._2
(vid,vdata)
}).cache()
updatedVerts.count()
labeledVerts.unpersist(blocking = false)
communtiyUpdate.unpersist(blocking=false)
communityMapping.unpersist(blocking=false)
val prevG = louvainGraph
louvainGraph = louvainGraph.outerJoinVertices(updatedVerts)((vid, old, newOpt) => newOpt.getOrElse(old))
louvainGraph.cache()
// gather community information from each vertex's local neighborhood
val oldMsgs = msgRDD
msgRDD = louvainGraph.mapReduceTriplets(sendMsg, mergeMsg).cache()
activeMessages = msgRDD.count() // materializes the graph by forcing computation
oldMsgs.unpersist(blocking=false)
updatedVerts.unpersist(blocking=false)
prevG.unpersistVertices(blocking=false)
// half of the communites can swtich on even cycles
// and the other half on odd cycles (to prevent deadlocks)
// so we only want to look for progess on odd cycles (after all vertcies have had a chance to move)
if (even) updated = 0
updated = updated + louvainGraph.vertices.filter(_._2.changed).count
if (!even) {
println(" # vertices moved: "+java.text.NumberFormat.getInstance().format(updated))
if (updated >= updatedLastPhase - minProgress) stop += 1
updatedLastPhase = updated
}
} while ( stop <= progressCounter && (even || (updated > 0 && count < maxIter)))
接下来，我们再回到Louvain方法的循环体中，也就是代码块9。之前的讲解中，我们已经详细阐述了louvainVertJoin(louvainGraph,msgRDD,totalGraphWeight,even)方法如何完成模块度增益的计算，并为每一个节点都重新分配了社区。接下来在代码块9中，也就是源码文件的第88行-113行，则基于节点新分配的社区，重新计算该社区的内外边的权重和（ ∑ t o t \sum_{tot} ∑tot）。和之前一样，也是属于该社区所有节点的内外部边的权重和（nodeWeight + internalWeight) 作为社区总的权重和（ ∑ t o t \sum_{tot} ∑tot）。最后，将重新计算后的社区信息都更新进入节点的VertexState对象中。
在代码块9中，我们可以发现Louvain方法的循环体默认是在两次循环后结束，至少是所有的节点都进行了社区的变化（留在当前社区或者加入其它节点的社区）。一次循环是奇循环，只允许高社区加入低社区，另外一次是偶循环，只允许低社区加入高社区，社区的高低由社区id的大小来判断，这样设计是为了保证节点加入社区时有一个人为设计的优先顺序，避免产生节点之间来回加入对方社区的死循环。具体的代码实现体现在源码文件的第124行-第133行，以及louvainVertJoin方法中的第212行。
上面的源码讲解对应的就行Louvain算法的Step1，计算模块度增益，然后为节点分配新的社区。
Louvain算法Step2源码解读
接下来我们讲解Step2的源码，这个阶段的部分实现在AbstractLouvainRunner.scala类中，另一部分实现来自LouvainCore.scala中的compressGraph方法。
代码块10:
do {
compressionLevel += 1
println(s"\nStarting Louvain level $compressionLevel")
// label each vertex with its best community choice at this level of compression
val (currentQModularityValue, currentGraph, numberOfPasses) = louvainCore.louvain(sc, louvainGraph, minimumCompressionProgress, progressCounter)
louvainGraph.unpersistVertices(blocking = false)
louvainGraph = currentGraph
// If modularity was increased by at least 0.001 compress the graph and repeat
// halt immediately if the community labeling took less than 3 passes
//println(s"if ($passes > 2 && $currentQ > $q + 0.001 )")
//if (numberOfPasses > 2 && currentQModularityValue > q_modularityValue + 0.001) {
if (currentQModularityValue > q_modularityValue + 0.001) {
saveLevel(sc, compressionLevel, currentQModularityValue, louvainGraph)
q_modularityValue = currentQModularityValue
louvainGraph = louvainCore.compressGraph(louvainGraph)
}
else {
halt = true
}
} while (!halt)
finalSave(sc, compressionLevel, q_modularityValue, louvainGraph)
}
上面的代码块10是一个循环体，它的作用是让图谱一直根据模块度增益为节点分配新的社区。在这个过程中，相同社区的节点合成新的节点和新的关系，并生成新的压缩图，接着用压缩图进行新一轮的模块度增益计算，当前总的模块度相比上一轮迭代的模块度增加幅度小于0.001时，则社区切分的过程就会停止。
代码块11:
/**
* Compress a graph by its communities, aggregate both internal node weights and edge
* weights within communities.
*/
def compressGraph(graph:Graph[VertexState,Long],debug:Boolean=true) : Graph[VertexState,Long] = {
// aggregate the edge weights of self loops. edges with both src and dst in the same community.
// WARNING can not use graph.mapReduceTriplets because we are mapping to new vertexIds
val internalEdgeWeights = graph.triplets.flatMap(et=>{
if (et.srcAttr.community == et.dstAttr.community){
Iterator( ( et.srcAttr.community, 2*et.attr) ) // count the weight from both nodes // count the weight from both nodes
}
else Iterator.empty
}).reduceByKey(_+_)
// aggregate the internal weights of all nodes in each community
var internalWeights = graph.vertices.values.map(vdata=> (vdata.community,vdata.internalWeight)).reduceByKey(_+_)
// join internal weights and self edges to find new interal weight of each community
val newVerts = internalWeights.leftOuterJoin(internalEdgeWeights).map({case (vid,(weight1,weight2Option)) =>
val weight2 = weight2Option.getOrElse(0L)
val state = new VertexState()
state.community = vid
state.changed = false
state.communitySigmaTot = 0L
state.internalWeight = weight1+weight2
state.nodeWeight = 0L
(vid,state)
}).cache()
// translate each vertex edge to a community edge
val edges = graph.triplets.flatMap(et=> {
val src = math.min(et.srcAttr.community,et.dstAttr.community)
val dst = math.max(et.srcAttr.community,et.dstAttr.community)
if (src != dst) Iterator(new Edge(src, dst, et.attr))
else Iterator.empty
}).cache()
// generate a new graph where each community of the previous
// graph is now represented as a single vertex
val compressedGraph = Graph(newVerts,edges)
.partitionBy(PartitionStrategy.EdgePartition2D).groupEdges(_+_)
// calculate the weighted degree of each node
val nodeWeightMapFunc = (e:EdgeTriplet[VertexState,Long]) => Iterator((e.srcId,e.attr), (e.dstId,e.attr))
val nodeWeightReduceFunc = (e1:Long,e2:Long) => e1+e2
val nodeWeights = compressedGraph.mapReduceTriplets(nodeWeightMapFunc,nodeWeightReduceFunc)
// fill in the weighted degree of each node
// val louvainGraph = compressedGraph.joinVertices(nodeWeights)((vid,data,weight)=> {
val louvainGraph = compressedGraph.outerJoinVertices(nodeWeights)((vid,data,weightOption)=> {
val weight = weightOption.getOrElse(0L)
data.communitySigmaTot = weight +data.internalWeight
data.nodeWeight = weight
data
}).cache()
louvainGraph.vertices.count()
louvainGraph.triplets.count() // materialize the graph
newVerts.unpersist(blocking=false)
edges.unpersist(blocking=false)
return louvainGraph
}
代码块11是整个Louvain算法step2部分的核心。
当step1为每个节点分配好新社区后，compressGraph方法会将相同社区内的节点和边合并成新的节点和新的边，并生成新的压缩图。
首先compressGraph方法会计算新社区内部边的权重和（ ∑ i n \sum_{in} ∑in），在压缩图中会将新社区内所有节点进行合并生成新的节点。新节点的internalWeight属性就是用来存储该新社区内部边的权重和（ ∑ i n \sum_{in} ∑in）。在源码中，我们可以看到内部边的权重和由weight1和weight2组成。weight1就是连接了相同社区节点的边的权重 × 2 \times2 ×2，对应源码中的“2*et.attr”这部分。这里为什么乘以2大家可以看回我们之前模块度第一版公式的定义中。weight2就是将相同社区中所有节点的internalWeight的属性值相加，其实这里就是上一轮迭代中加入该社区的其它社区的内部边权重和（ ∑ i n \sum_{in} ∑in）。
当更新完internalWeight后，就开始生成合并后的新节点、合并后新的边，以及计算合并后新社区总的权重和（ ∑ t o t \sum_{tot} ∑tot）。最终完成新的压缩图生成后，新图会加入到下一轮模块度计算的迭代中。
至此，整个Louvain算法的代码实现部分就讲解完了，本质就是利用GraphX框架中的消息传递机制，将每个节点周围的社区信息进行聚合，并以此来计算模块度的增益来决定节点是否需要加入新的社区。大家只要理解代码中对应的 k i , i n k_{i,in} ki,in、 ∑ t o t \sum_{tot} ∑tot、 k i k_i ki等核心参数的计算方法，就不难理解其背后的实现逻辑。