摘要:置信传播(BP)算法作为推断概率图模型的主流算法是求解随机块模型中联合概率分布的重要方法之一. 但现有的方法要么在处理核边结构问题上存在精度不足问题, 要么在理论的推导上存在近似太多, 导致求解过程复杂且难以理解问题, 或两个问题均存在. 当然, 精度不足也是由近似多造成的. 导致理论近似多且推导复杂的主要原因, 是随机块模型推断过程中求解联合概率分布并不是直接套用BP算法, 即处理的图(网络)与概率图模型的图不统一. 因此, 本文利用平均场近似修正联合概率分布, 使其完全匹配BP算法的迭代公式, 这样使得在理论推导上简单易懂. 最后通过实验验证, 该方法是有效的.
关键词: 随机块模型/
置信传播算法/
联合概率分布/
平均场近似

English Abstract

全文HTML

--> --> -->

复杂系统可以用复杂网络进行建模, 而复杂网络主要由点与边组成. 在现实中, 点与边的信息是可以直接或间接^[1]获取的, 如在航空网络系统中^[2], 节点为机场, 如果机场之间具有航班, 就存在一条边; 科学家合作网络^[3], 节点表示科学家, 如果两个科学家有合作就存在一条边等. 这些简单的点与连边关系可以帮助我们揭示复杂系统的隐藏信息, 如节点的属性, 可以对应着节点在系统中属于哪个功能块或组织^[4]. 如航空网络中的核心与外围(边缘)机场^[5], 科学家合作网络中的特定研究领域的小组^[6]等. 因此, 如何通过节点的连边信息去探测这些中尺度结构, 如社团结构、核边结构(核心-边缘), 是一个重要的科学问题, 这些问题得到了来自各个领域****的研究.
社团结构的描述为社团内部紧密相连, 而社团之间连接稀疏^[6]. 社团结构的探测方法已经有很多, 主要有基于模块度的方法^[7-9]、谱划分的方法^[4,10]、矩阵分解的方法^[11,12]、基于动力学的方法^[13,14]、随机块模型的方法^[15,16]以及其他方法^[17,18]; 核边结构的描述为核与核紧密相连, 核与边紧密相连, 边与边连接稀疏^[19,20]. 相较于社团结构, 对核边结构探测方法的研究相对较少. 主要有给定某种中心性指标, 然后通过截断选取最重要的一部分节点作为核心节点^[21,22]; 定义核边结构的指标, 然后利用最优化方法求解^[23,24]; 通过随机块模型探测核边结构^[25,26]以及其他方法^[5,27-29].
可以看出, 随机块模型在解决社团结构与核边结构的探测问题上都是行之有效的^[25]. 随机块模型是一种生成模型, 可以生成具有模块结构的网络^[15]. 即首先给定n个节点, 随机分配到k个模块中. 假设向量$ {\boldsymbol {\gamma}} $表示每个模块的规模, 即$ \gamma_r $表示第r个模块的大小占总节点个数的比例, 则每个节点就以$ \gamma_r $的概率分配到第r模块中. 然后给定一个$ k\times k $的关联矩阵p (对称矩阵), 任意两个节点(一个节点属于模块s, 一个节点属于模块r)以$ p_{rs} $的概率相连. 这样就可以得到一个以参数$ {\boldsymbol {\gamma}} $与p生成的随机网络. 如, 当$ k = 2 $时, 如果$ p_{11} > p_{12}, \;p_{22} > $$ p_{12} $, 生成的网络具有社团结构; 如果$p_{11}\geqslant p_{12} \geqslant p_{22}$, 生成的网络具有核边结构.
给定特定参数$ {\boldsymbol {\gamma}} $与p, 可以生成具有社团结构或核边结构的网络. 也就是说, 如果一个网络具有社团或核边结构, 那么可以通过学习该模型(即模型参数), 来确定它是怎么生成的, 以此来探测哪些节点属于哪个社团, 或者哪些节点属于核, 哪些节点属于边. 因此, 随机块模型不仅可以用来探测社团结构, 还可以用来探测核心边缘结构. 其中EM算法是学习该模型的重要方法^{[16,25,26,30,31]}, 而通过联合概率公式, 确定边际概率又是EM算法的关键部分与难点. 基于马尔科夫链的蒙特卡罗方法(MCMC)是解决该难点的一种较好的方法^[30], 但是由于采样空间非常大, 会造成时间复杂度过高、精度低的问题. 置信传播(BP)算法^[32]可以很好地解决上述问题^{[16,25,26,31]}, 因此被广泛用来解决社团或核边结构的探测问题. BP算法用来推断马尔科夫随机场中的边际概率, 当马尔可夫网络是树状图时, 具有精确解; 当马尔可夫网络是稀疏的, 会得到近似精确解^[33].
应用BP算法解决社团结构的探测问题时, 由于网络的稀疏性, 在公式的推导中关于节点邻居可以忽略不计的近似是合理的^[16,31], 但是处理核边结构探测问题时, 这种近似却是致命的, 因为核节点具有天然的大度性质. 因此, 在另外一篇文章中, 我们对BP算法进行了稍微的改进(IBP), 就可以很好地解决此问题^[25]. 但是伴随而来的就是IBP算法的推导过程中, 需要进行大量的近似处理, 以及各种模糊的推导, 并不是可以套用BP算法公式^[34]直接得到的. 同样的道理, BP算法解决社团结构的探测问题时, 即使其中各种近似是合理的, 也会存在上述不能套用BP算法公式直接得到的情况. 这使得即使BP或IBP算法在社团划分或核边结构探测中能得到很好的结果, 但是推导过程却是复杂且难于理解, 或者说是不那么精确的. 因此, 本文从平均场近似^[35]的角度出发, 先应用平均场近似对联合概率进行处理, 得到可以完全匹配BP算法的形式, 然后再套用BP算法公式(而不是以往的先套用公式, 然后在公式里近似处理), 在不需要任何修改、近似以及假设的情况下直接得到求解边际概率的迭代公式. 使得随机块模型解决社团结构或核边结构探测问题时, 既可以保证结果的精度, 又可以保证理论的简单明了、有理有据.

随机块模型最直接的用处就是按照某种需求生成具有某种模块结构的人工网络. 例如在社团探测中, 就可以用随机块模型生成具有社团结构的人工网络, 以此来评估社团挖掘算法的好坏. 但是本文将用到它的另外一个重要应用, 就是通过学习随机块模型中的参数使其最好地匹配给定网络, 从而揭示给定网络的结构. 显然, 这是一个参数估计问题, 可以采用最大似然估计的方法进行解决^[16,25,26].
2

-->

2.1.最大似然估计

最大似然估计研究的是给定什么样的参数可以使得似然函数的值最大, 也就是使得已经发生的事件概率最大. 假设给定一个无权无向网络的邻接矩阵A, 其中$ A_{ij} = 1 $表示第i个节点与第j个节点相连, $ A_{ij} = 0 $表示不相连. 在随机块模型中, 当给定参数p与$\boldsymbol {\gamma}$, 则似然函数就可以用$ P({\boldsymbol{A}}|{\boldsymbol{p}}, {\boldsymbol{\gamma}}) $来表示. 因此, 需要最大化$ P({\boldsymbol{A}}|{\boldsymbol{p}}, {\boldsymbol{\gamma}}) $, 以此求解p与$ {\boldsymbol{\gamma}} $, 从而来揭示网络的结构. 但是, 通过观察可以发现, 似然函数里包含一个隐形变量g, 其中$ g_i $表示节点i所属的模块. 因此有

