社区发现（子区划分）：模块度及模块度增量计算公式推导

社区发现：模块度$Q$及模块度增量$\Delta Q$计算公式

重点参考该文章（基本上是全网最详细的模块度计算说明，非常值得一看）：

Modularity的计算方法——社团检测中模块度计算公式详解 | 雅乐网 (yalewoo.com)

一、模块度计算公式

1.1、相关概念

1.1.1、图

图是一种由节点、边组成的复杂网状结构

• 按照边的有无方向可以分为有向图、无向图；按照边有无权重可以分为有权图、无权图；

• 相邻节点：两个节点通过一条边相连（图中节点1、2相邻）；

• 节点的度：与节点相连的边的数量（图中节点1度为2）

节点邻接矩阵$A$：

• $A_{ij}$表示节点$i、j$是否相连，不相连则$A_{ij}=0$，否则$A_{ij}=1$（如果是加权图，则为权重值）

$$\begin{array}{c}{A}_{ij}=\{\begin{array}{rl}1& i,j相邻\\ 0& i,j不相邻\end{array}\end{array}$$

节点	1	2	3	4	5
1	0	1	0	0	1
2	1	0	1	1	0
3	0	1	0	1	0
4	0	1	1	0	1
5	1	0	0	1	0

1.1.2、社区发现

社区：社区是图的普遍属性，一个图中可能有多个社区，同一社区内节点与节点之间关系紧密，社区与社区之间的关系稀疏

社区发现（Community Detection）：可以当作图中结构紧密的节点的聚类，即将各个节点划分为不同的社区，使得社区内节点之间连接紧密，社区间连接稀疏

• 与聚类算法同属非监督学习，核心都是节点间距离度量、节点所属类别判断、“聚类”效果评价；
• 经典社区发现算法以模块度作为社区紧密程度的度量指标，以模块度最大为“聚类”的目标，以模块度增量最大确定节点所属的社区
• 工程应用包括交通领域的控制子区划分：将交叉口作为节点，路段作为边，将包含几十个交叉口的大路网划分为若干个“小区”，分别进行信号控制。

模块度$Q$：是社区内节点紧密程度的度量指标，是衡量社区划分效果最常用的评价指标，值越大表明社区划分效果越好

社区关联矩阵$E$：

• $e_{ii}$表示子区$i$内边数与图中边数的比值；

• $e_{ij},i\ne j$表示子区$i$与子区$j$相连的边数与图中边数的比值的$1/2$（因为同一条边会被$e_{ij}、e_{ji}$重复计算）；

• $a_i={\sum }_j{e}_{ij}$，表示子区i内的节点度数和与图中所有节点度数和的比值；

• 图中共6条边（度数和为12），划分为3个子区，其中子区$I$内有1条边，子区$I、II$之间有1条边，子区$I、III$之间有2条边则$e_{11}=1/6=2/12，e_{22}=1/12，e_{22}=2/12，a_1=1/6+2/12+1/12=5/12$

子区	I	II	III
I	1/6	1/12	2/12
II	1/12	0	1/12
III	2/12	1/12	1/6

1.2、公式及证明

说明：以无权无向图为例进行证明

1.2.1、公式形式一

$$\begin{gathered} Q=\sum _i(e_{ii}-a_i^2)=tr(E)-||E^2|| \\ 其中，e_{ii}为社区关联矩阵E中对角线上的元素，表示子区i内边数与图中边数的比值 \\ {a}_i=\sum _j{e}_{ij}，表示子区i内的节点度数和与图中所有节点度数和的比值 \\ ||E^2||表示矩阵E^2的各个元素之和 \end{gathered}$$

证明：${\sum }_i(e_{ii}-a_i^2)=tr(E)-||E^2||$
$$\begin{gathered} \sum _i(e_{ii}-a_i^2)=\sum _i{e}_{ii}-\sum _i{a}_i^2 \\ 其中，\sum _i{e}_{ii}=tr(E)，仅需证明\sum _i{a}_i^2=||E^2|| \\ 对于下图所示E矩阵，计算E^2时： \\ 对于E_{11},与第一行的元素相乘形成3个多项式，分别为：e_{11}{e}_{11},e_{11}{e}_{12},e_{11}{e}_{13} \\ 对于E_{21},与第二行的元素相乘形成3个多项式，分别为：e_{21}{e}_{11},e_{21}{e}_{12},e_{21}{e}_{13} \\ 对于E_{31},与第三行的元素相乘形成3个多项式，分别为：e_{31}{e}_{11},e_{31}{e}_{12},e_{31}{e}_{13} \\ 则e_{11}{e}_{11}+e_{11}{e}_{12}+e_{11}{e}_{13}+e_{21}{e}_{11}+e_{21}{e}_{12}+e_{21}{e}_{13}+e_{31}{e}_{11}+e_{31}{e}_{12}+e_{31}{e}_{13} \\ =e_{11}(e_{11}+e_{12}+e_{13})+e_{21}(e_{11}+e_{12}+e_{13})+e_{31}(e_{11}+e_{12}+e_{13}) \\ =(e_{11}+e_{21}+e_{31})(e_{11}+e_{12}+e_{13}) \\ 因为E为对称矩阵，所以(e_{11}+e_{21}+e_{31})(e_{11}+e_{12}+e_{13})=(e_{11}+e_{12}+e_{13}{)}^2 \\ 即第一行的元素之和的平方，由此可知||E^2||为各行元素和的平方累加 \\ 即，||E^2||=\sum _i(\sum _j{e}_{ij}{)}^2=\sum _i{a}_i^2 \end{gathered}$$

1.2.2、公式形式二

$$\begin{gathered} Q=\frac{1}{2m}\sum _{ij}[A_{ij}-\frac{k_i{k}_j}{2m}]\delta (i,j) \\ 其中，m为图中边数，A_{ij}为邻接矩阵A中的元素，k_i、k_j分别为节点i、j的度 \\ \begin{array}{c}\delta (i,j)=\{\begin{array}{rl}1& i,j属于同一个社区\\ 0& i,j不属于同一个社区\end{array}\end{array} \end{gathered}$$

对于该公式形式，可以将模块度理解为各个社区内连边数与随机期望的差值之和，当实际的边数越高于随机的期望时，这个图的节点就越集中。

对于图中的节点1、2，实际连边数$A_{12}=1$，期望连边数=$k_1\frac{k_2}{2m}=2\frac{3}{12}=\frac{1}{2}$

即，期望连边数=$\frac{k_1{k}_2}{2m}$，证明详见：Modularity的计算方法——社团检测中模块度计算公式详解 | 雅乐网 (yalewoo.com)

1.2.3、公式等价性

证明：$Q={\sum }_i(e_{ii}-a_i^2)=tr(E)-||E^2||=\frac{1}{2m}{\sum }_{ij}[A_{ij}-\frac{k_i{k}_j}{2m}]\delta (i,j)$

其中：

对于左式：${\sum }_i(e_{ii}-a_i^2)={\sum }_i{e}_{ii}-{\sum }_i{a}_i^2$

对于右式：$\frac{1}{2m}{\sum }_{ij}[A_{ij}-\frac{k_i{k}_j}{2m}]\delta (i,j)=\frac{1}{2m}{\sum }_r{\sum }_{mn}{A}_{mn}^r-\frac{1}{(2m{)}^2}{\sum }_r{\sum }_{mn}{k}_m^r{k}_n^r$

${\sum }_{mn}$表示计算属于同一子区的节点之间的指标和，${\sum }_r$表示计算各个子区的指标和

分别证明${\sum }_i{e}_{ii}=\frac{1}{2m}{\sum }_r{\sum }_{mn}{A}_{mn}^r以及{\sum }_i{a}_i^2=\frac{1}{(2m{)}^2}{\sum }_r{\sum }_{mn}{k}_m^r{k}_n^r$

• 对于${\sum }_i{e}_{ii}=\frac{1}{2m}{\sum }_r{\sum }_{mn}{A}_{mn}^r$

$$\begin{gathered} \frac{1}{2m}\sum _r\sum _{mn}{A}_{mn}^r=\sum _r\sum _{mn}\frac{A_{mn}^r}{2m} \\ 对于子区r内的节点m、n,A_{mn}^r=A_{nm}^r \\ 因此\sum _r\sum _{mn}\frac{A_{mn}^r}{2m}=\sum _r\frac{子区内边数}{m}=\sum _i{e}_{ii} \end{gathered}$$

两个表达式都表示各个子区内的边数与总边数的比值之和，因此相等

• 对于${\sum }_i{a}_i^2=\frac{1}{(2m{)}^2}{\sum }_r{\sum }_{mn}{k}_m^r{k}_n^r$

$$\begin{gathered} a_i=\sum _j{e}_{ij}，表示子区i内的节点度数和与图中所有节点度数和的比值 \\ 即，a_i=\frac{度数{和}_i}{2m} \\ 因此，\sum _i{a}_i^2=\frac{{\sum }_i(度数{和}_i{)}^2}{(2m{)}^2} \\ 对于某一个子区，(度数{和}_i{)}^2=(k_1+k_2+...{)}^2=\sum _{mn}{k}_m^r{k}_n^r \\ 即子区内任意两点度数乘积之和 \\ 因此，\sum _i{a}_i^2=\frac{{\sum }_i(度数{和}_i{)}^2}{(2m{)}^2} \end{gathered}$$

1.3、模块度计算示例

1.3.1、对于公式形式一

$$\begin{gathered} e_{11}=1/6，e_{22}=0，e_{33}=1/6 \\ {a}_1=1/6+1/12+2/12=5/12，a_2=1/12+0+1/12=2/12，a_3=2/12+1/12+1/6=5/12 \\ Q=(1/6+0+1/6)-(5/12{)}^2-(2/12{)}^2-(5/12{)}^2=-1/24 \end{gathered}$$

1.3.2、对于公式形式二

令$B_{ij}=A_{ij}-\frac{k_i{k}_j}{2m}$，矩阵$B$如下：

节点	1	2	3	4	5
1	-1/3	1/2	-1/3	-1/2	2/3
2	1/2	-3/4	1/2	1/4	-1/2
3	-1/3	1/2	-1/3	1/2	-1/3
4	-1/2	1/4	1/2	-3/4	1/2
5	2/3	-1/2	-1/3	1/2	-1/3

根据$\delta (i,j)$的含义以及图中各个节点是否相邻确定$\delta (i,j)$的值，如下：

节点	1	2	3	4	5
1	1	1	0	0	0
2	1	1	0	0	0
3	0	0	1	0	0
4	0	0	0	1	1
5	0	0	0	1	1

计算得到模块度$Q=-1/24$

二、模块度计算公式的矩阵形式

矩阵形式公式参考该文章，其中原论文中也有类似的公式形式

• http://t.csdn.cn/9Q7Su
• http://www.pnas.org/content/103/23/8577.full.pdf

2.1、证明过程

$$\begin{gathered} Q=\frac{1}{2m}\sum _{ij}[A_{ij}-\frac{k_i{k}_j}{2m}]\delta (i,j)=\frac{1}{2m}tr(S^{T}BS) \\ 其中，B_{ij}=A_{ij}-\frac{k_i{k}_j}{2m} \end{gathered}$$

公式推导：

对于左式：

${\sum }_{ij}\delta (i,j)B_{ij}$表示对$B$中所属同一个子区的元素求和（包含对角线上的元素）

对于右式：

• $S$为$N\ast R$的矩阵，其中的元素$S_i^j$表示节点$i$是否属于子区$j$

$$\begin{gathered} \begin{array}{c} \\ {S}_i^j=\{\begin{array}{rl}1& i属于子区j\\ 0& i不属于子区j\end{array}\end{array} \end{gathered}$$

• 根据矩阵迹的交换律（章节2.2），可知$tr(S^{T}BS)=tr(S{S}^{T}B)$

• 根据矩阵$S$的特性可知$S{S}^T$为$N\ast N$的矩阵，并且其中的元素$S{S}_{ij}^T$表示节点$i,j$是否属于同一个子区

$$\begin{gathered} \begin{array}{c} \\ S{S}_{ij}^T=\{\begin{array}{rl}1& i,j属于同一个子区\\ 0& i,j不属于同一个子区\end{array}\end{array} \end{gathered}$$

• $S{S}^{T}B$为$N\ast N$的矩阵

$$\begin{gathered} tr(S{S}^{T}B)=\sum _{i=1}^N{S{S}^{T}B}_{ii}=\sum _{i=1}^N\sum _{k=1}^{N}S{S}_{ik}^T{B}_{ki} \\ 其中， \\ \begin{array}{c} \\ S{S}_{ik}^T{B}_{ki}=\{\begin{array}{rl}{B}_{ki}& i,k属于同一个子区\\ 0& i,k不属于同一个子区\end{array}\end{array} \end{gathered}$$

则$tr(S{S}^{T}B)$即为对$B$中所属同一个子区的元素求和（包含对角线上的元素）

因此，$tr(S{S}^{T}B)={\sum }_{ij}\delta (i,j)B_{ij}$

$Q=\frac{1}{2m}tr(S^{T}BS)$，得证

2.2、矩阵迹

（1）定义

在线性代数中，一个$N\ast N$的矩阵$A$的迹（或迹数），是指$A$的主对角线（从左上方至右下方的对角线）上各个元素的总和，一般记作$tr(A)$或$sp(A)$，$tr(A)=A_{11}+A_{22}+...+A_{NN}={\sum }_{i=1}^N{A}_{ii}$

（2）示例
$$A=\left[\begin{array}{ccc}3& 5& 1\\ 0& 9& 2\\ 7& 6& 4\end{array}\right]$$
则，$tr(A)=A_{11}+A_{22}+A_{33}=3+9+4=16$

（3）性质
$$\begin{gathered} tr(A+B)=tr(A)+tr(B) \\ tr(k\ast A)=k\ast tr(A),k为常数 \\ tr(A)=tr(A^T) \\ tr(AB)=tr(BA) \end{gathered}$$
其中，前3个性质显然成立，同时求解模块度计算公式的矩阵形式时会用到$tr(AB)=tr(BA)$，证明过程见2.3。

2.3、矩阵迹交换律证明

已知$A$为$N\ast M$的矩阵，$B$为$M\ast N$的矩阵，求证$tr(AB)=tr(BA)$

证明：
$$\begin{gathered} C=AB,C为N\ast N的矩阵 \\ 其中，C_{ij}=\sum _{k=1}^M{A}_{ik}{B}_{kj} \\ 则，tr(C)=C_{11}+C_{22}+...+C_{NN}=\sum _{i=1}^N\sum _{k=1}^M{A}_{ik}{B}_{ki} \\ 其中, \\ {C}_{11}=A_{11}\ast B_{11}+A_{12}\ast B_{21}+A_{13}\ast B_{31}+...+A_{1M}\ast B_{M1} \\ {C}_{22}=A_{21}\ast B_{12}+A_{22}\ast B_{22}+A_{23}\ast B_{32}+...+A_{2M}\ast B_{M2} \\ {C}_{33}=A_{31}\ast B_{13}+A_{32}\ast B_{23}+A_{33}\ast B_{33}+...+A_{3M}\ast B_{M3} \\ ... \\ {C}_{NN}=A_{N1}\ast B_{1N}+A_{N2}\ast B_{2N}+A_{N3}\ast B_{3N}+...+A_{NM}\ast B_{MN} \\ 对以上N个式子对应项按列相加， \\ 记D_{11}=B_{11}\ast A_{11}+B_{12}\ast A_{21}+B_{13}\ast A_{31}+...+B_{1N}\ast A_{N1}=\sum _{k=1}^N{B}_{1k}{A}_{k1} \\ {D}_{22}=B_{21}\ast A_{12}+B_{22}\ast A_{22}+B_{23}\ast A_{32}+...+B_{2N}\ast A_{N2}=\sum _{k=1}^N{B}_{2k}{A}_{k2} \\ {D}_{33}=B_{31}\ast A_{13}+B_{32}\ast A_{23}+B_{33}\ast A_{33}+...+B_{3N}\ast A_{N3}=\sum _{k=1}^N{B}_{3k}{A}_{k3} \\ ... \\ {D}_{MM}=B_{M1}\ast A_{1M}+B_{M2}\ast A_{2M}+B_{M3}\ast A_{3M}+...+B_{MN}\ast A_{NM}=\sum _{k=1}^N{B}_{Mk}{A}_{kM} \\ 则tr(C)=C_{11}+C_{22}+...+C_{NN}可写作以下形式： \\ tr(C)=D_{11}+D_{22}+...+D_{MM}=\sum _{k=1}^N{B}_{1k}{A}_{k1}+\sum _{k=1}^N{B}_{2k}{A}_{k2}+\sum _{k=1}^N{B}_{3k}{A}_{k3}+...+\sum _{k=1}^N{B}_{Mk}{A}_{kM}=\sum _{i=1}^M\sum _{k=1}^N{B}_{ik}{A}_{ki}=tr(BA) \\ tr(AB)=tr(BA)得证 \end{gathered}$$
矩阵迹交换律的证明过程可参考该链接：矩阵迹的交换律以及加和号的交换律 - 鸿羽的文章 - 知乎 https://zhuanlan.zhihu.com/p/452962118

其中该性质似乎可以使用交换求和证明，即$tr(AB)={\sum }_{i=1}^N{\sum }_{k=1}^M{A}_{ik}{B}_{ki}={\sum }_{k=1}^M{\sum }_{i=1}^N{B}_{ki}{A}_{ik}=tr(BA)$

交换求和的相关内容可参考该链接：恐怖技巧——交换求和 - Westlifer的文章 - 知乎 https://zhuanlan.zhihu.com/p/499839696

三、模块度增量计算

3.1、公式推导

增量公式（子区$i、j$合并的模块度增量）：
$$\Delta Q_{ij}=e_{ij}+e_{ji}-2a_i{a}_j$$
基于模块度计算的公式（形式一）进行推导：

$Q={\sum }_i(e_{ii}-a_i^2)={\sum }_i{e}_{ii}-{\sum }_i{a}_i^2$

合并前有$R$个子区，模块度记为$Q={\sum }_{k=1}^R{e}_{kk}-{\sum }_{k=1}^R{a}_k^2$

合并子区$i、j$，假设$i<j$，则合并后有$R-1$个子区（其中第$m$个子区为合并后的子区），合并后的模块度记为$Q^{\prime }={\sum }_{k=1}^{R-1}{e}_{kk}^{\prime }-{\sum }_{k=1}^{R-1}{a^{\prime }}_k^2$

则模块度增量为：
$$\Delta Q=Q^{\prime }-Q=(\sum _{k=1}^{R-1}{e}_{kk}^{\prime }-\sum _{k=1}^R{e}_{kk})+(\sum _{k=1}^{R-1}{a^{\prime }}_k^2-\sum _{k=1}^R{a}_k^2)$$

• 计算${\sum }_{k=1}^{R-1}{e}_{kk}^{\prime }-{\sum }_{k=1}^R{e}_{kk}$

$$\begin{gathered} \sum _{k=1}^{R-1}{e}_{kk}^{\prime }=\sum _{k\ne m}{e}_{kk}^{\prime }+e_{mm}^{\prime } \\ 其中，e_{mm}^{\prime }=e_{ii}+e_{jj}+e_{ij}+e_{ji} \\ 则，\sum _{k=1}^{R-1}{e}_{kk}^{\prime }=\sum _{k=1}^R{e}_{kk}+e_{ij}+e_{ji} \\ \sum _{k=1}^{R-1}{e}_{kk}^{\prime }-\sum _{k=1}^R{e}_{kk}=e_{ij}+e_{ji} \end{gathered}$$

• 计算${\sum }_{k=1}^{R-1}{a^{\prime }}_k^2-{\sum }_{k=1}^R{a}_k^2$

$$\begin{gathered} \sum _{k=1}^{R-1}{a^{\prime }}_k^2=\sum _{k\ne m}{a^{\prime }}_k^2+{a^{\prime }}_m^2 \\ 其中，{a^{\prime }}_1^2=(e_{11}+e_{12}+...+e_{1m}+...+e_{1(R-1)}{)}^2 \\ 其中，e_{1m}=e_{1i}+e_{1j} \\ 则，{a^{\prime }}_1^2=(e_{11}+e_{12}+...+e_{1i}+...+e_{1j}+...+e_{1(R-1)}{)}^2=a_1^2 \\ \sum _{k\ne m}^{R-1}{a^{\prime }}_k^2=\sum _{k\ne i,j}^R{a}_k^2 \\ 此外，{a^{\prime }}_m^2=(e_{m1}+e_{m2}+...+e_{mm}+...+e_{m(R-1)}{)}^2= \\ [(e_{i1}+e_{i2}+...+e_{ii}+...+e_{ij}+...+e_{i(R-1)})+(e_{j1}+e_{j2}+...+e_{ji}+...+e_{jj}+...+e_{j(R-1)}){]}^2= \\ {a}_i^2+a_j^2+2\ast a_i\ast a_j \\ 则，\sum _{k=1}^{R-1}{a^{\prime }}_k^2=\sum _{k\ne m}{a^{\prime }}_k^2+{a^{\prime }}_m^2=\sum _k^R{a}_k^2+2\ast a_i\ast a_j \\ \sum _{k=1}^{R-1}{a^{\prime }}_k^2-\sum _{k=1}^R{a}_k^2=2\ast a_i\ast a_j \end{gathered}$$

• 计算$\Delta Q$

$$\begin{gathered} \Delta Q=Q^{\prime }-Q=(\sum _{k=1}^{R-1}{e}_{kk}^{\prime }-\sum _{k=1}^R{e}_{kk})+(\sum _{k=1}^{R-1}{a^{\prime }}_k^2-\sum _{k=1}^R{a}_k^2)=e_{ij}+e_{ji}-2\ast a_i\ast a_j \\ 其中，e_{ij}=e_{ji} \\ 则，\Delta Q=2\ast (e_{ij}-a_i\ast a_j) \end{gathered}$$

3.2、社区合并示例

合并前：3个子区，模块度$Q$为$-1/24$

合并后：2个子区，模块度$Q^{\prime }=2/6+1/6-(7/12{)}^2-(5/12{)}^2=-1/72$

节点	I	II	$a_i$
I	2/6	3/12	7/12
II	3/12	1/6	5/12

模块度增量：$\Delta Q=Q^{\prime }-Q=-1/72-(-1/24)=1/36$

按照公式计算：$\Delta Q=2(e_{12}-a_1\ast a_2)=2(1/12-2/12\ast 1/12)=1/36$

四、相关链接

论文：

https://arxiv.org/pdf/cond-mat/0309508.pdf

http://www.pnas.org/content/103/23/8577.full.pdf

模块度计算相关：

Modularity的计算方法——社团检测中模块度计算公式详解 | 雅乐网 (yalewoo.com)

http://t.csdn.cn/G8lnM

矩阵迹的交换律：

以及加和号的交换律 - 鸿羽的文章 - 知乎 https://zhuanlan.zhihu.com/p/452962118