$ \begin{split}\;& P({\boldsymbol{A}}|{\boldsymbol{p}},{\boldsymbol{\gamma}}) = \sum\limits_{\boldsymbol{g}}P({\boldsymbol{A}},{\boldsymbol{g}}|{\boldsymbol{p}},{\boldsymbol{\gamma}}) \\ =\;& \sum\limits_{\boldsymbol{g}}P({\boldsymbol{A}}|{\boldsymbol{g}},{\boldsymbol{p}},{\boldsymbol{\gamma}})P({\boldsymbol{g}}|{\boldsymbol{\gamma}})\\ =\;& \sum\limits_{\boldsymbol{g}}\left[\prod\limits_{i < j}p_{g_ig_j}^{A_{ij}}(1-p_{g_ig_j})^{1-A_{ij}}\prod\limits_i\gamma_{g_i}\right]. \end{split} $

可以看出, 最大化公式(1)是一个含有隐变量的最大似然估计问题, 因此可以应用EM算法进行求解^[36].
2 -->

2.2.EM算法求解随机块模型

EM算法在理论推导上是根据詹森不等式不断求解对数似然函数的下界函数最优值, 以此逐渐收敛到原始函数的最优值. 因此, 对(1)式取对数, 然后根据詹森不等式有^[25,26]:

$\begin{split} &\log \sum\limits_{\boldsymbol{g}} {P\left( {{\boldsymbol{A}},{\boldsymbol{g}}|{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right)} \\\geqslant\;&\sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\log \frac{{P\left( {{\boldsymbol{A}},{\boldsymbol{g}}|{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right)}}{{q\left( {\boldsymbol{g}} \right)}}} \overset{\Delta} = L\left( {{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right), \end{split}$

其中

$\begin{split} q\left( {\boldsymbol{g}} \right) =\;& \frac{{P\left( {{\boldsymbol{A}},{\boldsymbol{g}}|{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right)}}{{\displaystyle\sum\limits_{\boldsymbol{g}} {P\left( {{\boldsymbol{A}},{\boldsymbol{g}}|{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right)} }}\\ =\;& \frac{\displaystyle{\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j}}}} \right)}^{1 - {A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}} }}{{\displaystyle\sum\limits_{\boldsymbol{g}} {\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j}}}} \right)}^{1 - {A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}} } }}. \end{split}$

可以看出, $ q({\boldsymbol{g}}) = P({\boldsymbol{g}}|{\boldsymbol{A}}, {\boldsymbol{p}}, {\boldsymbol{\gamma}}) $, 表示给定一个网络, 在参数p与$ {\boldsymbol{\gamma}} $给定的条件下, 该网络每个节点分配为g的概率, 这是一个联合概率.
当(3)式成立的情况下, (2)式等号成立, 所以最大化公式(1)问题可以转化为最大化$L\left( {\boldsymbol{p}}, {\boldsymbol{\gamma}} \right)$问题. 对$ L\left( {{\boldsymbol{p}}, {\boldsymbol{\gamma}} } \right) $进一步推导有:

$ \begin{split}\;& L\left( {{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right) = \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\log \frac{{P\left( {{\boldsymbol{A}},{\boldsymbol{g}}|{\boldsymbol{p}},{\boldsymbol{\gamma}} } \right)}}{{q\left( {\boldsymbol{g}} \right)}}} \\ =\;& \sum\limits_{\boldsymbol{g}} q\left( {\boldsymbol{g}} \right)\log \Bigg[ \prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j}}}} \right)}^{1 - {A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}} \Bigg]\\& - \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\log q\left( {\boldsymbol{g}} \right)} \\ =\;& \sum\limits_{\boldsymbol{g}} q\left( {\boldsymbol{g}} \right)\Bigg\{ \sum\limits_{i < j} \Big[ {A_{ij}}\log {p_{{g_i}{g_j}}} + \left( {1 - {A_{ij}}} \right) \\&\times\log \left( {1 - {p_{{g_i}{g_j}}}} \right) \Big] + \sum\limits_i {\log {\gamma _{{g_i}}}} \Bigg\}\\ &- \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\log q\left( {\boldsymbol{g}} \right)}\\ =\;& \frac{1}{2}\sum\limits_{i \ne j} \sum\limits_{{\boldsymbol{rs}}} \Big[ {A_{ij}}q_{{\boldsymbol{rs}}}^{ij}\log {p_{{\boldsymbol{rs}}}} + \left( {1 - {A_{ij}}} \right)\\ &\times q_{{\boldsymbol{rs}}}^{ij}\log \left( {1 - {p_{{\boldsymbol{rs}}}}} \right) \Big] \\& + \sum\limits_{i{\boldsymbol{r}}} {q_{\boldsymbol{r}}^i\log {\gamma _{\boldsymbol{r}}}}- \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\log q\left( {\boldsymbol{g}} \right)}, \\[-15pt]\end{split} $

其中$ q({\boldsymbol{g}}) $的边际概率$ q^i_{\boldsymbol{r}} $表示节点i属于群组r的概率:

$ q_{\boldsymbol{r}}^i = \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right){\delta _{{g_i},{\boldsymbol{r}}}}}, $

$ q({\boldsymbol{g}}) $的边际概率$ q^{ij}_{{\boldsymbol{rs}}} $表示节点i属于群组r的概率且节点j属于群组s的概率:

$ q_{{\boldsymbol{rs}}}^{ij} = \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right){\delta _{{g_i},{\boldsymbol{r}}}}{\delta _{{g_j},{\boldsymbol{s}}}}}. $

如果$ g_i = {\boldsymbol{r}} $, 则$ \delta _{{g_i}, {\boldsymbol{r}}} = 1 $, 否则$ \delta _{{g_i}, {\boldsymbol{r}}} = 0 $.
对(4)式求解最大值(满足约束条件$\displaystyle\sum\nolimits_{\boldsymbol{r}} {{\gamma _{\boldsymbol{r}}}} = $$ 1$), 应用拉格朗日乘子法可以得到:

$ {p_{{\boldsymbol{rs}}}} = \frac{{\displaystyle\sum\limits_{i \ne j} {{A_{ij}}q_{{\boldsymbol{rs}}}^{ij}} }}{{\displaystyle\sum\limits_{i \ne j} {q_{{\boldsymbol{rs}}}^{ij}} }}, $

$ {\gamma _{\boldsymbol{r}}} = \frac{1}{n}\sum\limits_i {q_{\boldsymbol{r}}^i} . $

根据上述推导过程, 可以得到一个迭代公式. 即首先初始化p与$ {\boldsymbol{\gamma}} $, 根据(3)式计算边际概率$ q^i_{\boldsymbol{r}} $与$ q^{ij}_{{\boldsymbol{rs}}} $, 然后根据(7)式与(8)式计算新的p与$ {\boldsymbol{\gamma}} $, 重复上述过程直到收敛为止. 最后可以通过收敛的边际概率$ q^i_{\boldsymbol{r}} $的值对每个节点进行划分. 其中, 当需要探测社团结构时, 只需初始化p为具有社团结构形式的关联矩阵, 如果需要探测核边结构时, 则需要初始化p为具有核边结构形式的关联矩阵.
需要注意的是, 在应用(7)式求解$ p_{{\boldsymbol{rs}}} $时, 需要计算任意节点对之间的$ q^{ij}_{{\boldsymbol{rs}}} $值, 这将花费大量的计算时间. 因此可以化简为

$\begin{split} \sum\limits_{ij} {q_{{\boldsymbol{rs}}}^{ij}} =\;& \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\sum\limits_i {{\delta _{{g_i},{\boldsymbol{r}}}}} \sum\limits_j {{\delta _{{g_j},{\boldsymbol{s}}}}} } \\=\;& \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right){n_{\boldsymbol{r}}}{n_{\boldsymbol{s}}}} = \left\langle {{n_{\boldsymbol{r}}}{n_{\boldsymbol{s}}}} \right\rangle \approx \left\langle {{n_{\boldsymbol{r}}}} \right\rangle \left\langle {{n_{\boldsymbol{s}}}} \right\rangle, \end{split}$

其中$n_{\boldsymbol{r}} = \displaystyle\sum\nolimits_i {{\delta _{{g_i}, {\boldsymbol{r}}}}}$, $ \left\langle \cdots \right\rangle $表示期望. 于是有

$\begin{split} \left\langle {{n_{\boldsymbol{r}}}} \right\rangle =\;& \sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right)\sum\limits_i {{\delta _{{g_i},{\boldsymbol{r}}}}} } \\= \;&\sum\limits_i {\sum\limits_{\boldsymbol{g}} {q\left( {\boldsymbol{g}} \right){\delta _{{g_i},{\boldsymbol{r}}}}} } = \sum\limits_i {q_r^i}. \end{split}$

所以(7)式可以写成:

$ {p_{{\boldsymbol{rs}}}} = \frac{{\displaystyle\sum\limits_{i \ne j} {{A_{ij}}q_{{\boldsymbol{rs}}}^{ij}} }}{{\displaystyle\sum\limits_i {q_{\boldsymbol{r}}^i} \displaystyle\sum\limits_j {q_{\boldsymbol{s}}^j} }}. $

应用(11)式代替(7)式只需要计算有边的节点对的两节点边际概率, 这将很大程度上减少计算量, 特别是在稀疏网络. 而且在BP算法中这种处理更显得尤为必要, 后面将进一步介绍.
到这里并没结束, 通过遍历所有构型($ k^n $种构型)的方法根据(3)式求解全概率$ q({\boldsymbol{g}}) $以此求解边际概率$ q^i_{\boldsymbol{r}} $与$ q^{ij}_{{\boldsymbol{rs}}} $是显然不可行的. 因此如何根据(3)式求解边际概率$ q^i_{\boldsymbol{r}} $与$ q^{ij}_{{\boldsymbol{rs}}} $是EM算法的难点, 也是关键部分. 下面将介绍两种求解方法: MC采样和BP算法, 并分析他们在处理社团与核边结构探测中存在的不足之处.
2 -->

2.3.马尔科夫链蒙特卡罗方法(MCMC)

这里应用MCMC模拟的方法按照联合概率公式(3)对网络的划分构型(即$g = [g_1, g_2, \cdots, g_n]^{\rm T}$)进行采样. 本文将采用单分量 Metropolis-Hastings (M-H) 方法^[37]. 设$g^t = [g_1^t,\; g_2^t, \;\cdots, \;g_n^t]^{\rm T}$表示在马尔科夫链中第t轮迭代后的状态, 其中每一轮迭代都需要对每个分量的状态进行一次迭代更新.
例如, 在$ t+1 $轮迭代中, 已经迭代好的分量(第i个分量)的状态记为$ g_i^{t+1} $. 假设要对第j个分量进行采样, 首先设$ g_{-j}^{t+1} $表示在$ t+1 $轮迭代中在对第j个分量进行处理时, 其他分量的状态. 特别地, 当对每一轮的各个分量进行迭代是按照顺序处理时, 有:

$ g_{-j}^{t+1} = [g_1^{t+1},\;g_2^{t+1},\;\cdots,\;g_{j-1}^{t+1},\;g_{j+1}^{t},\;\cdots,\;g_n^t]^{\rm T}. $

把第j个分量的状态$ g_j^{t} $随机变为一种状态$ \hat{g}_j^{t+1} $, 因为有k个模块, 所以建议转移概率$ q(g_j^{t}, \hat{g}_j^{t+1}) = $$ q(\hat{g}_j^{t+1}, g_j^{t}) = 1/k $. 根据(3)式, 按照 M-H 规则接收这种状态(从状态$ g_j^{t} $变为状态$ \hat{g}_j^{t+1} $)的概率:

$\begin{split} &\alpha(g_j^{t},\hat{g}_j^{t+1})= \\ \;& {\rm {min}}\left\{1,\frac{\displaystyle\bigg[{\prod\limits_{i\neq j} {p_{{g_i}{g_j^t}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j^t}}}} \right)}^{1 - {A_{ij}}}}}\bigg] {\gamma _{g_j^t}} }} {\displaystyle \bigg[\prod\limits_{i\neq j} {p_{{g_i}{\hat{g}_j^{t+1}}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{\hat{g}_j^{t+1}}}}} \right)}^{1 - {A_{ij}}}}}\bigg] {\gamma _{\hat{g}_j^{t+1}}}}\right\}, \end{split}$

式中$ g_i(i\neq j) $为$ g_{-j}^{t+1} $中分别对应的量. 在区间$ (0, 1) $内按均匀分布取出一个数记为$ \mu $, 如果$\mu\leqslant $$ \alpha(g_j^{t}, \hat{g}_j^{t+1})$, 则$ {g}_j^{t+1} = \hat{g}_j^{t+1} $, 否则$ {g}_j^{t+1} = {g}_j^{t} $. 当每一轮所有结点的状态更新完以后(也可以不按顺序更新), 再进行下一轮的更新. 在更新的轮数足够大的情况下, 采样的分布收敛到(3)式的概率分布, 因此可以通过统计估计出边际概率$ q^i_{\boldsymbol{r}} $与$ q^{ij}_{{\boldsymbol{rs}}} $.
但是, 由于所有构型共用$ k^n $种可能, 在网络较大的情况下, 采样空间巨大, 很难收敛到(3)式的概率分布, 从而造成计算结果产生误差, 并且不稳定. 下面介绍一种快速求解边际概率的方法, 即BP算法.
2 -->

2.4.BP算法

3 -->

2.4.1.BP算法简介

BP算法可以用来推断马尔科夫随机场(又称为概率图模型)中的边际概率^[33]. 在本文中, 考虑一种特殊的概率图模型—成对图模型. 成对图模型定义在一个无向图$G({\boldsymbol{V}}, {\boldsymbol{E}})$上, 每一个顶点$v_i\in {\boldsymbol{V }}$表示X中的一个随机变量$ x_i $, 每一条边$(v_i, v_j)\in {\boldsymbol{E}}$表示随机变量$ x_i $与$ x_j $之间的某种依赖关系. 给定每个节点$ v_i $一个点位势$ \phi_i(x_i) $, 以及每一条边$(v_i, v_j)\in {\boldsymbol{E}}$一个边位势$ \psi_{ij}(x_i, x_j) $, 则在该成对马尔科夫网上的分布为

$ P\left( {\boldsymbol{X}} \right) = \frac{1}{Z}\mathop \prod \limits_{\left( {i,j} \right) \in{\boldsymbol{ E}}} {\psi _{ij}}\left( {{x_i},{x_j}} \right)\mathop \prod \limits_{i \in {\boldsymbol{V}}} {\phi _i}\left( {{x_i}} \right), $

其中Z为归一化常数, 通常称为配分函数.
求解(14)式这个联合概率的边际概率可以通过BP算法迭代公式求得, 迭代过程如下^[34]:

$\begin{split} &{v_{i \to j}}\left( {{x_i}} \right) \propto {\phi _i}\left( {{x_i}} \right)\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{{x_k}} {{\psi _{ik}}\left( {{x_i},{x_k}} \right){v_{k \to i}}\left( {{x_k}} \right)} }, \end{split}$

$\begin{split} &{b_i}\left( {{x_i}} \right) \propto {\phi _i}\left( {{x_i}} \right)\prod\limits_{k \in N\left( i \right)} {\sum\limits_{{x_k}} {{\psi _{ik}}\left( {{x_i},{x_k}} \right){v_{k \to i}}\left( {{x_k}} \right)} }, \end{split}$

其中, $ N(i) $表示节点$ v_i $的邻居, $ v_{i \to j}\left( {{x_i}} \right) $为概率图修改后(即删除节点$ v_j $) $ x_i $的边际概率, $ b_i(x_i) $为$ x_i $的边际概率. 另外还可以得到一条边上两个节点的联合概率:

$ {b_{ij}}\left( {{x_i},{x_j}} \right) \propto {\psi _{ij}}\left( {{x_i},{x_j}} \right){v_{i \to j}}\left( {{x_i}} \right){v_{j \to i}}\left( {{x_j}} \right). $

3 -->

2.4.2.BP算法求解随机块模型

对于联合概率公式(3), 可以看成是一个马尔科夫网上的概率分布. 在这里, 马尔科夫网是一个全连通图, 且随机变量为$g = [g_1,\; g_2, \;\cdots, \;g_n]^{\rm T}$. 对于一对节点$ (v_i, v_j) $, 如果$ A_{ij} = 1 $, 其边位势为$ \psi_{ij}(g_i, g_j) = $$ P_{g_ig_j} $, 如果$ A_{ij} = 0 $, 其边位势为$ \psi_{ij}(g_i, g_j) = 1- $$ P_{g_ig_j} $, 节点$ v_i $的点位势为$ \phi_i(g_i) = \gamma_{g_i} $. 所以根据BP 算法的迭代公式(15), 当j节点是i节点的邻居时, 可以得到^[26,31]:

$ \eta _{\boldsymbol{r}}^{i \to j} = \frac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_{i \to j}}}}\prod\limits_{k \in {\boldsymbol{V}}/{N^ * }\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\left( {1 - {p_{{\boldsymbol{rs}}}}} \right)} }\prod\limits_{k \in {\boldsymbol{N}}\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }, $

其中“消息”$ \eta _{\boldsymbol{r}}^{i \to j} $对应于(15)式中的$ v_{i \to j}( {{x_i}} ) $, 在这里可以理解为移除节点$ v_j $, $ v_i $属于模块r的概率. $ Z_{i \to j} $表示配分函数, $ N(i) $表示节点$ v_i $的邻居, $ N^*(i) = N(i)\bigcup \{i\} $. 为了更方便更快捷地计算, 需要做两个假设, 会得到两个近似.
假设1　如果第j个节点不是第i节点的邻居, 移除节点$ v_j $对节点$ v_i $属于模块r的概率没有影响, 即有 $ \eta _{\boldsymbol{s}}^{k \to i} \approx q^k_{\boldsymbol{s}} $;
假设2　当网络规模很大且很稀疏的时候, 对集合$ {\boldsymbol{V}}/{N^ * } $中所有元素做运算近似于对集合V中所有元素做运算, 即有$ f({\boldsymbol{V}}/{N^ * } )\approx f({\boldsymbol{V}}) $, 其中$ f(*) $代表某种运算.
根据上述两个假设, 有:

$\begin{split}& \prod\limits_{k \in {\boldsymbol{V}}/{N^ * }\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\left( {1 - {p_{{\boldsymbol{rs}}}}} \right)} } \\\approx\,& \prod\limits_{k \in {\boldsymbol{V}}/{N^ * }\left( i \right)}{\sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k\left( {1 - {p_{{\boldsymbol{rs}}}}} \right)} } \\\approx\,& \prod\limits_{k \in {\boldsymbol{V}}} {\sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k\left( {1 - {p_{{\boldsymbol{rs}}}}} \right)} } \\= \,& \prod\limits_{k \in {\boldsymbol{V}}} {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)}, \end{split}$

这是一个外场量. 所以(18)式可以写成:

$\begin{split} \eta _{\boldsymbol{r}}^{i \to j} =\;& \dfrac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_{i \to j}}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \\ &\times\prod\limits_{k \in N\left( i \right)/j}{\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }, \end{split}$

其中配分函数为

$\begin{split} {Z_{i \to j}} =\;& \sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \\ &\times\prod\limits_{k \in N\left( i \right)/j}{\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }. \end{split}$

根据(16)式可以得到边际分布$ q^i_{\boldsymbol{r}} $:

$ q_{\boldsymbol{r}}^i = \dfrac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_i}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \prod\limits_{k \in N\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }, $

其中配分函数为

$ {Z_i} = \sum\limits_{\boldsymbol{r}} {{\gamma _{\boldsymbol{r}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \prod\limits_{k \in N\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} } }. \;\; \;\;\qquad$

根据(17)式可以得到一条边上两个点的联合概率:

$ q_{{\boldsymbol{rs}}}^{ij} = \frac{{\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}{p_{{\boldsymbol{rs}}}}}}{{\displaystyle\sum\limits_{{\boldsymbol{rs}}} {\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}{p_{{\boldsymbol{rs}}}}} }}. $

通过不停迭代(20)式和(22)式, 直到收敛, 然后通过(22)式与(24)式即可计算出$ q^i_{\boldsymbol{r}} $与$ q^{ij}_{{\boldsymbol{rs}}} $.
上述算法在网络规模比较大且稀疏的时候是适用的, 但是应用此方法探测核边结构, 就会存在问题. 众所周知, 一个具有核边结构的网络的核节点不仅和核节点紧密相连, 而且和边节点紧密相连, 也就是说核节点与整个网络中的节点都是紧密相连的, 具有天然的大度性质. 所以在假设2中, 对集合$ {\boldsymbol{V}}/{N^ * } $中所有元素做运算近似于对集合V中所有元素做运算, 这种假设在核边结构的探测中是不合理的, 有可能是错误的. 针对于此, 我们在另外一篇文章中对BP 算法进行了修正^[25], 可以很好地解决这个问题.
3 -->

2.4.3.修正后的BP算法(IBP)

首先对联合概率分布公式(3)稍作处理, 可得:

$\begin{split} \;& q\left( {\boldsymbol{g}} \right) = \frac{{\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j}}}} \right)}^{1 - {A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}} }}{{\displaystyle\sum\limits_{\boldsymbol{g}} {\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j}}}} \right)}^{1 - {A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}} } }} \\=\;& \frac{{\prod\limits_{i < j} {{{\left( {\frac{{{p_{{g_i}{g_j}}}}}{{1 - {p_{{g_i}{g_j}}}}}} \right)}^{{A_{ij}}}}\left( {1 - {p_{{g_i}{g_j}}}} \right)} \prod\limits_i {{\gamma _{{g_i}}}} }}{{\displaystyle\sum\limits_{\boldsymbol{g}} {\prod\limits_{i < j} {{{\left( {\frac{{{p_{{g_i}{g_j}}}}}{{1 - {p_{{g_i}{g_j}}}}}} \right)}^{{A_{ij}}}}\left( {1 - {p_{{g_i}{g_j}}}} \right)} \prod\limits_i {{\gamma _{{g_i}}}} } }}, \end{split}$

则每条边的信使为^[25,38]

$\begin{split} \eta _{\boldsymbol{r}}^{i \to j} =\;& \dfrac{{{\gamma _r}}}{{{Z_{i \to j}}}}\prod\limits_{k \in {\boldsymbol{V}}/\left\{ {i,j} \right\}} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\left( {1 - {p_{{\boldsymbol{rs}}}}} \right)} } \\&\times\prod\limits_{k \in N\left( i \right)/j}{\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\dfrac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }. \end{split}$

给定一个假设:
假设3　当网络很大且很稀疏的时候, 对集合$ {\boldsymbol{V}}/\{i, j\} $中所有元素做运算近似于对集合V中所有元素做运算. 显然, 这个假设对结果的影响可以忽略不计.
根据假设1与假设3可以得到下面近似:

$\begin{split} &\prod\limits_{k \in {\boldsymbol{V}}/\left\{ {i,j} \right\}} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\left( {1 - {p_{{\boldsymbol{rs}}}}} \right)} } \approx \\&\prod\limits_{k \in {\boldsymbol{V}}/\{ {i,j} \}}{\bigg(1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} \bigg)} \approx \prod\limits_{k \in {\boldsymbol{V}}} {\bigg( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \bigg)}, \end{split}$

这是一个与(19)式一样的外场量, 但是这里的近似与(19)式的近似相比, 也就是假设3与假设2相比不会依赖节点度的大小, 所以更适用于具有核边结构的网络. (26)式可以写成:

$\begin{split} \eta _{\boldsymbol{r}}^{i \to j} =\;& \frac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_{i \to j}}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \\&\times\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{rs}}}}} }, \end{split}$

其中配分函数为

$\begin{split} {Z_{i \to j}} = \;&\sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \\ &\times\prod\limits_{k \in N\left( i \right)/j}{\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }, \end{split}$

则边际分布$q^i_{\boldsymbol{r}}$为^[25,38]

$\begin{split} q_{\boldsymbol{r}}^i =\;& \frac{{{\gamma _r}}}{{{Z_i}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)}\\ &\times\prod\limits_{k \in N\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{rs}}}}} }, \end{split}$

其中

$\begin{split} {Z_i} =\;& \sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\prod\limits_k {\left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} \\ &\times\prod\limits_{k \in N\left( i \right)}{\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }. \end{split}$

根据(17)式可以得到一条边上两个点的联合概率:

$ q_{{\boldsymbol{rs}}}^{ij} = \frac{{\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}\dfrac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}}}{{\displaystyle\sum\limits_{{\boldsymbol{rs}}} {\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }}. $

通过(28)式, (30)式, (32)式与(20)式, (22)式, (24)相比, 可以发现算法改进后, 最终的结果仅仅是$ \dfrac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}} $代替了$ p_{{\boldsymbol{rs}}} $, 没有增加任何多余的计算复杂度, 但是修改后的方法已经验证了更适用于含有大度节点的网络^[25], 如核边结构网络.
改进后的算法虽然在结果上表现出来强大的优势, 但是在理论推导上存在一定的问题. 根据BP算法的计算信使公式(15), 从联合概率公式(25)并不能直接推导得到的信使迭代公式(26)(这种直接写法是参考文献[35]). 这是因为概率分布公式(25)相比(3)式虽然做了变形处理, 但是两个公式在数值上依然相等, 即在把一个全连通图当做一个马尔科夫网的情况下, 对于一对节点$ (v_i, v_j) $, 如果$ A_{ij} = 1 $, 其边位势还是$ \psi_{ij}(g_i, g_j) = P_{g_ig_j} $, 如果$ A_{ij} = 0 $, 其边位势也还是$\psi_{ij}(g_i, g_j) = 1- P_{g_ig_j}$. 所以说这里得到的信使迭代公式应该还是(18)式而不是(26)式.
在这里如果想从概率分布公式(25)得到相比信使迭代公式(18)更好的(26)式, 换句话说, 如果想得到比BP算法更好的IBP 算法, 可以按照下面的思路去理解. 这里的马尔科夫网不仅仅是一个全连通图, 还需要在全连通图上再加上M (原网络上边的个数)条额外边. 对于一对节点$ (v_i, v_j) $, 如果$ A_{ij} = 1 $, 马尔科夫网上的边就要看成两条边, 其边位势就分解成$ \psi_{ij}(g_i, g_j) = 1-P_{g_ig_j} $与$\psi_{ij}(g_i, g_j) = \dfrac{P_{g_ig_j}}{1-P_{g_ig_j}}$, 如果$ A_{ij} = 0 $, 其边位势还是$ \psi_{ij}(g_i, g_j) = 1-P_{g_ig_j} $. 这样理解以后, 就可以按照BP算法公式(15)得到信使迭代公式(26). 即使这样, 根据假设3与假设2还是得不到最终的信使迭代公式(28), 在这里还应满足相比假设2更严格的假设, 即: 在马尔科夫网中所有边位势为$ \psi_{ij}(g_i, g_j) = {1-P_{g_ig_j}} $的节点对$ {(v_i, v_j)} $(包括原网络中不存在边以及额外增加的存在边), 有$ q^i_{\boldsymbol{r}} \approx \eta _{\boldsymbol{r}}^{i \to j} $. 意思就是存在一些少部分节点对$ {(v_i, v_j)} $, 即使$ v_j $是$ v_i $的邻居也要近似处理, 即移除节点$ v_j $对节点$ v_i $属于模块r的概率的影响忽略不计. 总的来说, 虽然改进后的BP算法在实验结果上会表现出很好的效果, 但是存在理论推导上比较复杂、近似处理的过程太多、不易理解等问题.
3 -->

2.4.4.小　结

BP算法可以对联合概率公式(3)与(25)进行求解, 以此得到两种迭代公式(前者用BP算法表示, 后者用IBP算法表示)求解边际概率. 但是都存在一些问题: 1)在BP 算法中, 要满足假设1与假设2, 因此不适用核边结构的探测; 2)在IBP算法中, 虽然在实验中已经验证了比BP算法有着更好的性能, 但是在理论上处理得太过复杂, 还存在一些在理论上不易理解的假设.
而且, 通过观察可以看出, 无论是BP算法还是IBP算法所处理的马尔科夫网(至少也是全连通图)都是与要处理的原始网络不一致的. 虽然原始网络是稀疏的, 但是马尔科夫网络却是全连通的, 这也造成了即使BP算法与IBP算法表现得效果好, 却有所疑惑. 这是因为我们都知道: 当马尔科夫网是树状网络, BP 算法是精确的, 当马尔科夫网是稀疏的时候, BP算法是近似精确的.
所以, 下面将从另外一个角度出发, 即首先对联合概率公式(3)与(25)进行近似处理(而不是在BP算法的推理中近似处理), 得到一个新的概率公式, 这个概率公式要满足两个条件: 1)由这个公式导出的马尔科夫随机网与原始网络一致; 2)这个公式在形式上与(14)式完全一致. 前者可以保证马尔科夫网的稀疏性, 后者可以保证在应用BP算法的时候简洁明了且便于理解.

首先对联合概率公式(3)进行平均场近似处理. 概率分布(3)式可以表示为

$\begin{split} &q\left( {\boldsymbol{g}} \right) \propto \prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}{{\left( {1 - {p_{{g_i}{g_j}}}} \right)}^{1 - {A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}} =\\ &\bigg[ {\exp \bigg( {\sum\limits_{i < j} {( {1 - {A_{ij}}} )\ln ( {1 - {p_{{g_i}{g_j}}}} )} } \bigg)} \bigg]\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}} \prod\limits_i {{\gamma _{{g_i}}}}. \end{split}$

令${h_i} = \displaystyle\sum\nolimits_{j \ne i} {\left( {1 - {A_{ij}}} \right)\ln \left( {1 - {p_{{g_i}{g_j}}}} \right)}$, 通过平均场近似有

$\begin{split} h_i^{\rm MF} =\;& \sum\limits_{j \ne i} {\sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^j\left( {1 - {A_{ij}}} \right)\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)} }\\ \approx\;& \sum\limits_{\boldsymbol{s}} {\left[ {\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)\sum\limits_j {q_{\boldsymbol{s}}^j\left( {1 - {A_{ij}}} \right)} } \right]}, \end{split}$

再应用一次平均场近似, 即

$ \sum\limits_j {q_{\boldsymbol{s}}^j{A_{ij}}} \approx \frac{{{d_i}}}{n}\sum\limits_j {q_{\boldsymbol{s}}^j}, $

其中${d_i} = \displaystyle\sum\nolimits_j {{A_{ij}}}$, 表示节点$ v_i $的度. 则(34)式可以写成:

$ h_i^{\rm MF} = \sum\limits_{\boldsymbol{s}} {\left[ {\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)\left( {1 - \dfrac{{{d_i}}}{n}} \right)\sum\limits_j {q_{\boldsymbol{s}}^j} } \right]}. $

因此联合概率公式(3)经过两次平均场近似可以写成:

$\begin{split}&q\left( {\boldsymbol{g}} \right) \propto \left[ {\exp \left( {\sum\limits_i {h_i^{\rm MF}} } \right)} \right]\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}} \prod\limits_i {{\gamma _{{g_i}}}}\\=\;& \prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}} \prod\limits_i \Bigg\{ \exp \Bigg[ \sum\limits_{\boldsymbol{s}} \ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)\\&\times\left( {1 - \dfrac{{{d_i}}}{n}} \right)\sum\limits_j {q_{\boldsymbol{s}}^j} \Bigg]{\gamma _{{g_i}}} \Bigg\}. \end{split}$

把此联合概率分布定义为一个成对马尔科夫网上的概率分布, 而且马尔科夫网与原始网络一致. 其中, 对于一条边$ (v_i, v_j) $, 其边位势为$ \psi_{ij}(g_i, $$ g_j) = P_{g_ig_j} $, 对于一个节点$ v_i $, 其点位势为$\phi_i(g_i) = $$ {\exp \left[ {\displaystyle\sum\nolimits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)\left( {1 - \dfrac{{{d_i}}}{n}} \right)\displaystyle\sum\nolimits_j {q_{\boldsymbol{s}}^j} } } \right]{\gamma _{{g_i}}}}$. 所以直接根据BP算法的公式(15)—(17)式, 可以不经过任何近似地得到:

$\begin{split} \eta _{\boldsymbol{r}}^{i \to j} = \;&\frac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_{i \to j}}}}\exp \Bigg( \sum\limits_{\boldsymbol{s}} \ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\left( {1 - \frac{{{d_i}}}{n}} \right)\\&\times\sum\limits_k {q_{\boldsymbol{s}}^k} \Bigg)\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }, \\[-15pt]\end{split}$

$\begin{split} \eta _{\boldsymbol{r}}^{i \to j} = \;&\frac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_{i \to j}}}}\exp \Bigg( \sum\limits_{\boldsymbol{s}} \ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\left( {1 - \frac{{{d_i}}}{n}} \right)\\&\times\sum\limits_k {q_{\boldsymbol{s}}^k} \Bigg)\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }, \\[-15pt]\end{split}$

$ q_{{\boldsymbol{rs}}}^{ij} = \dfrac{{\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}{p_{{\boldsymbol{rs}}}}}}{{\displaystyle\sum\limits_{{\boldsymbol{rs}}} {\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}{p_{{\boldsymbol{rs}}}}} }}, $

其中配分函数$ Z_{i\rightarrow j} $与$ Z_i $可以分别表示为

$\begin{split} {Z_{i \to j}} =\;& \sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\exp \left[ {\sum\limits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\left( {1 - \frac{{{d_i}}}{n}} \right)\sum\limits_k {q_{\boldsymbol{s}}^k} } } \right]\\&\times\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }, \\[-10pt]\end{split}$

$\begin{split} {Z_i} = \;&\sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\exp \left[ {\sum\limits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\left( {1 - \frac{{{d_i}}}{n}} \right)\sum\limits_k {q_{\boldsymbol{s}}^k} } } \right]\\&\times\prod\limits_{k \in N\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}{p_{{\boldsymbol{rs}}}}} }. \\[-10pt]\end{split}$

上述迭代过程我们记为MFBP算法.
类似地, 针对IBP算法处理过程, 也可以对联合概率公式(3)先进行变形、再平均场近似处理. 即概率分布(3)式可以表示为

$\begin{split} q\left( {\boldsymbol{g}} \right) \propto\;& \prod\limits_{i < j} {{{\left( {\frac{{{p_{{g_i}{g_j}}}}}{{1 - {p_{{g_i}{g_j}}}}}} \right)}^{{A_{ij}}}}\left( {1 - {p_{{g_i}{g_j}}}} \right)} \prod\limits_i {{\gamma _{{g_i}}}} \\=\;& \left[ {\exp \left( {\sum\limits_{i < j} {\ln \left( {1 - {p_{{g_i}{g_j}}}} \right)} } \right)} \right]\\&\times\prod\limits_{i < j} {{{\left( {\frac{{{p_{{g_i}{g_j}}}}}{{1 - {p_{{g_i}{g_j}}}}}} \right)}^{{A_{ij}}}}} \prod\limits_i {{\gamma _{{g_i}}}}, \\[-15pt]\end{split}$

令${\hat h_i} = \displaystyle\sum\nolimits_{j \ne i} {\ln \left( {1 - {p_{{g_i}{g_j}}}} \right)}$, 通过平均场近似有

$\begin{split} \hat h_i^{\rm MF}{\rm{ = }}\;&\sum\limits_{j \ne i} {\sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^j\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)} }\\ \approx\;& \sum\nolimits_{\boldsymbol{s}} {\left[ {\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)\sum\nolimits_j {q_{\boldsymbol{s}}^j} } \right]}. \end{split}$

因此联合概率公式(3)经过一次平均场近似可以写成

$\begin{split} &q\left( {\boldsymbol{g}} \right) \propto \left[ {\exp \left( {\sum\limits_i {\tilde h_i^{\rm MF}} } \right)} \right]\prod\limits_{i < j} {p_{{g_i}{g_j}}^{{A_{ij}}}} \prod\limits_i {{\gamma _{{g_i}}}}\\ =\;& \prod\limits_{i < j} {{{\left( {\frac{{{p_{{g_i}{g_j}}}}}{{1 - {p_{{g_i}{g_j}}}}}} \right)}^{{A_{ij}}}}} \\&\times \prod\limits_i \Bigg\{ \exp \Bigg[ \sum\limits_s \ln \Bigg( {1 - {p_{{g_i}s}}} \sum\limits_j {q_s^j} \Bigg){\gamma _{{g_i}}} \Bigg\}. \\[-10pt]\end{split}$

把此联合概率分布定义为一个成对马尔科夫网上的概率分布, 而且马尔科夫网与原始网络一致. 其中, 对于一条边$ (v_i, v_j) $, 其边位势为$\psi_{ij}(g_i, g_j) = $$ {\dfrac{{{p_{{g_i}{g_j}}}}}{{1 - {p_{{g_i}{g_j}}}}}}$, 对于一个节点$ v_i $, 其点位势为$\phi_i(g_i) = $$ \exp \left[ {\displaystyle\sum\nolimits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{g_i}{\boldsymbol{s}}}}} \right)\displaystyle\sum\nolimits_j {q_{\boldsymbol{s}}^j} } } \right]{\gamma _{{g_i}}}$. 所以直接根据BP 算法的公式(15)—(17)式, 可以得到:

$\begin{split} \eta _{\boldsymbol{r}}^{i \to j} =\;& \frac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_{i \to j}}}}\exp \left( {\sum\limits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\sum\limits_k {q_{\boldsymbol{s}}^k} } } \right)\\&\times\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }, \end{split}$

$\begin{split} q_{\boldsymbol{r}}^i =\;& \frac{{{\gamma _{\boldsymbol{r}}}}}{{{Z_i}}}\exp \left( {\sum\limits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\sum\limits_k {q_{\boldsymbol{s}}^k} } } \right)\\&\times\prod\limits_{k \in N\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }, \end{split}$

$ q_{{\boldsymbol{rs}}}^{ij} = \dfrac{{\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}\dfrac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}}}{{\displaystyle\sum\limits_{{\boldsymbol{rs}}} {\eta _{\boldsymbol{r}}^{i \to j}\eta _{\boldsymbol{s}}^{j \to i}\dfrac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }}, $

其中配分函数$ Z_{i\rightarrow j} $与$ Z_i $可以分别表示为

$\begin{split} {Z_{i \to j}} =\;& \sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\exp \left[ {\sum\limits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\sum\limits_k {q_{\boldsymbol{s}}^k} } } \right]\\&\times\prod\limits_{k \in N\left( i \right)/j} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }, \end{split}$

$\begin{split} {Z_i} = \;&\sum\limits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}}\exp \left[ {\sum\limits_{\boldsymbol{s}} {\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right)\sum\limits_k {q_{\boldsymbol{s}}^k} } } \right]\\&\times\prod\limits_{k \in N\left( i \right)} {\sum\limits_{\boldsymbol{s}} {\eta _{\boldsymbol{s}}^{k \to i}\frac{{{p_{{\boldsymbol{rs}}}}}}{{1 - {p_{{\boldsymbol{rs}}}}}}} }. \end{split}$

上述迭代过程记为MFIBP算法.
对比BP算法与MFBP算法、IBP算法与MFIBP算法的最终迭代公式, 发现不同的只是外场量. 在节点个数足够多, 且网络比较稀疏的情况下, BP与IBP算法的外场都可以表示为

$\begin{split} \prod\limits_{k \in {\boldsymbol{V}}} {\left[ {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right]} =\;& \exp \left[ {\sum\limits_k {\ln \left( {1 - \sum\limits_{\boldsymbol{s}} {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } \right)} } \right] \\=\;& \exp \left( {\sum\limits_{\boldsymbol{s}} {\sum\limits_k {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } } \right). \end{split}$

又因为, $ {{{d_i}}}/{n} \to 0 $且$\ln \left( {1 - {p_{{\boldsymbol{rs}}}}} \right) \to $$ {p_{{\boldsymbol{rs}}}}$, 则MFBP算法的外场量可以近似为

$\begin{split} &\exp \left[ \sum\limits_{\boldsymbol{s}} \ln \Big( {1 - {p_{{\boldsymbol{rs}}}}} \Big)\Big( 1 - \dfrac{{{d_i}}}{n} \Big)\sum\limits_k {q_{\boldsymbol{s}}^k} \right] \\&\to \exp \left( {\sum\limits_{\boldsymbol{s}} {\sum\limits_k {q_{\boldsymbol{s}}^k{p_{{\boldsymbol{rs}}}}} } } \right), \end{split}$

MFIBP算法的外场量可以近似为

$ \exp \Big[ \sum\limits_{\boldsymbol{s}} \ln \left(\right. {1 - {p_{{\boldsymbol{rs}}}}} \left.\right)\sum\limits_k {q_{\boldsymbol{s}}^k} \Big] \to \exp \left( {\sum\limits_{\boldsymbol{s}} {\sum\limits_k {q_{\boldsymbol{s}}^k p_{{\boldsymbol{rs}}}} } } \right). $

可以看出: 理论上, 经过平均场近似后得到的MFBP算法和MFIBP算法在一定条件下分别与原有的BP算法和IBP算法是等价的. 但是在理论推导上, 平均场近似的方法要更简单合理且便于理解. MFIBP算法与MFBP算法相比较, 少了一次平均场近似.
下面将对比BP算法和MFBP算法, 以及IBP算法和MFIBP算法实验上的效果.

下面将在社团结构及核边结构探测上进行实验验证. 在社团结构探测中, 为了验证不同算法探测结果的好坏, 将采用NMI指标进行衡量, 即^[39]:

$ {{\rm NMI}\left( {X,Y} \right) = \frac{{2I\left( {X,Y} \right)}}{{H\left( X \right) + H\left( Y \right)}}}, $

其中, X与Y分别为算法探测的结果以及真实的结果; $ I(X, Y) $表示X与Y之间的互信息; $ H(X) $与$ H(Y) $分别表示X与Y的信息熵.
对于核边结构的探测结果的评价, NMI并不是一个合理的指标, 这是因为核边结构的核与边的地位并不是等价的. 所以将采用$ F1 $指标进行衡量, 即^[40]:

$ {F1 = {{2PR}}/({{P + R}}}), $

其中, P表示预测结果中, 预测为正(本文中为核)的样本中预测正确的概率, R表示数据样本中, 正样本中预测正确的概率.
2

-->

4.1.社团结构探测

首先采用随机块模型生成具有社团结构的基准网络, 以此来验证MFBP算法和MFIBP算法分别与BP算法和IBP算法的区别. 假设有k个模块, 则每个节点就以$ \gamma_{\boldsymbol{r}} $的概率分配到第r模块中. 社团之间节点的连边概率定义为$ p_{{\boldsymbol{rs}}} = c_{\rm in}/n ({\boldsymbol{r}} = {\boldsymbol{s}}) $和$p_{{\boldsymbol{rs}}} = c_{\rm out}/n ({\boldsymbol{r}} \ne {\boldsymbol{s}})$, 因此网络的平均度为$c = $$ \displaystyle\sum\nolimits_{{\boldsymbol{r}} < {\boldsymbol{s}}} {2{\gamma _{\boldsymbol{r}}}{\gamma _{\boldsymbol{s}}}{c_{\rm out}}} + \displaystyle\sum\nolimits_{\boldsymbol{r}} {\gamma _{\boldsymbol{r}}^2{c_{\rm in}}}$.
如果每个模块的大小一致, 即$ \gamma_{\boldsymbol{r}} = 1/k $, 则网络的平均度为$ c = \left[ {{c_{\rm in}} + q\left( {{c_{\rm out}} - 1} \right)} \right]/q $. 定义模块外部与社团内部的连接概率比为$ \varepsilon = {c_{\rm out}}/{c_{\rm in}} $. 所以, 当$ \varepsilon < 1 $时生成网络具有社团结构, 且越小社团结构越明显; 当$ \varepsilon = 1 $时生成网络为ER随机图. 当$ n \to \infty $时, 网络检测社团结构会存在一个阈值^[38]:

$ {\varepsilon ^ * } = ({{\sqrt c - 1}})/({{\sqrt c - 1 + k}}), $

即, 如果$ \varepsilon > \varepsilon ^ * $时, 各种社团算法都无法检测出社团结构; 当$ \varepsilon < \varepsilon ^ * $时, BP算法可以探测出社团结构.
从理论上, 前面的分析可以看出, 当网络节点个数足够多, 网络连接比较稀疏时, BP算法、MFBP算法、IBP算法以及MFIBP算法是近似等价的, 通过图1(a)可以验证, 所有结果为10 次结果的平均值(下同). 通过图1(b)—图1(d)可以发现, 当网络连接稀疏且规模较小时, BP算法、MFBP算法、IBP算法以及MFIBP算法的结果也是一致的. 但是, 在网络连接稠密的情况下(如图2), 可以发现MFBP算法要优于BP算法, 且与MFIBP算法、IBP算法一致.
图 1 稀疏基准网络社团结构探测　(a) $n = 10000$, $c = 3$, $\gamma = {[0.5,\; 0.5]^{\rm T}}$, $\varepsilon^* = 0.27$; (b) $n = 200$, $c = 10$, $\gamma = {[0.5,\; 0.5]^{\rm T}}$, $\varepsilon^* = 0.520$; (c) $n = 200$, $c = 10$, $\gamma = {[0.3,\; 0.7]^{\rm T}}$; (d) $n = 200$, $c = 10$, $\gamma = {[1/3,\; 1/3,\; 1/3]^{\rm T}}$, $\varepsilon^* = 0.419$
Figure1. Community structure detection in sparse benchmark networks: (a) $n = 10000$, $c = 3$, $\gamma = {[0.5,\; 0.5]^{\rm T}}$, $\varepsilon^* = 0.27$; (b) $n = 200$, $c = 10$, $\gamma = {[0.5,\; 0.5]^{\rm T}}$, $\varepsilon^* = 0.520$; (c) $n = 200$, $c = 10$, $\gamma = {[0.3,\; 0.7]^{\rm T}}$; (d) $n = 200$, $c = 10$, $\gamma = {[1/3,\; 1/3, \;1/3]^{\rm T}},$ $\varepsilon^* = 0.419$

图 2 稠密基准网络社团结构探测　(a) $n = 200$, $c = 50$, $\gamma = {[0.5,\; 0.5]^{\rm T}}$, $\varepsilon^* = 0.75$; (b) $n = 200$, $c = 50$, $\gamma = {[0.3, \;0.7]^{\rm T}}$, (c) $n = 200$, $c = 30$, $\gamma = {[1/3,\; 1/3,\; 1/3]^{\rm T}}$, $\varepsilon^* = 0.599$
Figure2. Community structure detection in dense benchmark networks: (a) $n = 200$, $c = 50$, $\gamma = {[0.5,\; 0.5]^{\rm T}}$, $\varepsilon^* = 0.75$; (b) $n = 200$, $c = 50$, $\gamma = {[0.3,\; 0.7]^{\rm T}}$; (c) $n = 200$, $c = 30$, $\gamma = {[1/3,\; 1/3,\; 1/3]^{\rm T}}$, $\varepsilon^* = 0.599$

2 -->

4.2.核边结构探测

生成具有核边结构特性的基准网络来验证MFBP算法和MFIBP算法分别与BP算法和IBP算法的区别. 给定3个参数: 核与核的连边概率$ P_{\rm CC} $, 核与边的连边概率$ P_{\rm CP} $, 以及边与边的连边概率$ P_{\rm PP} $, 当设置$ {P_{\rm CC}} > {P_{\rm CP}} > {P_{\rm PP}} $时, 可以生成一个具有核边结构的网络^[28]. 在本文中, 设$ P_{\rm PP} = $$ 0.05 $, $ P_{\rm CC} = \theta $, $ P_{\rm CP} = 0.6\theta $. 网络的大小$ n = 200 $, 其中核节点的个数为50, 边节点的个数为150. 图3所示为4种算法在不同参数基准网络的核边结构的探测情况, 其中实验结果为10次结果的平均值. 同样可以发现, MFBP算法要优于BP算法, 且与MFIBP算法、IBP算法一致. 在BP 算法中, 当核边结构足够明显时, 却表现出了较差的结果, 这是因为核边结构明显, 意味着网络结构稠密, 特别是核节点的度很大, 这会导致(19)式关于用所有节点集合近似邻居的补集变得极其不合理, 而MFBP算法、IBP算法以及MFIBP算法没有这方面的缺陷.
图 3 基准网络的核边结构探测
Figure3. CP structure detection in benchmark networks.

接着, 在美国航空网络上验证本文中的算法, 其中每个节点代表一个机场, 共有332个节点; 每条边表示两个机场具有航班, 共有2126 条边^[5]. 应用BP 算法, 求解结果有27 个核; 应用MFBP算法, 求解结果有34个核; 应用IBP算法, 求解结果有47 个核; 应用MFIBP算法, 求解结果有47个核. 其中对应的参数结果见(54)式—(57)式. 为了验证上述4个算法求得结果的精确性, 采用MCMC采样求解EM算法, 共采样1000 轮(每轮n次), 得到的结果共有47 核, 具体参数如(58)式. 如果把MCMC结果作为真实结果, 则BP算法、MFBP算法、IBP算法以及MFIBP算法结果的$ F1 $值分别为0.730, 0.840, 1.000以及1.000. 从上述结果可以看出, 在美国航空网络的核边结构实验中, MFBP算法结果要优于BP算法, MFIBP算法结果与IBP算法一致, 且接近精确解. 而且MFIBP算法要优于MFBP算法, 这是因为MFBP算法比MFIBP算法多一次平均场近似.

$ {\gamma _{\rm BP}} = \left[ \begin{array}{l} 0.081\\ 0.919 \end{array} \right],\;{p_{\rm BP}} = \left[ {\begin{array}{*{20}{c}} {0.873}&{0.151}\\ {0.151}&{0.012} \end{array}} \right], $

$ {\gamma _{\rm MFBP}} = \left[ \begin{array}{l} 0.102\\ 0.898 \end{array} \right],\;{p_{\rm MFBP}} = \left[ {\begin{array}{*{20}{c}} {0.836}&{0.120}\\ {0.120}&{0.010} \end{array}} \right], $

${\gamma _{\rm IBP}} = \left[ \begin{array}{l} 0.142\\ 0.858 \end{array} \right],\;{p_{\rm IBP}} = \left[ {\begin{array}{*{20}{c}} {0.711}&{0.074}\\ {0.074}&{0.008} \end{array}} \right], $

${\gamma _{\rm MFIBP}} = \left[ \begin{array}{l} 0.142\\ 0.858 \end{array} \right],\;{p_{\rm MFIBP}} = \left[ {\begin{array}{*{20}{c}} {0.713}&{0.074}\\ {0.074}&{0.008} \end{array}} \right], $

${\gamma _{\rm MCMC}} = \left[ \begin{array}{l} 0.142\\ 0.858 \end{array} \right],\;{p_{\rm MCMC}} = \left[ {\begin{array}{*{20}{c}} {0.709}&{0.074}\\ {0.074}&{0.008} \end{array}} \right]. $

本文详细描述了随机块模型探测社团结构、核边结构的理论, 以及MCMC和BP算法求解该理论的方法, 并分析其优缺点. 通过理论分析以得到: 1)原始的BP算法不适用于核边结构的探测; 2)修正后得到的IBP算法, 虽然适用于核边结构, 但在理论上推导近似多且不明朗; 3)无论是BP算法还是IBP算法, 所处理的马尔科夫网(联合概率分布), 既不是原始网络, 也不满足网络稀疏条件. 最后一条是造成随机块模型中BP算法和IBP算法理论难以理解的主要原因. 针对于此, 本文首先对联合概率分布进行平均场近似, 得到一个好的联合分布函数, 然后可以直接套用BP 算法理论, 且在这一步不需要任何近似, 使得BP算法求解随机块模型这一套方法, 不仅在结果上效果显著, 而且在理论上通俗易懂. 这种先通过平均场近似再套用BP 算法还适用于利用最小化哈密顿量(关于模块度的函数)探测社团结构的方法中^[38], 和本文一样, 将使得理该论推导部分变得简单易懂.
本文的计算工作得到了安徽大学高性能计算平台的支持.